# Thiessen & Isohyetal Precip

1.) THIESSEN POLYGONS METHOD FOR PRECIPITATION

Fig 1. Thiessen polygons constructed from a set of weather station points (precipitation measurements). Each point gets it own polygon and the point value is distributed throughout the entire polygon. So the precip value, initially measured at the station point, expands to become the value for the entire polygon.

Mapping point data usually involves some sort of interpolation. There are many point-interpolation methods out there (contours, IDW, TIN, tension spline, kriging, trend surface, etc.) as well as non-interpolation methods for displaying point values (proportionally-sized circles). The Thiessen Polygon method is an interpolation method commonly used for precipitation, but can be used on other point datasets.

Thiessen Polygons are Voronoi Cells, a geometric means of dividing up an area given a set of known values at a relatively small number of points. This interpolation method was first applied to weather station data by A.H. Thiessen (1872-1956), an American meteorologist for the Weather Bureau, now called NOAA National Weather Service (NWS). The Thiessen polygon method is one of 5 different ways station precip measurements are extrapolated. The others are Arithmetic Mean, Isohyetal, Distance-Weighted Grid, and MapX grid.

By the way, “Theissen” is pronounced THEE-sun (surname has Norman roots).

LESSON
The goal in this lesson is to find the area-weighted precipitation for a set of Watersheds given a set of point measurements (precipitation measurements at a set of stations).

```Instructors: If you are setting up your own lesson, you will need a.) a shapefile containing a set of numbered station points with precip values and b.) a shapefile of named watershed polygons. Stations should be distributed both inside and outside watershed polygons (see Fig. 2). There should be more watersheds than station points. Provide only one precip measurement per station to keep things simple. An eexaggerated range of precip values makes it easier to find errors in calculations. Students should report the area-weighted precip for each Watershed. I usually have students report values at 3 scales: Watershed, Watershed Zone (sets of adjacent watersheds; add a ZONE field to the attribute table -> A,B,C,D), and for the entire Basin (all watersheds). Discuss results with respect to scale and measurement density. ```

– Start with Watershed polygons and Weather_Station points (shapefiles). Each station should have a unique precip value in the attribute table.

– Add labels to Weather_Station points using the Precip field.

File Naming Terminology
Watershed polygons = WS_polys
Thiessen polygons = T_polys
Intersected polygons = Int_polys

Fig 2. Watersheds A,B,C,D in blue. Weather_Stations with precipitation values in red. In this example, every watershed has one station and one precip value. Click for larger image.

– Construct Thiessen polys from the station points layer (Analysis tools > Proximity > Create Thiessen polygons). Once created, the station precip point value is now distributed everywhere throughout the T_poly.

– Calculate the area of each of the WS_polys. Do this in the attribute table. Create a New Field called WSArea_km (FLOAT format) and Calculate Geometry in km2 units.

** Note: Make sure your coordinate system and display units are correct (Data Frame Properties > General tab and Coordinate System tab).

Fig 3. Thiessen polygons are constructed from the station points (not the precip values). The map area is divided into polygons based on Voronoi diagram rules. Each precip value is now  distributed throughout the entire Thiessen polygon. Click for larger image.

– Intersect WS_polys with T_polys (Analysis tools > Overlay > Intersect). This gives you Int_polys.

– Add a New Field to Int_polys layer called IntPolyID (Short Integer). Use Field Calculator to duplicate FID in this new field. See Fig 4. The reason for doing this simply for visualization. We want to symbolize by Int_poly ID number. But because FID field (or ObjectID field) is unavailable in the Symbology > Categories choices, we have to duplicate the values in a new field so we can use it to color each intersect poly uniquely (next step).

– Symbolize Int_polys on this new field (IntPolyID) by Categories > Unique Values, so that all intersected polygons have their own color.

Fig 4. The intersection of Watershed polys and Thiessen polys creates many new polygons (ID numbers 1 through 23). The next step is to calculate the area of each intersect polygon. Recognize that the set of polygons generated by the Intersect tool may be numbered in such a way that one ID# may represent more than one entity (for example see polygons #8, #17, #19, #20). Intersecting complicates the attribute table. Click for larger image.

