WRIA-1 ETc Estimation

The Development Of A GIS Based Procedure

Leon Basdekas

GIS in Water Resources




Approach Decision

Exporting Data Procedures

Calculating ETo

Combining Land Themes

Assigning Kc values

Calculating ETc







The Utah Water Research Water Lab (UWRL) at Utah State University (USU) has been contracted to provide technical assistance in the development of a management plan for Water Resource Inventory Area 1(WRIA-1) in far northwestern Washington Sate.   The Washington State Watershed Management Act allows citizens and local governments to develop a plan specifically related to local water-related circumstances and needs.  It requires all participating local governments to address water quantity and optionally water quality, instream flows and fish habitat all of which are being addressed in WRIA-1.   A memorandum of understanding was signed between Lummi Nation, Nooksack Tribe, City of Bellingham, Public Utility District 1 (PUD 1), and Whatcom County to participate in the effort.


In order to estimate the consumptive use within the study area an estimation of evapotranspiration (ETc) is needed.  This estimate will comprise a portion of a preliminary water balance.  The intent of this term project is to develop a procedure for estimating the ETc within the study area.  This process will later be included in a decision support system (DSS), which will be used to model management alternatives within the WRIA.   Hargreaves equation was selected based on review of available climatic data and recommendations in FAO-56 when temperature is the only climatic data available.  Future efforts may make use of the Penman-Monteith equation while estimating the missing climatic data as outlined in FAO-56.  This decision will be made after examining the Penman-Monteith equations sensitivity to those missing data given the terrain and climate diversity within the study area.  The increased level of effort to obtain a more accurate estimate of Etc may not yield much true increased accuracy due to the physical variability of terrain, lack of representative crop coefficients (Kc) for much of the upland areas and other factors limiting accurate estimations.




The following data sets are being utilized:


Drainage boundaries for WRIA-1, PUD 1

PRISM temperatures mean monthly maximum and minimum for 1961-1990, Oregon Climate Services

Land cover state wide 1991 GAP, Washington State Department of Fish and Wildlife


Figure 1

 Study Area Location (Bar scale applies to drainage detail)


Grid verses Polygon


Two approaches to this analysis were considered.  The first was a grid-based analysis that had the advantage of being a familiar process of using the map calculator for grid-by-grid calculations.  The second process was a polygon approach, which was not nearly as familiar and the particulars of how to conduct the analyses, was not known. 

Despite the unfamiliarity with polygon analysis, it has several advantages for this project.  Polygons have more representative boundaries instead of a grid approximation of a polygon boundary, the greatly reduced size of polygon data files versus grid data files and more versatility with respect to attributing the theme at least as far as known methods are concerned. 


The polygon method was attempted after the discovery of the field calculator.  The field calculator was discovered during the resolution of a problem related to a string versus number data format, the particulars of which will be discussed later.  It was upon this discovery that the polygon method appeared feasible.  Procedures and analyses were conducted with ArcView 3.2 and Excel 2000.


Exporting GIS data:


An avenue script was obtained and used to determine the center of mass (centroid) of each polygon to obtain the northing and easting of the centroid point.  This script also generated a unique identifier for each centroid point.  This unique drainage id was used later for joining processes described later.


Drainage names were assigned to the centroid points.   A spatial join was used and was successful in all but about a half dozen cases.  In the failed instances some polygons had two centroids assigned due to the shape of the polygon being such that the centroid was actually outside its physical boundary (Figure 2).  These were manually edited taking advantage of the summary operations available in ArcView.  Name and count were summarized from the joined tables.  The resulting summary table was linked back to the feature attribute table (FAT) so dynamic selection could take place by selecting those Drainages with more than one centroid.


A second avenue script was used to export the coordinates of the centroid points theme to an Excel file.  This data was saved as an ASCII text file for use in Corpscon (U.S. Army Corps of Engineers coordinate conversion software).  Corpscon was used to convert from UTM NAD 83 meters to NAD 83 geographic coordinates in decimal degrees.  This procedure was required to obtain the latitude for each centroid for later use in the ETo calculations.  The latitudes were later imported into an Excel spreadsheet.


