Landsat Calculating Surface Reflectance from METRIC Mapping Evapotranspiration

Landsat: Calculating Surface Reflectance (from: METRIC: Mapping Evapotranspiration at high Resolution using Internalized Calibration) University of Idaho, Kimberly, Idaho University of Nebraska-Lincoln

Landsat – Polar Orbiting A new image each 16 days for a specific location

Landsat – Polar Orbiting

What Landsat sees Transmissivity of atmosphere Visible 12 3 Near Infrared 4 5 7 0 0. 4 0. 6 0. 8 1. 2 1. 6 2. 0 2. 4 : (Band 6 is the surface temperature band (not shown)) Wavelength in Microns Land Surface

Obtaining a Landsat image Ø Ø Ø The Landsat images are available from USGSEROS, terrain-corrected (L 1 T) at no cost Images can be previewed and downloaded at http: //glovis. usgs. gov/ The images can also be downloaded from http: //earthexplorer. usgs. gov/

Download from glovis. usgs. gov/

Download from glovis. usgs. gov/

Download from earthexplorer. usgs. gov/

Download from earthexplorer. usgs. gov/ 4. 2. 3. 1.

Download from earthexplorer. usgs. gov/ The images can be downloaded from here as zipped files In some cases the images are not available for immediate download. The images can be ordered and will be available for download after 1 -2 days

Image file name format Different file names are used for the downloaded (zipped) file and for the unzipped files Zipped file: LSDPPPRRRYYYYDOY_. tar. gz Indv. bands: LSDPPPRRRYYYYDOY_BXX. TIF L = Landsat A = T=>TM (Thematic mapper, Landsat 5); E=>ETM (Enhanced TM, Landsat 7); C=> Landsat 8 D = 5 (landsat 5); 7(landsat 7); 8 (landsat 8) PPP = WRS Path RRR = WRS Row (note, that the row information is duplicated in the file name) YYYY = Year; DOY = Day of the Year S = Satellite (5 or 7) B = Band XX = Band number (1 = band 1; 61 = low gain band 6; 62 = high gain band 6) Example: LC 80330372013173 LGN 00_B 1. TIF

Downloaded files The unzipped downloaded file The original downloaded file (from Glovis or Earth Explorer) The unzipped files (individual bands) Header file The files can be unzipped using e. g. 7 -zip from www. 7 -zip. org/

Preparing the Image Ø Ø Ø The Landsat images are available from USGS-EROS at no cost A layered spectral band image is created from downloaded and unzipped files using ERDAS Imagine software. A subset image is created if a smaller area is to be studied.

Locating the Header File Header file with information to be used with METRIC

