Estimating Evapotranspiration using Landsat 5 TM and ArcGIS9.3

Prepared by Parikshit Ranade and Ayse Irmak

GIS in Water Resources

Fall 2009

Purpose

The purpose of this exercise is to utilize raw Landsat 5 TM data in order to display a false color composite (FCC) bands of the Landsat scene (path 29, Row 32) and calculate vegetation index (NDVI) in ArcGIS 9.3 environment. The FCC will be used to explore various features  in the study area.  The Normalized Difference Vegetation index (NDVI) will be calculated from the Landsat 5 TM bands in order to compute  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 also be studied.

Computer and Data Requirements

To carry out this exercise, you need to have a computer, which runs the ArcInfo version of ArcGIS 9.3. There is no need for a remote sensing software. The data are provided in the accompanying zip file, located at  http://www.ce.utexas.edu/prof/maidment/giswr2009/ex6/ex6.zip  (138MB)   These files can also be accessed at civil.ce.utexas.edu\class\class_files\MAIDMENT\giswr2009\Ex6.

Downloading Landsat data

Landsat data can be freely downloaded from http://edcsns17.cr.usgs.gov/EarthExplorer/ or http://glovis.usgs.gov/. USGS global visualization viewer facilitates viewing and downloading of satellite imagery. USGS provides Landsat 5 TM data in TIFF format. This data can be converted into other formats using ArcGIS or any image processing software (i.e. ERDAS-IMAGINE, ENVI).

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 southeastern part of Nebraska and Landsat passes over this area at approximately 10:15 CST. Each Landsat scene has a unique identifier. Our dataset has identifier LT50290322005251EDC00. This data is 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 comprises of Landsat scene with 0% cloud cover. All the aforementioned information regarding Landsat data can be found in Header file of Landsat data.  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 Landsat scene is also provided in the header file. A screen shot of attributes for this Landsat scene is provided below.

1.      Converting raw Landsat 5 TM data (TIFF format) to ERSI GRID 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 designation 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

Various composites can be mapped to facilitate viewing of the raw imagery. A false color composite of near infrared is useful for observing spatial distribution of the vegetation in the Landsat scene by combining bands 2, 3 & 4 (RGB). We are going to generate the color composite by using spectral bands 2 (green), 3 (red), and 4 (near infrared).

Now we are going to import raw Landsat data (TIFF format) into ArcGIS and convert it into GRID format. Open ArcGIS and click add data button to select Landsat data. Browse to the location where you have saved the raw Landsat 5 TM data. Open band 4, band 3 and band 2. These are the files ending with B40.TIF, B30.TIF, B20.TIF for band 4, band 3 and band 2, respectively.

Do not create pyramids and click NO for following message box. Since we are not going to use this data again, building pyramids is unnecessary for this exercise. Pyramids are useful if you are going to explore the data and display it again in varying resolution.

You can see the three Landsat bands (similar to the following image) in your ArcGIS window. This data is in TIFF format and cannot be utilized in ArcGIS to generate near infrared RGB composite. We need to convert this data to ESRI GRID format.

Right click on the Landsat band 4 data and click Data/Export to export TIFF data into ESRI GRID format.

Select the Location where you want to save this dataset as ESRI GRID format. Keep the other settings the same. This means leave the Extent, Spatial Reference and Cell Size (30 by 30m) the same as original dataset (default). Name your dataset as “Band4” and select Format as GRID.  Click Save. Below is a screen shot to export band 4.

Now you have exported the TIFF data into GRID data. Click YES to display this data as a map layer.

Repeat the same procedure to export Band 3 and Band 2 data into GRID format and rename the files as Band3, Band2, respectively. Display them in Arcmap. Make sure you save all the files (band 4, Band 3, and Band 2 grids) at the same location.

2.      Creating false color composite (FCC) of Landsat imagery

We will now create a false near infrared RGB composite to study the vegetation distribution of the Landsat scene. Since it is false color composite, vegetation is red due to its high reflectance in the near infrared spectral band, which was used as the red component in the color composite. Chlorophyll reflects energy in NIR (near infrared) band and vegetation is vividly depicted as varying shades of red.