It should be noted that a USGS HUC type system was developed for this theme however its ability to create a unique identifier had not been fully tested by (PUD 1) therefore it was not used for this preliminary analysis.






















Figure 2

Map showing how more than one centroid was incorrectly applied to one drainage.


PRISM data grids (these crossed the border into Canada and need to be considered provisional) were clipped using the SINMAP extension and the rectangular grid cutter.  This was done to reduce the size of the grid theme to a more manageable size.  The grid was then summarized by drainage using the drainage id field joined from the centroid theme.  These processes generated summary tables that were then exported to Excel by using a slightly customized avenue script.  The script was altered to write selected fields to an existing Excel sheet and not a new workbook as originally designed.  The export process was repeated for all 24-temperature grid themes.  An Excel macro was created to take the summary temperature data and populate the calculation spreadsheet.   It should be noted that the temperature grids were integer grids and the values were degrees centigrade times 100.  This correcting adjustment was made in Excel in addition to the other necessary temperature calculations required by the equations in the following section. 


Calculating ETo with an Excel spreadsheet:                                                 Top


Excel was used to conduct the remaining calculations necessary for ETo.  The formulas used in Excel are as follows:


Hargreaves ETo equation where:


ETo = 0.0023(Tmean + 17.8)(Tmax - Tmin)^0.5*Ra


Temperature based equation with solar radiation

The extraterrestrial radiation, Ra, for each day of the year and for different latitudes can be estimated from the solar declination and the time of the year by:

Ra=15.392*dr[ws sin(f)sin(d)+cos(f)sin(d)sin(ws)]

Where, D-denotes variation by drainage, M-denotes variation by month, B-Both

Ra extraterrestrial radiation [mm day-1],                                  B
dr inverse relative distance Earth-Sun                                     M
w s sunset hour angle (Equation 25) [rad],                              D
f latitude [rad],                                                                       D
d solar decimation (Equation 24) [rad].                                  M

The relative distance Earth-Sun, dr, and the solar declination, d, are given by:

The sunset hour angle, w s, is given by:

w s = arccos [-tan (f) tan (d)]

where, J is the Julian day or the number of the day in the year between 1 (1 January) and 365 or 366 (31 December).

Values for J along with the seasonal effects in d and dr are shown in Table 1 as entered and calculated in Excel.





























Julian Day













Solar Dec d













Rel Dist dr














Table 1


Table 2 is a sample of ETo data after all the calculations have been performed.  ETo has been converted to inches/month.  ID is the field corresponding to a specific drainage.




Et Jan

Et Feb

Et Mar

Et Apr

Et May

Et June

Et July

Et Aug

Et Sept

Et Oct

Et Nov

Et Dec








































































Table 2



After the somewhat complex process of calculating 181 drainages by 12-months, the ETo summary sheet was saved in dbf format to be joined to the drainage theme by the drainage id field.  The calculation spreadsheet can be examined by clicking on the following link WRIA1_ET.xls.


A new drainage polygon theme was created with just the unique id and drainage name.  Using the unique id the new drainage theme was joined with the newly created dbf file containing the ETo data.  Now ETo attributes could be displayed as polygons by drainage.  The attributes of this new theme are the 12 monthly ETo values.  In a similar manner mean temperature data were exported from Excel and then joined to the centroid point theme.  This was done so temperature and ETo could be viewed simultaneously with ETo data.  August ETo (inches) and mean temperatures (oC) are shown in Figure 3.

Figure 3

August ETo (inches) and mean temperatures (oC)


Combining land themes                                                                                  Top


For the developed areas a modified tax parcel theme was used.  This data set is used to find location of agriculture lands at a finer resolution than possible with the GAP data.  For the remaining areas within the WRIA GAP data was used.  An example can be seen in Figure 4 where the tax parcels identified as agriculture are shown over different land cover types in the GAP data.



