Estimating Evapotranspiration using Landsat 5 TM and ArcGIS 10.1

 

Prepared by Ayse Kilic

Modified by David Tarboton

GIS in Water Resources, Fall 2012

Due date: November 6, 2012

 

Purpose

The purpose of this exercise is to utilize raw Landsat 5 TM data in order to display a false color composite (FCC) of the Landsat bands from a scene (path 29, Row 32), and calculate vegetation index (NDVI) and estimate evapotranspiration (ET) in the ArcGIS environment. The FCC will be used to explore various features of the study area.  The Normalized Difference Vegetation index (NDVI) will be calculated from the Landsat 5 TM bands in order to compute the fraction of vegetation cover.  We will also create a gridded reference evapotranspiration map using data from automatic weather stations. The gridded reference evapotranspiration and fraction of vegetation cover will be used to estimate actual evapotranspiration on a regional scale. Zonal statistics of NDVI and evapotranspiration will be studied.

 

It should be noted that the NDVI tool of ArcGIS uses the ‘digital numbers’ of bands 3 and 4 as a ‘default.’  We will do this same thing in this exercise to save time.  However, the more accurate way to calculate NDVI is to use the ‘reflectances’ of bands 3 and 4.  Reflectances are defined as the fraction of incoming radiation that is reflected back to the surface.  ArcGIS 10.1 now includes a new tool called “apparent reflectance.”  The apparent reflectance tool can convert the digital numbers of Landsat bands into the apparent surface reflectances.  Those reflectances can then be used in the NDVI computation to produce a more consistent value for NDVI that does not change with time of year (and sun angle).  You do not have to use these reflectances in your calculation of NDVI, but you should consider using surface reflectance in the future to compute all vegetation indices, when serious estimates are needed. Instructions on using the Apparent Reflectance function in ArcGIS 10.1 are at:

http://resources.arcgis.com/en/help/main/10.1/index.html#/Apparent_Reflectance_function/009t0000023s000000/

 

A manual process for converting digital number to surface reflectance in ArcGIS is included at the end of this exercise as an addendum.  However, again, we will only use digital number in this exercise to reduce time requirements.

Computer and Data Requirements

To carry out this exercise, you need to have a computer, which runs the ArcInfo version of ArcGIS 10.1 The data are provided in the accompanying Ex5 zip file.

 

The  following data will be used for this exercise:

1.      Landsat data, LT50290322005251EDC00. The data can be freely downloaded from http://edcsns17.cr.usgs.gov/EarthExplorer/ or http://glovis.usgs.gov/.