To create FC,  first, enable spatial analyst extension (Tools/Extensions/Spatial Analyst). On the Spatial Analyst toolbar, click on Options. Under the General tab, select a working Directory and hit OK button. The working directory should be the same as the folder where you have saved your Band4, Band3, Band2 GRID data.

Now click Spatial Analyst/Raster Calculator. We will use MAKESTACK command (ArcInfo command) to create a GRID stack of three bands. The syntax for the command is following:

MAKESTACK newstack LIST [Band1] [Band2] [Band3]

Where newstack is “the name of the output grid stack”. The Band 1, Band 2 and Band3 are the names of raster files for the three bands.

You can paste following expression in Raster Calculator window.

MAKESTACK FCC LIST [Band4] [Band3] [Band2]

Now RGB composite with a name entitled “FCC” was created in your working directory. Click YES to display this in ArcMap. Following image will be displayed in your ArcMap window. Explore this RGB composite closely. If you zoom in the image, you will see Red and Greenish circles. These are the center pivot fields. Dark sheds of red color indicates high percentage of green vegetation while green indicates lack of green vegetation. The green color fields may be fallow land. It is also possible that the crops are harvested/reached maturity in these fields. You will also be able to see some other features in the image such as urban settlement, roads, streams/rivers etc.

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

3.      Exploring the RGB composite of a Hydrologic Unit

After exploring the RGB image, you will have an idea about the geographical area. This is southeastern part of Nebraska where majority of landcover is agriculture. Now we will try to locate some agricultural fields irrigated with center pivot on the Upper Sand Creekwatershed.

Add Hydrologic Unit shapefile to the ArcMap. You can see the HUC 12 level hydrologic watersheds. We will select a particular watershed (Sand creek) and identify the center pivot locations in the watershed.

HUC 12 level hydrologic watersheds on Path 29 Row 32

Open the Attribute Table of Hydrologic Unit shapefile.

Go to Options/Select by Attributes. Select the Upper Sand Creek watershed by its FID. In the window paste the expression "FID" = 28 and click on “Apply”. This will select Upper Sand Creekwatershed.

Now zoom in the Upper Sand Creekwatershed to explore the agricultural fields and center pivot locations within the area. You can see some urban settlement (bright shed of green) in the middle of the watershed. You can also see a highway passing through the center of the watershed from east to west. Try to locate the fields with center pivot irrigation system. The fields should have circular shapes with different shades/tones of blue and red color. Try to pick up the field with clear circular shape.

Now you know where your center pivot fields are in the Upper Sand Creekwatershed.  Add cities, landuse on the arcmap document to explore the study area. This is a landuse map for 2005 developed by CALMIT at UNL  (calmit.unl.edu)

4.      Uploading GPS Locations for Center Pivot Fields in the Upper Sand Creekwatershed

We will create a shapefile for center pivot locations. First, let’s create an empty shapefile and then add the center pivot data into it. Open ArcCatalog. Locate the folder where you have saved your exercise data (your working directory). Right click and select New/Shapefile.

Name the shapefile as “Center_Pivot”.  Select feature type as Point. Click Edit to select the co-ordinate system. We are going to use the coordinate system of the landuse raster. Click Import, browse for the folder where you have land use (Land_Use). Now your shapefile will have NAD 1983 UTM Zone 14 co-ordinate system.  Now drag/add the “Center_Pivot” shapefile to the ArcMap.

We have already provided GPS locations for some center pivot fields within the Landsat scene. It is provided as shape file with a name entitled  “GPS_locations” which shows x, y coordinates for nine center pivot irrigated corn and soybean fields. Add the GPS_locations  shapefile to the ArcMap. We will create the polygons for these fields. First, we will digitize the center points of the nine fields and then will create buffer around them in order to create an approximate field boundary for each of them.

Locations of agricultural fields with center pivot irrigation system.

Next, bring “editor” toolbar to the ArcGIS window.  Click on View/Toolbars/Editor.