Header File for Landsat 8 image Date image was processed (important) Path, row LMAX and LMIN for Bands 1 -11 (for METRIC: Band 2 -7, 10 are important

Header File for Landsat 8 image Qcmax = 65535; Qcmin = 1 Radiance multipliers and add Radiance_b = DN* RADIANCE_MULT_BAND_b + RADIANCE_ADD_BAND_b K 1 and K 2 constants for thermal bands Sun azimuth, elevation and relative distance earth-sun

Header File for Landsat 8 image For L 8 images acquired after about May 2016, the co-registration of TIRS with OLI on L 8 is ‘preliminary’ for the first 12 to 16 days following the overpass. After 12 to 16 days, the image is reprocessed using updated TIRS registration information to improve spatial accuracy. The “TIRS SSM Model” flag is changed from “Preliminary” to “Final”.

Solar Radiation and Reflectance “non-reflected” radiation is what is absorbed at the surface and part of the energy balance. Therefore it is important to calculate accurately.

Satellite Sensor rt, b Esunb (varies by band) Top of Atmosphere tin, b tout, b (varies by band) rs, b (varies by band) Land Surface

Disposition of Solar Radiation in the Atmosphere H 2 O, O 2, O 3, N 2 O Indirect H 2 O, O 2, O 3, N 2 O Direct Solar

Radiance Equation for Landsat Lb = Spectral Radiance in band ‘b’, W/m 2/sr/μm 2 LMAX = Maximum W/m /sr/μm in calibration LMIN = Minimum W/m 2/sr/μm in calibration DN = “digital number” (0 -255) Qcalmax = Maximum DN (associated with LMAX (255)) Qcalmin = Minimum DN (associated with LMIN (0 or 1)) or for older Landsat 7 images (and some 5 images): Lb = (Gain × DN) + Bias

Landsat 8 Landsat 7 LMAX and LMIN for Bands 1 -8 (for METRIC: Band 1 -7 are important ) Landsat 7 has two flavors of band 6; use 6_VCID_2 (low gain) for METRIC Landsat 5 and 7: Qcal max =255, Qcal min =1 Landsat 8: Qcal max =65535, Qcal min =1

(BG) Alternative way to get radiance Landsat 8 Landsat 7 Landsat 5

At satellite reflectances ρt, b ESUNb cosθ dr = the t subscript means top-of-atmosphere (i. e. , atsatellite) and the b subscript indicates the specific band number = Potential Solar Radiation in band b (Table 6. 3) = cosine of solar angle from nadir = inverse of square of relative distance from sun to earth For June 25, 2015: Sun elevation angle ( ) = 66. 310, q = (90 - ) = 32. 690 DOY = 176, dr =0. 967

ESUNb for Landsat 5 and 7 W/m 2/μm New* (post-2008 LPGS L 1 T) *from Chander et al. 2009, Remote Sens. Environ. 113: 893 -903 based on Thuillier et al. , 2003, Solar Physics 214: 1 -22 (SOLSPEC – ATLAS 123/EURECA missions)

ESUNb for Landsat 8 • For Landsat 8 ESUN can be calculated from the METADATA file • ESUN values are automatically calculated in the METRIC spreadsheet after reading the METADATA file

ESUNb for Landsat 8 • For Landsat 8 ESUNb can be calculated from the METADATA file: d=Earth Sun Distance All the values are available in the metadata file:

At Surface Reflectance ρs, b METRIC Approach ρt, b at-satellite reflectance for band “b” ρa, b “path” reflectance for band “b” τin, b and τout, b are narrowband transmittances for incoming solar radiation and for surface reflected shortwave radiation

Incoming Transmissivity C 1 -C 5 Pair W Kt qh = Generalized Coefficients fitted to MODTRAN and SMARTS 2 models by Tasumi, Allen and Trezza (2008) = mean atmospheric pressure, k. Pa (= f(elevation)) = precipitable water in atmosphere (= f(near surface vapor pressure from weather station)) = turbidity (clearness) coefficient (default = 1. 0) = solar angle from nadir of horizontal surface Eq. has similar form to broadband t equation of FAO-56, ASCE-EWRI

Outgoing Transmissivity C 1 -C 5 Pair W Kt qh = Generalized Coefficients fitted to MODTRAN model = mean atmospheric pressure, k. Pa (= f(elevation)) = precipitable water in atmosphere (= f(near surface vapor pressure from weather station)) = turbidity (clearness) coefficient (default = 1. 0) = satellite angle from nadir of horizontal surface (0 for Landsat)

Transmittance Coefficients

Path Reflectance

Narrow band Transmittances: Figure 2. Comparison of transmittance for Landsat bands estimated by Equation 6 (Tasumi et al, 2005) vs. transmittance simulated by SMARTS 2 for 100 combinations of sun angle, precipitable water, and land elevation.

Broadband Surface Albedo Wb = weighting coefficient that considers fraction of all potential solar energy at the surface over range represented by specific band. (Wb’s sum to 1. 0) Range for W 5 0 0. 4 0. 6 0. 8 1. 2 1. 6 2. 0 2. 4 Wavelength in Microns Band: 1 2 3 4 5 7 weighting coefficients by Allen et al. 2006 for Landsat 5 and 7 0. 103

Broadband Surface Albedo Comparison with MODTRAN

Broadband transmittance – used for Rso where τB is the transmissivity index for direct beam radiation [unitless] and τD is the transmissivity index for diffuse radiation This is the ASCE (2005) formulation used with the Standardized Penman. Monteith

– Vegetation Indices Normalized difference vegetation index NDVI = (r. NIR - rred) / (r. NIR + rred) Normalized difference water index NDWI = (r. SWIR-1 - rblue) / (r. SWIR-1 + rblue) Normalized difference snow index NDSI = (rblue - r. SWIR-1) / (rblue + r. SWIR-1) Soil adjusted vegetation index SAVI = (1+L)(r. NIR - rred) / (L+r. NIR + rred)

Vegetation Indices: Landsat 5 and 7 used to estimate aerodynamic roughness and thermal emissivity NDVI = (r 4 - r 3) / (r 4 + r 3) NDWI = (r 5 - r 2) / (r 5 + r 2) (Normalized Difference VI) (Normalized Difference Water Index) SAVI = (1 + L) (r 4 - r 3) / (L + r 4 + r 3) (Soil Adjusted VI) For Southern Idaho: L = 0. 1 SAVIID = 1. 1(r 4 - r 3) / (0. 1 + r 4 + r 3) Leaf Area Index (LAI): LAI = 11(SAVIID)3 Current Function: LAI = 11(NDVIs)3 We limit LAI 6. 0

Vegetation Indices: Landsat 8 used to estimate aerodynamic roughness and thermal emissivity NDVI = (r 5 - r 4) / (r 5 + r 4) NDWI = (r 6 - r 3) / (r 6 + r 3) (Normalized Difference VI) (Normalized Difference Water Index) SAVI = (1 + L) (r 5 - r 4) / (L + r 5 + r 4) (Soil Adjusted VI) For Southern Idaho: L = 0. 1 SAVIID = 1. 1(r 5 - r 4) / (0. 1 + r 5 + r 4) Leaf Area Index (LAI): LAI = 11(SAVIID)3 Current Function: LAI = 11(NDVIs)3 We limit LAI 6. 0

(BG) Warning!! NDVI = (r 4 - r 3) / (r 4 + r 3) (Normalized Difference VI) Please Note! that NDVI (and SAVI) are calculated using reflectances and not digital numbers and not radiances. The variables in the equations must be ‘normalized’ reflectances, by definition. Many novices and nonthinkers commonly compute NDVI using DN or radiance. DN is improper because its scale can change over time. In addition, both DN and radiance magnitudes will change with time of year as the sun angle changes. It also changes with time of day. Reflectance is much more constant and consistent. One can use surface reflectance or top-of-atmosphere reflectance in the calculations. Results are usually similar since atmospheric attenuation is similar for both bands 3 and 4. We choose to use top-of-atmosphere in METRIC NDVI to be consistent with many other uses. However, using surface reflectance is probably slightly more consistent. Note also that NDVI computed from different satellite systems like MODIS will not be the same as from Landsat because of differences in band widths and centers.

– Surface Temperature Ø Surface Emissivity Ø Surface Temperature

Surface Emissivity ( o and NB) (for thermal (infrared) radiation) METRIC uses two surface emissivities: 0 is the emissivity for the broad band spectrum and is used to compute the outgoing longwave radiation via the Stephan-Boltzmann function. NB is the emissivity for the narrow band spectrum and is used to compute the surface temperature via the Planck function.

(BG) Surface Emissivity NB = 0. 97 + 0. 0033 LAI; for LAI < 3 0 = 0. 95 + 0. 01 LAI; for LAI < 3 NB = 0. 98 and 0 = 0. 98 when LAI 3 For water; NDVI < 0 and a < 0. 47, NB = 0. 99 and o = 0. 985 For snow; NDVI < 0 and a 0. 47, NB = 0. 99 and o = 0. 985 Note that some bare rock may have emissivity as low at 0. 90. The user can consult various emissivity libraries or measure.

Longwave (Infrared) Radiation in the Atmosphere

Thermal Radiance (Rc) Rp is the path radiance in the 10. 4 – 12. 5 μm band Rsky is the narrow band downward thermal radiation for a clear sky (units are W/(m 2 sr μm) ) (we consider the 1 - NB component that reflects from the surface) For low aerosol conditions Rp=0. 91, τNB=0. 866 and Rsky=1. 32, based on comparisons with MODTRAN in southern Idaho (Allen et al. 2007)

Surface Temperature (Ts) (Planck’s Law) NB is the emissivity for the narrow band spectrum of the Landsat thermal band (10. 45 -12. 42 μm on Landsat 5 and 10. 31 -12. 36 μm on Landsat 7) K 1 and K 2 are constants Rc is thermal radiance emitted from the surface in the narrow band, W/(m 2 sr μm) Ts is surface temperature in K For Landsat 8 check metadata