2.      The daily reference evapotranspiration data from 12 Automated Weather Data Network (AWDN) stations within and surrounding the study area has been downloaded from the High Plains Regional Climate Center (HPRCC) in Nebraska (http://www.hprcc.unl.edu/). The locations of the weather station, and data are prepared to create a point shape file.

3.      A detailed landuse data for Nebraska has been created at CALMIT at the University of Nebraska. The landuse data for the entire state is available at http://calmit.unl.edu/2005landuse/statewide.php

4.      Base data:  The HUC 12 level hydrologic watersheds, NHD flowlines, and Natural resources districts (NRD) boundaries were obtained from the Nebraska Department of Natural Resources (http://dnr.ne.gov/databank/spat.html). There are 23 NRDs in Nebraska with leadership responsibilities for protecting ground water from overuse and pollution.

 

Downloading Landsat data

The USGS provides Landsat 5 TM data in geoTiff (TIFF) format. As the Landsat satellite passes over Earth, the area can be identified by the path and row combination. Our study area is under the path 29 row 32. This is the southeastern part of Nebraska and Landsat passes over this area at approximately 10:15 CST. Each Landsat scene (date) has a unique identifier. Our dataset has identifier LT50290322005251EDC00. This Landsat data was acquired on 2005/9/8.

 

Landsat 5 TM data is acquired every 16 days. Landsat scene might contain some cloud cover and data for that part can be erroneous. Our dataset is a Landsat scene with 0% cloud cover. All the aforementioned information regarding Landsat data can be found in the Header or metadata file of Landsat.  These files will end with _GCP and _MTL. These files also contain the information regarding each band pertinent to that scene. Useful information like center and corner coordinates of the Landsat scene is provided in the header file. A screen shot of attributes for this Landsat scene is provided below.

 

1. VIEWING BASE MAP DATA

 

We’ll begin by getting the input data for our study area in Nebraska. Open ArcGIS desktop, and add the following layers from base data.gdb

 

 

Arrange the layer order and change the symbology of your data layers as you wish (see example below). Right click on NRD_Districts, and label features (select NRD_name).  Save your project as Ex5.

 

To turn in:

1.      Prepare an ArcMap layout showing labeled Natural Resources Districts (NRD) within NE. You should include NHD_flowlines, and HUC-12 units. Use a layout view, include a scale bar to indicate distance, a north arrow to indicate direction and labels or legends with units. A snapshot of ‘data view’ of the map is depicted above.

 

2. CREATING COLOR COMPOSITES FROM LANDSAT 5 TM GRIDS IN TIFF FORMAT

 

Landsat 5 TM has 7 different bands. These bands are useful for extracting various information related to vegetation, temperature, clouds, soil moisture, biomass, rocks, minerals and so on.  Below are details of band designations for Landsat 5 TM. More detailed information can be found at http://egsc.usgs.gov/isb/pubs/factsheets/fs02303.html.

 

 

Spectral Bands1

Wavelength  (micrometers)

Potential Information Content

Resolution (meters)

Band 1

Blue

0.45 - 0.52

Discriminates soil and rock surfaces from vegetation. Provides increased penetration of water bodies

30

Band 2

Green

0.52 - 0.60

Useful for assessing plant vigor

30

Band 3

Red

0.63 - 0.69

Discriminates vegetation slopes

30

Band 4

Near IR

0.76 - 0.90

Biomass content and shorelines

30

Band 5

Mid  IR

1.55 - 1.75

Discriminates moisture content of soil and vegetation; penetrates thin clouds.

30

Band 6

Thermal IR

10.40 - 12.50

thermal mapping and estimated soil moisture

120

Band 7

Mid  IR

2.08 - 2.35

Mapping hydrothermally altered rocks associated with mineral deposits

30

1: IR stands for infrared

 

Create a folder “imagery” under your Ex5 folder to save imagery results.  To do this open Catalog and right click the Home Ex5 location (this appears when you have saved the map document file) and select new folder.  Change the name to Imagery.

 

 

Color composites are used to facilitate viewing of imagery.  We will create a color composite from the Landsat band geoTIFF files.  Open the composite bands tool located under data management\raster\raster processing\composite bands. (You can also search for composite bands under the search tool).

 

 

 

For Input Rasters select all the bands from LT50290322005251EDC00 which contains TIFF images for the seven landsat bands. The trick is to select the first image L5029032_03220050908_B10.TIF, hold down the shift key and select the others.

 

Make sure all bands are added in the order (from L5029032_03220050908_B10.TIF to L5029032_03220050908_B70.TIF).  You can use the up and down arrows to reorder if necessary.

Save the output file as Composite under Imagery folder.

 

 

Zooming in, the resulting image looks like this.

 

 

We are going to recolor the color composite by using spectral bands 2 (green), 3 (red), and 4 (near infrared), but mapped onto colors red for band 4, green for band 3 and blue for band 2.  Click on the Red color component of Composite in the table of contents and adjust the band assigned to red to be the band labeled "compositec4".  This is band 4, the 4th geoTIFF grid loaded in to the composite. 

 

Similarly assign Green to band 3 "compositec3" and Blue to band 2 "compositec2".  Notice as you do this how the color display changes as different underlying data is mapped onto a different color used to construct the image. 

 

Zooming in the image should appear similar to the below

 

 

RGB color composite of bands 4, 3, and 2 for Path 29 Row 32.

 

Next we are going to zoom in on the Upper Sand Creek hydrologic unit. 

 

Open the Attribute Table of the Hydrologic Unit layer.

 

 

You will see that there are 100 HUC-12 watersheds in this table.

 

Go to Select by Attributes. Double click on "HU_12_NAME".  Click on =.  Click on Get Unique Values and double click on "Upper Sand Creek".  The result should be the query below.  Click Apply. This will select Upper Sand Creek watershed.

 

 

Click on zoom to selected features to zoom in on Upper Sand Creek

 

 

 

To turn in:

2.      A Landsat false color composite map (bands 4, 3, 2) of Upper Sand Creek HUC 12 watershed.  Use the layout view and include a north arrow, scale bar and title.  Indicate the NRD District that Upper Sand Creek is in.

 

3. EVALUATING NDVI

 

NDVI, the normalized difference vegetation index, is a quantity used to assess the presence of live green vegetation.  NDVI is computed using the formula:    

The RED and NIR stand for the spectral reflectance measurements acquired in the red and near-infrared regions of electromagnetic spectrum, respectively. NDVI takes values from -1 to 1. The higher the NDVI, higher the fraction of live green vegetation present. Landsat band 4 (0.77-0.90 µm) measures the reflectance in NIR region and Band 3 (0.63-0.69 µm) measures the reflectance in the Red region. 

 

In ArcMap select Windows/Image Analysis.  This opens the Image Analysis tool shown

 

 

We can set the ‘Image analysis options’ using the  button on the top left corner of the Image analysis window.

 

Click the button. Under NDVI, make sure Red Band is 3 and Infrared Band is 4. This band combination holds true for Landsat Images. You might have other hypersprectral or satellite images where layer numbers are different. Also if you have stacked only two or three bands in your Composite image, Layer numbers could change accordingly.  Also select Scientific Output.  Click OK

 

The selection of Scientific output ensures that the formula NDVI = ((NIR - Red)/(NIR + Red)) is applied to evaluate NDVI.  If you do not select this, ArcGIS uses a scaling of NDVI between 0 and 200.

 

Make sure the "Composite" layer is selected in the Image Analysis layer list (at the top) and click on the NDVI button in the Processing section near the bottom.

 

 

You will see that a ‘NDVI_Composite’ layer is created in ArcMap with range -1 to 1.  Close the Image Analysis window.

 

Adjust the NDVI_Composite Symbology.

 

To turn in:

3.      A NDVI map of Upper Sand Creek HUC 12 watershed. 

 

To work with this NDVI data we need to project it to the same projection as our Hydrologic Units and Extract it using the Hydrologic Unit layer as a mask. 

 

Open ArcToolbox/ Data Management Tools/ Projections and Transformations/ Raster/ Project Raster.

 

Set inputs as follows

 

 

For the output coordinate system select Layers and NAD_1983_UTM_zone_14N to select a coordinate system consistent with the Hydrologic Unit layer.

 

Clear selected features in ‘Hydrologic_Unit’ feature dataset to make sure you have no features selected (or else the extraction that follows will be for those features only). 

 

Open ArcToolbox/Spatial Analyst Tools/Extraction/Extract by Mask (You can also search for Extract by Mask under the search). Select Input Raster as NDVI-prj, Input Raster or feature mask data as Hydrologic_Unit. Save the Output Raster in the desired location with the name NDVI-Extract.

 

 

I have found that this function sometimes fails and have been able to work around this by closing and re-opening ArcMap and trying again.  The result is an NDVI grid masked by the Hydrologic Unit Layer.

 

 

4. ESTIMATION OF FRACTIONAL VEGETATION COVER

 

The method proposed by Brunsell and Gillies (2003)1 to obtain the fraction of vegetation cover will be used in this exercise. The method scales the NDVI to obtain the fraction of vegetation cover and then scales the fraction between the NDVI of bare soil and a full canopy.

 

N* = (NDVI – NDVI0)/(NDVImax – NDVI0)

 

Where NDVI0 is the bare soil NDVI value for the Landsat scene and NDVImax is the maximum NDVI of the scene corresponding to full cover dense vegetation.  The fraction of cover is then estimated as

 

1Brunsell, N.A., and R. R. Gillies. 2003. Length Scale Analysis of Surface Energy Fluxes Derived from Remote Sensing. Journal of Hydrometeorology, 4, 1212-1219.

 

Fr = (N*)2

 

We will estimate fraction of vegetation cover using the raster calculator. Click on Spatial Analyst/Map Algebra/Raster Calculator. In this image, the NDVI0 is set at 0.14 for bare soil, and NDVImax is estimated to be 0.75. Construct the following formula as an expression in the raster calculator to calculate the fraction of vegetation cover.

 

Square(("NDVI-Extract" - 0.14) / (0.75 -0.14))

 

Name the output raster “FofVeg”

 

 

In layer properties, change symbology to “classified” and use 7 classes with a Red to Green color ramp.

 

 

You can observe the difference in fractions of vegetation cover on the image. Differences are clearly visible between various land uses, and various fields. Also fraction of vegetation cover can be different within a single field (within-field variability).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 


To turn in:

4.      A Fractional Vegetation Cover map of Upper Sand Creek HUC 12 watershed. 

 

5. COMPUTATION OF REFERENCE EVAPOTRANSPIRATION

 

We will estimate the spatial distribution of the Reference Evapotranspiration (ET) for the Landsat overpass date (DOY=251). We have daily reference evapotranspiration data from 12 Automated Weather Data Network (AWDN) stations within and surrounding the study area. These weather stations are operated by the High Plains Regional Climate Center (HPRCC) in Nebraska. Data can be downloaded from the HPRCC website (http://www.hprcc.unl.edu/). We have done this for you and have placed the data in o the Ref_ET Feature Class.

 

Add the Ref_ET Feature Class from Base Map.gdb to ArcMap. Open the attribute table of Ref_ET. The Attribute table contains the longitude, latitude, elevation and daily reference evapotranspiration data for 12 AWDN stations for many days of the year.

 

 

 

We will create a continuous surface of reference evapotranspiration from this point data using the spline Spatial Interpolation technique. On the Spatial Analyst toolbar go to Spatial Analyst/Interpolation/Spline

 

 

 

Select Input Point features as Ref_ET layer. Z value field as DoY251 which represents the reference ET values for day 251. The Spline type should be set as Regularized, and Number of points as 5. For Output Cell Size select the FofVeg grid to inherit the same cell size as this (60 m).  Save the output raster with the name RefET (in the imagery folder). We selected day 251 because the Landsat data we are using was acquired on the September 9, 2005 which is 251st day of the year.

 

 

Following is the result.

 

 

This reference evapotranspiration is the ‘alfalfa reference’ ET, generally abbreviated as ETr.  The alfalfa reference ETr contrasts with the other type of reference ET, which is for clipped grass vegetation.  Grass reference ET is usually abbreviated as ETo.  Because alfalfa is taller and leafier than clipped grass, ETr is typically about 20 to 30% greater than ETo.  ETr represents the upper limit of ET as constrained by environmental energy to convert liquid water to vapor.  This energy comes mostly from solar radiation and air temperature.

 

6.  COMPUTATION OF ACTUAL EVAPOTRANSPIRATION

We will calculate the Actual Evapotranspiration by multiplying fraction of vegetation cover with the reference evapotranspiration. Here, we assumed that the fraction of reference ET, ETrF, is equal to the fraction of vegetation cover. Click on Spatial analyst/raster calculator

 

Type following formula as an expression in the raster calculator to calculate the actual evapotranspiration.

 

"RefET" * "FofVeg"

 

Name this raster ET, and click OK.

 

 

This method for estimating ET from fraction of cover is only approximate because it assumes that the soil surface between vegetation is relatively dry.  If it has recently rained, then there will be some evaporation from the soil between vegetation and the actual ET will be greater than what we have estimated.

 

Zoom in and observe some of the patterns of actual ET.

 

 

7.  EXPLORING ZONAL STATISTICS OF EVAPOTRANSPIRATION FOR CERTIFIED IRRIGATED FIELDS

 

By performing zonal statistics, the statistics for each zone of a zone dataset can be calculated based on the information in a value raster. Using zonal statistics we can explore the distribution of actual evapotranspiration at each center pivot field. Also we can see the statistical differences in ET in each hydrologic unit.

 

Add the subset_of_certified_fields Feature Class from Base Map.gdb to ArcMap.

 

To calculate zonal statistics, click on Spatial Analyst/Zonal Statistics as a Table and enter the following inputs

                                                                                                                         

 

The result is a table that gives ET statistics for each field.

 

 

Let's see how this field ET is related to NDVI.  Click on Spatial Analyst/Zonal Statistics as a Table and enter the following inputs.  The only difference from above is that we are now using NDVI as the value raster and a different output table name

 

 

The mean column in each table gives the mean ET and mean NDVI for each field respectively. 

Open the fieldet table.  Use FIELD_ID to join the fieldndvi table as follows

 

Click No to the following message.

 

Export the resulting table in dbf format.

 

 

Open the resulting DBF file in Excel and produce a plot of ET vs NDVI

 

 

To turn in:

5.      A plot of ET vs NDVI for the certified fields.  Discuss the cause for the pattern in this plot.  Why do the points roughly, but not exactly follow what looks like a parabola? 

6.      Zoom in to the triangle area shown on the ET image below. When you closely examine ET values, you will see that some of center pivot irrigated fields have higher ET (shown as darker blue) while a few of them have low ET (shown as yellow). Examine two center pivot irrigated fields with high and low NDVI values, and provide an approximate value of NDVI, the fraction of vegetation cover, RefET, and ET value for each of field.

 

 

 

 

SUMMARY OF ITEMS TO BE TURNED IN

1.      Prepare an ArcMap layout showing labeled Natural Resources Districts (NRD) within NE. You should include NHD_flowlines, and HUC-12 units. Use a layout view, include a scale bar to indicate distance, a north arrow to indicate direction and labels or legends with units. A snapshot of ‘data view’ of the map is depicted above.

2.      A Landsat false color composite map (bands 4, 3, 2) of Upper Sand Creek HUC 12 watershed.  Use the layout view and include a north arrow, scale bar and title.  Indicate the NRD District that Upper Sand Creek is in.

3.      A NDVI map of Upper Sand Creek HUC 12 watershed. 

4.      A Fractional Vegetation Cover map of Upper Sand Creek HUC 12 watershed. 

5.      A plot of ET vs NDVI for the certified fields.  Discuss the cause for the pattern in this plot.  Why do the points roughly, but not exactly follow what looks like a parabola? 

6.      Zoom in to the triangle area shown on the ET image below. When you closely examine ET values, you will see that some of center pivot irrigated fields have higher ET (shown as darker blue) while a few of them have low ET (shown as yellow). Examine two center pivot irrigated fields with high and low NDVI values, and provide an approximate value of NDVI, the fraction of vegetation cover, RefET, and ET value for each of field.

 

Addendum: Computing Spectral radiance,  Reflectance and NDVI using Landsat bands

 

A more exact method to calculate NDVI is to first calculate the Spectral Radiance and then Reflectance for band 3 (Red) and band 4 (NIR), rather than using the digital number (DN) in the NDVI function. The spectral radiance (L) must be computed for band 3 and band 4  based on the digital number of each individual pixel. Spectral radiance is the outgoing radiation energy of the band as observed at the top of the atmosphere by the satellite. The three steps in order to calculate NDVI for Landsat 5 is provided in details below:

 

 

Step 1. Conversion from Digital Number to Spectral Radiance

 

                                         

 

Where:

                L                             = Spectral radiance at the sensor aperture (watt m-2 ster-1 μm-1)

                Lmax                     = Spectral radiance scaled to Qcalmax (watt m-2 ster-1 μm-1)

                Lmin                      = Spectral radiance scaled to Qcalmin (watt m-2 ster-1 μm-1)

                Qcal                        = Quantized calibrated pixel value = DN

                Qcalmin                = Minimum quantized calibrated pixel value corresponding to Lmin

                Qcalmax               = Maximum quantized calibrated pixel value corresponding to Lmax

 

Lmax and Lmin are given in Table 1. Since our Landsat image is acquired in 2005, we are going to use calibration constants on the right side of colum on Table 1.

 

Qcalmax = 255, and Qcalmin = 0 for Landsat 5, so the equation becomes:

 

 (Eq. 1)

 

 

Table 1. Landsat  -5 TM POSTCALIBRATION DYNAMIC RANGES FOR U.S. PROCESSED NLAPS DATA

(IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 41, NO. 11, NOVEMBER 2003)

 

 

 

 

 

Add the Landsat Band 4 and Band 3 in ESRI GRID data (not the raw TIFF data) to the ArcMap.

On the spatial analyst toolbar, click Spatial Analyst/Option. Set the working directory the same as the folder where you have saved your Band 3 and Band 4 ESRI GRID files.

 

 

Now open the raster calculator. Spatial Analyst/Raster Calculator. Paste the following expression in the raster calculator in order to calculate spectral radiance of band 4,

 

(((221 + 1.51) / (255 - 1)) * ([band4] - 1)) - 1.51

 

Rename your output as Band4_Radiance.  

 

Similarly, paste the following expression in the raster calculator to compute spectral radiance of band 3

 

(((264 + 1.17) / (255 - 1)) * ([band3] - 1)) - 1.17

 

Make sure you rename your output files as Band4_Radiance and Band3_Radiance.  

 

Step 2. Conversion from radiance to reflectance (at-satellite reflectance)

 

 

Where:

            r           = Planetary reflectance (unitless)

            L          = Spectral radiance at the sensor aperture (watt m-2 ster-1 μm-1) (Eq. 1)

            dr         = Inverse square of earth-sun distance (astronomical unit)

            Esun    = Mean solar exoatmospheric irradiances (watt m-2 μm-1, Table 2)

            θ          = Solar zenith angle (degree)

 

 

dr =  1 + 0.033 Cos (DOY *2 * 3.141592654)/365)

DOY (day of year) = 251

 

θ (Theta) = (90-Bheta) where Beta is sun elevation angle provided in Landsat header file. This equation works fairly well for flat terrain (our case).

 

Table 2.  TM SOLAR EXOATMOSPHERIC SPECTRAL IRRADIANCES.

 

 

Now lets calculate the reflectance of each band in arcgis. Paste the following expression in the raster calculator to calculate reflectance values for  band 4,

 

(3.14 * [Band4_Radiance]) / (1036 * 0.758 * (1 + 0.033 * cos((251 * 2 * 3.14) / 365)))

Rename your output as Band4_ Reflectance.  

 

 

Similarly, paste the following expression in the raster calculator to get reflectance of band 3,

 

(3.14 * [Band3_Radiance]) / (1554 * 0.758 * (1 + 0.033 * cos((251 * 2 * 3.14) / 365)))

Rename your output as Band3_ Reflectance.  

 

 

 

Step 3. Computing NDVI using Landsat bands

 

The Normalized Difference Vegetation Index (NDVI) is a simple numerical index to assess the presence of live green vegetation. NDVI is computed using below formula:

 

 

RED and NIR stand for the spectral reflectance measurements acquired in the red and near-infrared regions of electromagnetic spectrum, respectively. NDVI takes the value from -1 to 1. The higher the NDVI, higher the fraction of live green vegetation present in the scene. Landsat band 4 (0.77-0.90 µm) measures the reflectance in NIR region and Band 3 (0.63-0.69 µm) measures the reflectance in Red region. 

 

We already have the refectances of Band 4 and Band3 calculated in step 2 in order to calculate NDVI. Now open the raster calculator. Spatial Analyst/Raster Calculator. Paste the following expression in the raster calculator to estimate NDVI.

 

([Band4_ Reflectance] - [Band3_ Reflectance]) / ([Band4_ Reflectance] + [Band3_ Reflectance])

 

 

 

Rename the output layer as NDVI. 

 

This calculation for NDVI will be more accurate and consistent (with varying sun angle) than the NDVI computed earlier using digital numbers (DN), only.  If you are curious, you can compare values between the two sets of NDVI calculations to view the differences.