Select Editor/Start Editing. Select Center_Pivot.shp for editing. This is the empty file that you created in Arc Catalog. We are going to create a point file (Center_Pivot.shp ) which will display the middle/center locations of nine fields irrigated with center pivot irrigation system. Therefore, we need to “Zoom in” as much as possible on the GPS_locations  file in order to identify middle/center point of nine agricultural fields. The existing GPS locations were taken randomly from the fields without considering the center points of fields.

On the editor toolbar, select sketch tool. Click on the middle of a center pivot field and souble click. This will create a point feature at the center of the field.

If you think you have misplaced the center/middle point, you can select the point feature again using Select Tool . And then delete it using DELETE key on keyboard. Once you are done with selecting few center pivot locations, click Editor/Stop Editing. Click Yes to save the edits to the shapefiles.

What you have obtained now is center point of the center pivot fields.

5.      Creating buffers around GPS _Locations to draw polygons for the center pivot locations

Lets create the polygons (field boundaries) for these center pivot fields. Open ArcToolbox/Analysis Tools/ Proximity/ Buffer. We will create a 400 m buffer around each location to obtain the field polygons for the center pivot locations.

Select input feature as Center_Pivot shapefile. Browse output location and save the new Polygon as Center_Pivot_Fields. Select the Linear Unit value as 400 meter.

You will see the circles (polygon feature) around the center pivot fields. Zoom in to see if the polygon boundary has matched well with the actual field boundary. You can delete few polygons if they are not matching with the actual fields. You can use Editor Toolbar to delete these polygons. We used 400 for the average radius of the fields in this area. We can select other size based on trial and error or using Measure Tool .

6.     Computing Spectral radiance,  Reflectance and NDVI using Landsat bands

In order to calculate NDVI, we first need to calculate the Spectral Radiance and then Reflectance for band 3 (Red) and band 4 (NIR). 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 ESRI GRID data (not the raw TIFF data) to the ArcMap.

On the spatial analyst toolbar, click Spatial Analyst/Option. Set the working directory 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.

Now we will re-project the NDVI raster. Open ArcTolbox/ Data Management Tools/ Projections and Transformations/ Raster/ Project Raster.

Select Input Raster as NDVI. Save the re-projected data in the desired location with a  name NDVI-project. Click to Output co-ordinate system. Click Import, Select Geographic Transformation as NAD_1983_To_WGS1984_1. Now our NDVI data has same co-ordinate system as other datasets.

You will see the following window after completion of the transformation.

The NDVI data that we have derived has very high and low values (anomaly). We will extract this NDVI data for our watersheds within the Landsat scene. Open ArcToolbox/Spatial AnalystTools/Extraction/Extract by Mask. Select Input Raster as NDVI-project, Input Raster or feature mask data as Hydrologic_Unit. Save the Output Raster in the desired location with the name NDVI_Extract. Add this data layer to your ArcMap.

Now you will see the neat NDVI map with reasonable range.

7.    Estimation of Fractional Vegetation Cover

The method proposed by Brunsell and Gillies (2002) 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 emissivity of bare soil and a full canopy.

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

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

Fr = N*2

We will estimate fraction of vegetation cover using raster calculator. Click on Spatial Analyst/Raster Calculator. Type following formula as an expression in the raster calculator to calculate the fraction of vegetation cover.

Sqr(([NDVI_Extract] - 0.14) / (0.81 - 0.14))

Rename the output raster as Fraction_of_Veg_Cover. Export this raster file (Fraction_of_Veg_Cover)  to the desired location with the name Veg_Frac.

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