– Calculate areas of each Int_poly layer in the attribute table just like you did for the watersheds layer (New Field called IntArea_km, FLOAT, Calculate Geometry, km2). OK any warning pop ups.

– Examine the attribute table for Int_polys:
FID = ID for Int polygons.
Precip_1 = Precip value for the intersect polygons.
WSArea_km = Area of each Watershed.
IntArea_km = Area of each intersect polygon.

– Create a new field called PropArea (FLOAT) and use Field Calculator to find the proportion of the watershed each intersect poly represents.
Formula example: [IntArea_km] / [WSArea_km]

– Create another new field called Weighted (FLOAT) and use Field Calculator to find the proportion of the watershed each intersect poly represents. This gives you the weighted precip for each Int_poly.
Formula example: ([IntArea_km] / [WSArea_km]) * [Precip_1)

– If you want to pretty up the table for printing, export the Int_polys attrribute table (select all/copy/paste or export as .txt or .dbf) for use in Excel. Do not add the new table to the map; open it in Excel and work with it there.

** Note: How to Open .txt files Exported from ArcGIS in Excel 8:
File > Open > navigate to the file –> Make sure you change file type to show “All Files”
Original Data Type = Delimited
Delimiters = Comma
Column Data Format = General

– Since several Int_polys comprise each Watershed, sum the weighted precip values generated in the previous step for each Watershed. This is easily done in Excel. This is your weighted average precip by Watershed via the Thiessen method.

– If you need to report weighted precip by Zone (groups of Watersheds), you will need the area for each Zone polygon. Substitute Zone_Area for WSArea_km in the proportion calculation (PropArea) and weight calculation (Weighted).

– If you need to report for entire Basin (all Watersheds), you will need total area for all Watersheds (or Zones). Substitute BArea_km for WSArea_km in the proportion calculation (PropArea) and weight calculation (Weighted).

Fig 5. Example Excel spreadsheet. This table shows a record for each Intersect polygon, which is overkill. An improvement would be to show only records for Watersheds, then merge cells in additional columns for Zones and Basin.

Refs for Thiessen Method:

NWS Overview: Arithmetic Mean, Thiessen, Isohyetal & MapX methods LINK
A Review of Spatial Interpolation Methods for Environmental Scientists PDF
Mair & Fares (2010) PDF
Thiessen, A.H. (1911) Precip averages for large areas. Monthly Weather Review 39 LINK
Weather Bureau historical note on Major Thiessen LINK

2.) ISOHYETAL METHOD FOR PRECIPITATION

Iso = single, as in a single value, the value along a contour line
Hyetal = pertaining to rain and rainfall distribution

Isohyetal maps are just contour maps of precipitation. Each isohyet (each contour line) represents equal precip along its length. Often one map will represent a snapshot of precip in time. In this example, we find the average precip for a watershed at one snapshot in time given station measurements. You could substitute your study area boundary (or limit of contouring) for the watershed boundary.

Start with precip measurements (rainfall) collected at several points (weather stations). There are several stations in the example watershed.

Contour the precip values using standard contouring rules. ArcGIS can create contours from point data via IDW, TIN, Interpolate to Raster, 3D Analyst. In a pinch, you could do it by hand. Google it. In the example, we use a contour interval of 0.5.

We need to determine the area represented between each isohyet. To do this, draw polygons. The blue polygon, for example, represents area between the 4.0 and 4.5 isohyet contours. Right-click > Properties to locate each graphic’s Area info. Make sure your Data Frame is set to the proper coordinate system and map display units. Click for larger image.

Area in square miles (the label on this figure is incorrect) for each of the inter-contour polygons shown in blue text.

Next, determine the Average Precip value for each isohyet zone, as shown by the dashed lines and black text. The values are halfway between the isohyet contour values. Don’t map them. Rather, enter the values into Column B in table below.

Set up a table like this. Include columns for Isohyet Zone, Average Precip Contour, Polygon Area (zone area), and 3 columns for calculations. Formulas are shown for calculations in red italic text. Click for larger image.

Refs for Isohyetal Method:
NWS Overview: Arithmetic Mean, Thiessen, Isohyetal & MapX methods LINK