Figure 4

Differences in land themes (GAP legend too complex for display) light green represents GAP identified agricultural lands other colors are other land covers (Agriculture parcels already removed by the process described next)


The agricultural lands needed to be separated from the rest of the tax parcel theme (total size ~124,000 records).  The ArcView query builder would not function with this large of a data set so an alternative selection process was used.  A summary table was created which contained values of land use codes and count.  This table was linked both ways with the original FAT.  The codes corresponding to agricultural lands were selected from the summary table and by link from the FAT.  These selected records were subsequently saved as a new theme (Figures 5 and 6).









Figure 5

 Map of Selected Agriculture parcels





Figure 6

  Linked tables and Example Kc table waiting to be joined


Assigning Kc values


For consumptive use calculations Kc values need to be applied to the ETo data to determine ETc by use of the following equation:


 ETc = Kc*ETo 


The Kc values need to be added to the land data sets.   Two new tables were created which contained Kc values for each land code (one table for each theme).  A problem appeared during this process for the newly created tax_ ag theme.  The FAT use_category field (the join field) was a string and could not be converted using Excel formatting commands.  As a remedy a new field was created and the field calculator was used to add the use category code as the same value but as a number and not a string format. (It was this use of the field calculator which provided insight into a polygon based analysis.   All the work done to this point would have been required for the grid analysis as well.)  The join processes were now performed to get the Kc values into the two land theme FATís.


The Geoprocessing wizard was used to create a union between the tax_ ag theme and the GAP data.  The resulting FAT contained two Kc fields one for each of the previous themes. When a Kc value was present for one field the other Kc field was zero so each polygon (record) had only one valid Kc. 


A second union was performed between the new composite land theme and the ETo theme.  Since the composite land use/cover theme did not extend across the Canadian border, there were polygons that did not have Kc values.  Those polygons were assigned a default Kc value of 1.00 by creating a new field and adding the value where both Kc values were absent for a particular record. 


A new field title, comp_ Kc, was created.  The field calculator was used to add the two Kc fields and the default Kc field which created a unique value for each polygon. 


Calculating ETc with ArcView


A drainage scale weighted average ETc needed to be determined.  A new field was added for the polygon area and the area calculated with the field calculation tool.  Total drainage area was added to a new field by use of join from a previously created summary table.  This now gave each record an attribute of the drainage area in which the polygon resides.    The next step was to multiply the polygon area and the comp_ Kc field and then divide that product by the total drainage area and then multiply by the drainage ETo (Part_cons=Polygon area*comp_Kc/drainage area*ETo).  This resulting value is what each polygon within a particular drainage contributes to ETc.  When a summarize operation was conducted on the modified FAT for the last union operation the name, count and ETc contributions were selected.  This resulted in a composite ETc value for every drainage.  A script has been written to automate this calculation process, which will streamline the process for the remaining eleven months.  See Figure 7 for ETc in August.





Figure 7

  August ETc in inches




This project is divided into five basic components.


1.               Summarizing and then exporting data from ArcView into Excel

2.               Calculating ETo within Excel and exporting summary tables in dbf format

3.               Combining land themes into a composite theme

4.               Assigning Kc values to various land use and land cover themes by using relational tables

5.               Calculating Etc within ArcView


This project of developing a process to calculate Etc by using GIS was a success.  By developing this process the future integration into the DSS will be smoother.  It will also aide in the development of a much needed preliminary water balance.  Any new data that arrives can now be incorporated fairly easily through relational tables and then running the script to perform the process described in Calculating ETc with ArcView .  Even if basic information should change such as the provisional PRISM temperature data the scripts and macros used in this process can be used to greatly reduce the time required to form a new product.   




Acknowledgments                                                                                           Top


Dr.  Tarboton-  Project definition and subtle motivation

Dr.  Pack-  Data reprojections (part of entire WRIA study effort)

Mark Winkelaar-  USU GIS specialist provided some scripts and tricks of the trade