8.    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 High Plains Regional Climate Center (HPRCC) in Nebraska. Data can be downloaded from HPRCC website (http://www.hprcc.unl.edu/).

Add the Ref_ET.shp to ArcMap. Open the attribute table. Attribute table contains the longitude, latitude, elevation and daily reference evapotranspiration data for 11 AWDN stations.

We will create continuous surface of reference evapotranspiration from this point data using Spatial Interpolation technigue. On Spatial Analyst toolbar go to Spatial Analyst/Interpolate to Raster/Inverse Distance Weighted…

Select Input Points as Ref_ET layer. Z value as DoY251 which represents the reference ET values for day 251. The Power should be set as 2, Search Radius type as Variable and Number of points as 5. Change the Output Cell Size to 30. You can save the output raster to the desired location with the name Ref_ET_Raster. Landsat data is acquired on the September 9, 2005 which is 251st day of the year. So we selected the DOY as 251 to create the reference evapotranspiration surface.

User can define interpolation parameters based on the input data. We are using commonly used power value 2 for IDW interpolation method. Selecting higher power will emphasis more on nearest points. Specifying a lower value for power will provide a bit more influence to surrounding points a little farther away. Higher power will lead to “Bull’s eye” effect, i. e. higher values for the locations surrounding the sample points. Kriging is more sophisticated method of interpolation available with spatial analyst. But selection of parameters is crucial and critical for accuracy of interpolation.

Observe the distribution of the reference evapotranspiration over the watersheds.

9.    Computation of Actual Evapotranspiration

We will calculate the Actual Evapotranspiration by multiplying fraction of vegetation cover with the reference evapotranspiration. Click on Spatial analyst/raster calculator

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

[Vegatation_Fraction] * [Ref_ET_Raster]

Click OK and then actual ET layer will be added to the ARC map. Rename this raster  as ET. Export this layer with name ET to the desired location.

Once you have the ET map, observe the patterns of actual evapotranspiration within different landuses by using landuse data layer, within different fields and within same fields. Do you see any similarity between pattern of fraction of vegetation and pattern of actual evapotranspiration?

10.                       Exploring Zonal Statistics of Evapotranspiration

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 we have marked. Also we can see the statistical differences in ET in each hydrologic unit. To calculate zonal statistics, click on Spatial Analyst/Zonal Statistics

Select Zonal Dataset as Center Pivot Buffer. Select Zonal Field as CPID. Select Value Raster as ET. Save the Output Table at desired location with a name CP_Zonal_Stat.

Now we will see the differences in the ET in each center pivot field.

For the second part, select Zonal Dataset as Hydrologic_Unit. Select Zonal Field as HUC_12_NAME. Select Value Raster as ET. Save the Output Table at desired location with name HUC_12_Zonal_Stat.

Now we will see the differences in the ET in each of our watershed (HUC12 unit).

12 Items to be turned in

Q1. Which is the biggest urban settlement in the study area?

Q2. Which is biggest city in the located in Beaver Creek watershed (note that Beaver Creek watershed has 5 sub watersheds?

Q3. Which center pivot field (specify the CPID) has the highest mean ET value?  Which center pivot field has the lowest ET mean values? What are the crop types in both fields? (Hint-Note that Value field in zonal statistics is same as the CPID field in Center_Pivot_Buffer attribute table)

Q4. Which crop type has higher average ET for the selected 9 center pivot fields? Try to speculate the reasons for that.

Q5. Run the Zonal Satistics for NDVI value raster with HUC 12 watershed as a Zonal dataset. Show the graphical results.

Q6. Run the Zonal Satistics for Fraction od vegetation cover value value raster with Center Pivot GPS locations as a Zonal dataset. Show the graphical results.

Q7. Plot the relationship between mean NDVI and Mean ET for all the 9 center pivot locations.

Q8. Plot the relationship between mean fractions of vegetation cover and mean ET for all the HUC 12 watersheds.

Q9. Observe the distribution of the ET for “City of York” watershed. What speculations can you make?. Use zonal statistics to support your speculations about distribution of ET in the watershed. What factors are contributing to spatial distribution of ET over the watershed?

Q10. Which watershed has the highest percentage of Range, Pasture, and Grass land? Can you support your conclusion using RGB False Color Composite? How?

Q11. Alfalfa is one of the reference crops to compute evapotranspiration. What do you observe about NDVI values of alfalfa crop in the entire study area? Are they high or low compare to surrounding areas? What does is imply?