Details:
Published: Aug. 4, 2017
Working with ASTER L1T Visible and Near Infrared (VNIR) Data in R
This tutorial demonstrates how to work with The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Level 1 Precision Terrain Corrected Registered At-Sensor Radiance ASTER L1T HDF-EOS (.hdf) files in R. ASTER L1T VNIR data can be converted from Radiance stored as Digital Numbers (DN) into Top of Atmosphere (TOA) reflectance, which can then be used to analyze the biophysical properties of the land surface.
In this example, we show how to convert DNs to TOA Reflectance and then calculate the Normalized Difference Vegetation Index (NDVI), a measure of vegetation greenness. We end the tutorial by analyzing seasonal and disturbance-related changes in NDVI by calculating difference maps between two ASTER L1T-derived NDVI images. By the end of the tutorial, the user will have GeoTIFF files for each of the 3 VNIR bands of TOA reflectance as well as RGB and NDVI images generated in the tutorial.
Topics Covered in this Tutorial
- Setting up your working environment.
- Defining functions for converting ASTER L1T data from Digital Number (DN) <- Radiance <- Top of Atmosphere (TOA) Reflectance.
- Reading ASTER L1T HDF-EOS file metadata.
- Georeferencing raster arrays and defining the Coordinate Reference System (CRS).
- Subsetting the Visible-Near Infrared (VNIR) datasets from the .hdf file.
- Applying functions on rasters (convert the data stored as DN to Radiance and convert Radiance to TOA Reflectance).
-
Visualizing and Plotting ASTER L1T data.
- Generating an RGB composite image.
- Calculating Normalized Difference Vegetation Index (NDVI).
- Visualizing Histograms, Calculating Statistics.
-
Exporting Results to Georeferenced Tagged Image File Format (GeoTIFF) files.
- Creating difference maps between two ASTER L1T observations.
- Exporting plots to image files.
Prerequisites
GDAL will need to be installed on your OS in order to execute this tutorial.
R Libraries:
- rgdal - 1.1.10
- raster - 2.5.2
- gdalUtils - 2.0.1.7
- tiff - 0.1.5
- colorspace
- rmarkdown
Disclaimer: The code for this tutorial was tested in the following environment:
- R: Version 3.3.0 (Anaconda 4.1.1) on Windows OS
- Geospatial Data Abstraction Library (GDAL): Versions 2.0.0 and 2.1.0
ASTER L1T Information
Additional information about the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Level 1 Precision Terrain Corrected Registered At-Sensor Radiance (AST_L1T) data product is located here: https://doi.org/10.5067/aster/ast_l1t.003.
The ASTER L1T files used in this tutorial can be downloaded from:
https://e4ftl01.cr.usgs.gov/ASTT/AST_L1T.003/2008.07.17/AST_L1T_00307172008185216_20150524180249_98171.hdf https://e4ftl01.cr.usgs.gov/ASTT/AST_L1T.003/2016.04.18/AST_L1T_00304182016185206_20160419084818_31976.hdf
1. Set Up Working Environment
First, we load the R packages necessary to run the tutorial. If you are missing any of the packages, you can download them using the install.packages('missing package') command in the console.
# Load necessary packages into R
library(rgdal); library(raster); library(gdalUtils); library(tiff); library(colorspace); library(rmarkdown)
require(rgdal)
Next, we point to the directory containing the ASTER L1T files and set it as our working directory. We recommend saving your R Markdown (.Rmd) file in this directory.
# Set input/current working directory
in_dir <- 'c:/ASTER_L1T/'
setwd(in_dir)
# Create and set output directory
out_dir <- paste(in_dir, 'ASTER_L1T_DEMO_output/', sep='')
suppressWarnings(dir.create(out_dir))
# Search and create a list for all ASTER L1T files contained in the working directory:
file_names <- list.files(pattern = 'AST_L1T_.*hdf$')
file_names
## [1] "AST_L1T_00304182016185206_20160419084818_31976.hdf"
## [2] "AST_L1T_00307172008185216_20150524180249_98171.hdf"
For this exercise we will be working with two ASTER scenes over Santa Barbara County, CA. We will start with the earlier file from July 17, 2008.
ASTER L1T file names are constructed as:
'L1T Short Name''Collection Version''Start Date-Time-Group''Production Date-Time-Group'_'Processing Random Number'.hdf.
The 'Start Date-Time-Group' portion of the file name identifies the acquisition date, in MMDDYYYY format.
# Here we set the earlier file (2008) as our file name for the first part of this exercise
file_name <- file_names[2]
file_name
## [1] "AST_L1T_00307172008185216_20150524180249_98171.hdf"
2. Define Functions for Calculations
ASTER L1T data are stored as Digital Numbers (DN). In order to convert DNs to Top of Atmosphere (TOA) reflectance, which is needed in order to generate NDVI, we must first convert to 'At-Sensor Radiance'. We then convert radiance to Top of Atmosphere TOA reflectance.
DN to Radiance (Abrams, 1999):
Radiance to TOA Reflectance (from Landsat 7 Users Handbook):
where:
For more information on the equations above, see the references at the end of this document.
UCC Values:
| Band No. | High gain | Normal gain | Low gain 1 |
|---|---|---|---|
| 1 | 0.676 | 1.688 | 2.25 |
| 2 | 0.708 | 1.415 | 1.89 |
| 3N | 0.423 | 0.862 | 1.15 |
| 3B | 0.423 | 0.862 | 1.15 |
First, we will set up a look-up table for the unit conversion coefficients by constructing a data frame.
# Construct a dataframe for the UCC values:
bands <- c('1', '2', '3N', '3B') # Row Names
gain_names <- c('Band', 'High Gain', 'Normal Gain', 'Low Gain 1') # Column Names
# Set up a matrix with four rows and three columns containing the ucc values
ucc_vals <- matrix( c(0.676, 1.688, 2.25,
0.708, 1.415, 1.89,
0.423, 0.862, 1.15,
0.423, 0.862, 1.15), nrow = 4, ncol = 3, byrow = TRUE)
# Append band names as first column and convert to data frame
ucc <- data.frame( bands, ucc_vals)
# Set gain names as column header names
names(ucc) <- gain_names
ucc
## Band High Gain Normal Gain Low Gain 1
## 1 1 0.676 1.688 2.25
## 2 2 0.708 1.415 1.89
## 3 3N 0.423 0.862 1.15
## 4 3B 0.423 0.862 1.15
We will also need to define the irradiance values used in the calculations. Spectral irradiance values used here are from Thome et al. 2001.
# From Thome et al. 2001, which uses spectral irradiance values from MODTRAN
# Create a list, ordered b1, b2, b3N, b4, b5...b9
irradiance <- c(1848,1549,1114,225.4,86.63,81.85,74.85,66.49,59.85)
# Remove unnecessary variables
rm(bands,gain_names,ucc_vals)
Below is where we define our functions, including functions for converting degrees to radians, DN values to radiance, and radiance to TOA reflectance. Later on, we will show how to pass raster arrays through these functions.
# Write a function to convert degrees to radians
calc_radians <- function(x) {(x * pi) / (180)}
# Write a function to convert DN values to Radiance
calc_radiance <- function(x){(x - 1) * ucc1}
# Write a function to convert Radiance to TOA Reflectance
calc_reflectance <- function(x){(pi * x * (earth_sun_dist^2)) / (irradiance1 * sin(pi * sza / 180))}
3. Read File Metadata
Next, we will use the file name to grab the date and day of year (DOY).
# grab date from the filename and convert to DOY
date <- paste(substr(file_name, 16, 19), substr(file_name, 12, 13), substr(file_name, 14, 15), sep = '-')
doy <- as.numeric(strftime(date, format = '%j'))
# Now that we have the DOY, we can calculate the earth-sun distance using the function 'calc_radians' defined earlier
# Calculate Earth Sun Distance (Achard and D'Souza 1994; Eva and Lambin 1998--see references for more information)
earth_sun_dist <- (1 - 0.01672 * cos(calc_radians(0.9856 * (doy - 4))))
earth_sun_dist
## [1] 1.016343
Here is where we use the Geospatial Data Abstraction Library (GDAL) gdalinfo() command to read the hdf file metadata.
# Use GDAL to read file metadata
md <- gdalinfo(file_name)
# This will give us a preview of the information contained in the very large string object of metadata
head(md)
## [1] "Driver: HDF4/Hierarchical Data Format Release 4"
## [2] "Files: AST_L1T_00307172008185216_20150524180249_98171.hdf"
## [3] "Size is 512, 512"
## [4] "Coordinate System is `'"
## [5] "Metadata:"
## [6] " ASTEROBSERVATIONMODE.1=VNIR1, ON"
# Convert object into a dataframe
md2 <- data.frame(md)
names(md2) <- 'File Metadata'
md2
## File Metadata
## 1 Driver: HDF4/Hierarchical Data Format Release 4
## 2 Files: AST_L1T_00307172008185216_20150524180249_98171.hdf
## 3 Size is 512, 512
## 4 Coordinate System is `'
## 5 Metadata:
## 6 ASTEROBSERVATIONMODE.1=VNIR1, ON
## 7 ASTEROBSERVATIONMODE.2=VNIR2, ON
## 8 ASTEROBSERVATIONMODE.3=SWIR, ON
## 9 ASTEROBSERVATIONMODE.4=TIR, ON
## 10 ASTEROPERATIONMODE=OBSERVATION
## 11 ASTERSCENEID=42, 103, 4
## 12 AVERAGEOFFSET10=-0.321720987558365, -0.131992995738983
## 13 AVERAGEOFFSET11=-0.321720987558365, -0.131992995738983
## 14 AVERAGEOFFSET12=-0.321720987558365, -0.131992995738983
## 15 AVERAGEOFFSET13=-0.321720987558365, -0.131992995738983
## 16 AVERAGEOFFSET14=-0.321720987558365, -0.131992995738983
## 17 AVERAGEOFFSET4=0.0, 0.0
## 18 AVERAGEOFFSET5=0.0, 0.0
## 19 AVERAGEOFFSET6=0.0, 0.0
## 20 AVERAGEOFFSET7=0.0, 0.0
## 21 AVERAGEOFFSET8=0.0, 0.0
## 22 AVERAGEOFFSET9=0.0, 0.0
## 23 AVGCORRELCOEF4=0.0
## 24 AVGCORRELCOEF5=0.0
## 25 AVGCORRELCOEF6=0.0
## 26 AVGCORRELCOEF7=0.0
## 27 AVGCORRELCOEF8=0.0
## 28 AVGCORRELCOEF9=0.0
## 29 BANDSUSED=01023N3B0405060708091011121314
## 30 CALENDARDATE=20080717
## 31 COARSEDEMVERSION=1.00, 1997-09-03, This data is generated from GTOPO30
## 32 CONUNIT1=W/m2/sr/um
## 33 CONUNIT10=W/m2/sr/um
## 34 CONUNIT11=W/m2/sr/um
## 35 CONUNIT12=W/m2/sr/um
## 36 CONUNIT13=W/m2/sr/um
## 37 CONUNIT14=W/m2/sr/um
## 38 CONUNIT2=W/m2/sr/um
## 39 CONUNIT3B=END_OBJECT
## 40 CONUNIT3N=W/m2/sr/um
## 41 CONUNIT4=END_OBJECT
## 42 CONUNIT5=END_OBJECT
## 43 CONUNIT6=END_OBJECT
## 44 CONUNIT7=END_OBJECT
## 45 CONUNIT8=END_OBJECT
## 46 CONUNIT9=END_OBJECT
## 47 CORINTEL1=N/A
## 48 CORINTEL10=Uncorrected Intertelescope Error
## 49 CORINTEL11=Uncorrected Intertelescope Error
## 50 CORINTEL12=Uncorrected Intertelescope Error
## 51 CORINTEL13=Uncorrected Intertelescope Error
## 52 CORINTEL14=Uncorrected Intertelescope Error
## 53 CORINTEL2=N/A
## 54 CORINTEL3B=END_OBJECT
## 55 CORINTEL3N=N/A
## 56 CORINTEL4=END_OBJECT
## 57 CORINTEL5=END_OBJECT
## 58 CORINTEL6=END_OBJECT
## 59 CORINTEL7=END_OBJECT
## 60 CORINTEL8=END_OBJECT
## 61 CORINTEL9=END_OBJECT
## 62 CORPARA1=N/A
## 63 CORPARA10=N/A
## 64 CORPARA11=N/A
## 65 CORPARA12=N/A
## 66 CORPARA13=N/A
## 67 CORPARA14=N/A
## 68 CORPARA2=N/A
## 69 CORPARA3B=END_OBJECT
## 70 CORPARA3N=N/A
## 71 CORPARA4=END_OBJECT
## 72 CORPARA5=END_OBJECT
## 73 CORPARA6=END_OBJECT
## 74 CORPARA7=END_OBJECT
## 75 CORPARA8=END_OBJECT
## 76 CORPARA9=END_OBJECT
## 77 CORRECTIONACHIEVED=Terrain+Precision
## 78 CTHLD4=0.0
## 79 CTHLD5=0.0
## 80 CTHLD6=0.0
## 81 CTHLD7=0.0
## 82 CTHLD8=0.0
## 83 CTHLD9=0.0
## 84 EASTBOUNDINGCOORDINATE=-119.349965967973
## 85 FLYINGDIRECTION=DE
## 86 FUTUREREVIEWDATE=20080705
## 87 GAIN.1=01, HGH
## 88 GAIN.10=09, NOR
## 89 GAIN.2=02, HGH
## 90 GAIN.3=3N, NOR
## 91 GAIN.4=3B, NOR
## 92 GAIN.5=04, NOR
## 93 GAIN.6=05, NOR
## 94 GAIN.7=06, NOR
## 95 GAIN.8=07, NOR
## 96 GAIN.9=08, NOR
## 97 GEOMETRICDBVERSION=03.01, 2003-10-19
## 98 HDFEOSVersion=HDFEOS_V2.17
## 99 IDENTIFIER_PRODUCT_DOI=10.5067/ASTER/AST_L1T.003
## 100 IDENTIFIER_PRODUCT_DOI_AUTHORITY=http://dx.doi.org
## 101 IMAGEDATAINFORMATION1=5731, 5125, 1
## 102 IMAGEDATAINFORMATION10=956, 855, 2
## 103 IMAGEDATAINFORMATION11=956, 855, 2
## 104 IMAGEDATAINFORMATION12=956, 855, 2
## 105 IMAGEDATAINFORMATION13=956, 855, 2
## 106 IMAGEDATAINFORMATION14=956, 855, 2
## 107 IMAGEDATAINFORMATION2=5731, 5125, 1
## 108 IMAGEDATAINFORMATION3B=5731, 5125, 0
## 109 IMAGEDATAINFORMATION3N=5731, 5125, 1
## 110 IMAGEDATAINFORMATION4=0, 0, 0
## 111 IMAGEDATAINFORMATION5=0, 0, 0
## 112 IMAGEDATAINFORMATION6=0, 0, 0
## 113 IMAGEDATAINFORMATION7=0, 0, 0
## 114 IMAGEDATAINFORMATION8=0, 0, 0
## 115 IMAGEDATAINFORMATION9=0, 0, 0
## 116 INCL1=0.675999999046326
## 117 INCL10=0.00688199978321791
## 118 INCL11=0.00677999993786216
## 119 INCL12=0.00658999988809228
## 120 INCL13=0.00569299980998039
## 121 INCL14=0.00522499997168779
## 122 INCL2=0.708000004291534
## 123 INCL3B=0.0
## 124 INCL3N=0.861999988555908
## 125 INCL4=0.0
## 126 INCL5=0.0
## 127 INCL6=0.0
## 128 INCL7=0.0
## 129 INCL8=0.0
## 130 INCL9=0.0
## 131 INPUTGRANULEID=ASTL1A 0807171852160010269001
## 132 INSTRUMENTSHORTNAME=ASTER
## 133 L1TREPROCESSINGACTUAL=not reprocessed
## 134 LOWERLEFT=34.3316816992043, -120.263823002482
## 135 LOWERLEFTM=3.8037600e+06, 1.997100e+05
## 136 LOWERRIGHT=34.3530214355242, -119.330457832703
## 137 LOWERRIGHTM=3.8037600e+06, 2.856600e+05
## 138 MAPORIENTATIONANGLE=0.0
## 139 MAPPROJECTIONNAME=Universal Transverse Mercator
## 140 MEANANDSTD1=104.647094726562, 20.8634662628174
## 141 MEANANDSTD10=1.529042846679688e+03, 137.137466430664
## 142 MEANANDSTD11=1.663153076171875e+03, 151.48616027832
## 143 MEANANDSTD12=1.782587158203125e+03, 159.71842956543
## 144 MEANANDSTD13=2.0595585937500e+03, 168.563079833984
## 145 MEANANDSTD14=2.126138183593750e+03, 163.313018798828
## 146 MEANANDSTD2=90.0983734130859, 26.018196105957
## 147 MEANANDSTD3B=0.0, 0.0
## 148 MEANANDSTD3N=86.3480682373047, 16.3615665435791
## 149 MEANANDSTD4=0.0, 0.0
## 150 MEANANDSTD5=0.0, 0.0
## 151 MEANANDSTD6=0.0, 0.0
## 152 MEANANDSTD7=0.0, 0.0
## 153 MEANANDSTD8=0.0, 0.0
## 154 MEANANDSTD9=0.0, 0.0
## 155 MEASUREMENTPOINTNUMBER10=188
## 156 MEASUREMENTPOINTNUMBER11=188
## 157 MEASUREMENTPOINTNUMBER12=188
## 158 MEASUREMENTPOINTNUMBER13=188
## 159 MEASUREMENTPOINTNUMBER14=188
## 160 MEASUREMENTPOINTNUMBER4=0
## 161 MEASUREMENTPOINTNUMBER5=0
## 162 MEASUREMENTPOINTNUMBER6=0
## 163 MEASUREMENTPOINTNUMBER7=0
## 164 MEASUREMENTPOINTNUMBER8=0
## 165 MEASUREMENTPOINTNUMBER9=0
## 166 MINANDMAX1=59, 255
## 167 MINANDMAX10=692, 1980
## 168 MINANDMAX11=771, 2158
## 169 MINANDMAX12=827, 2299
## 170 MINANDMAX13=1052, 2582
## 171 MINANDMAX14=1127, 2617
## 172 MINANDMAX2=39, 255
## 173 MINANDMAX3B=0, 0
## 174 MINANDMAX3N=29, 219
## 175 MINANDMAX4=0, 0
## 176 MINANDMAX5=0, 0
## 177 MINANDMAX6=0, 0
## 178 MINANDMAX7=0, 0
## 179 MINANDMAX8=0, 0
## 180 MINANDMAX9=0, 0
## 181 MODEANDMEDIAN1=91, 157
## 182 MODEANDMEDIAN10=1478, 1336
## 183 MODEANDMEDIAN11=1620, 1464
## 184 MODEANDMEDIAN12=1742, 1563
## 185 MODEANDMEDIAN13=2006, 1817
## 186 MODEANDMEDIAN14=2088, 1872
## 187 MODEANDMEDIAN2=70, 147
## 188 MODEANDMEDIAN3B=0, 0
## 189 MODEANDMEDIAN3N=78, 124
## 190 MODEANDMEDIAN4=0, 0
## 191 MODEANDMEDIAN5=0, 0
## 192 MODEANDMEDIAN6=0, 0
## 193 MODEANDMEDIAN7=0, 0
## 194 MODEANDMEDIAN8=0, 0
## 195 MODEANDMEDIAN9=0, 0
## 196 MPMETHOD1=UTM
## 197 MPMETHOD10=UTM
## 198 MPMETHOD11=UTM
## 199 MPMETHOD12=UTM
## 200 MPMETHOD13=UTM
## 201 MPMETHOD14=UTM
## 202 MPMETHOD2=UTM
## 203 MPMETHOD3B=END_OBJECT
## 204 MPMETHOD3N=UTM
## 205 MPMETHOD4=END_OBJECT
## 206 MPMETHOD5=END_OBJECT
## 207 MPMETHOD6=END_OBJECT
## 208 MPMETHOD7=END_OBJECT
## 209 MPMETHOD8=END_OBJECT
## 210 MPMETHOD9=END_OBJECT
## 211 NORTHBOUNDINGCOORDINATE=35.0236619068926
## 212 NUMBERGCPCHIPSCORRELATED=77
## 213 NUMBEROFBADPIXELS1=0, 0
## 214 NUMBEROFBADPIXELS10=0, 0
## 215 NUMBEROFBADPIXELS11=0, 0
## 216 NUMBEROFBADPIXELS12=0, 0
## 217 NUMBEROFBADPIXELS13=0, 0
## 218 NUMBEROFBADPIXELS14=0, 0
## 219 NUMBEROFBADPIXELS2=0, 0
## 220 NUMBEROFBADPIXELS3B=0, 0
## 221 NUMBEROFBADPIXELS3N=0, 0
## 222 NUMBEROFBADPIXELS4=0, 0
## 223 NUMBEROFBADPIXELS5=0, 0
## 224 NUMBEROFBADPIXELS6=0, 0
## 225 NUMBEROFBADPIXELS7=0, 0
## 226 NUMBEROFBADPIXELS8=0, 0
## 227 NUMBEROFBADPIXELS9=0, 0
## 228 NUMBEROFMEASUREMENTS10=299
## 229 NUMBEROFMEASUREMENTS11=299
## 230 NUMBEROFMEASUREMENTS12=299
## 231 NUMBEROFMEASUREMENTS13=299
## 232 NUMBEROFMEASUREMENTS14=299
## 233 NUMBEROFMEASUREMENTS4=0
## 234 NUMBEROFMEASUREMENTS5=0
## 235 NUMBEROFMEASUREMENTS6=0
## 236 NUMBEROFMEASUREMENTS7=0
## 237 NUMBEROFMEASUREMENTS8=0
## 238 NUMBEROFMEASUREMENTS9=0
## 239 OFFSET1=-0.675999999046326
## 240 OFFSET10=-0.00688199978321791
## 241 OFFSET11=-0.00677999993786216
## 242 OFFSET12=-0.00658999988809228
## 243 OFFSET13=-0.00569299980998039
## 244 OFFSET14=-0.00522499997168779
## 245 OFFSET2=-0.708000004291534
## 246 OFFSET3B=0.0
## 247 OFFSET3N=-0.861999988555908
## 248 OFFSET4=0.0
## 249 OFFSET5=0.0
## 250 OFFSET6=0.0
## 251 OFFSET7=0.0
## 252 OFFSET8=0.0
## 253 OFFSET9=0.0
## 254 ORBITNUMBER=45645
## 255 PCTIMAGEMATCH4=0
## 256 PCTIMAGEMATCH5=0
## 257 PCTIMAGEMATCH6=0
## 258 PCTIMAGEMATCH7=0
## 259 PCTIMAGEMATCH8=0
## 260 PCTIMAGEMATCH9=0
## 261 PGEVERSION=1.0
## 262 PLATFORMSHORTNAME=Terra
## 263 POINTINGANGLE.1=0.022
## 264 POINTINGANGLE.2=0.033
## 265 POINTINGANGLE.3=0.004
## 266 PROCESSEDBANDS=01023N3BXXXXXXXXXXXX1011121314
## 267 PROCESSINGCENTER=ASTER_LPDAAC
## 268 PROCESSINGFLAG10=1
## 269 PROCESSINGFLAG11=1
## 270 PROCESSINGFLAG12=1
## 271 PROCESSINGFLAG13=1
## 272 PROCESSINGFLAG14=1
## 273 PROCESSINGFLAG4=0
## 274 PROCESSINGFLAG5=0
## 275 PROCESSINGFLAG6=0
## 276 PROCESSINGFLAG7=0
## 277 PROCESSINGFLAG8=0
## 278 PROCESSINGFLAG9=0
## 279 PROCESSINGLEVELID=1T
## 280 PRODUCTIONDATETIME=2015-05-24T23:03:45.000Z
## 281 PROJECTIONPARAMETERS1=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 282 PROJECTIONPARAMETERS10=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 283 PROJECTIONPARAMETERS11=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 284 PROJECTIONPARAMETERS12=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 285 PROJECTIONPARAMETERS13=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 286 PROJECTIONPARAMETERS14=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 287 PROJECTIONPARAMETERS2=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 288 PROJECTIONPARAMETERS3B=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 289 PROJECTIONPARAMETERS3N=6.3781370e+06, 6.356752300e+06, 0.9996, 0.0, -2.04203522483336, 0.0, 5.000000e+05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 290 PROJECTIONPARAMETERS4=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 291 PROJECTIONPARAMETERS5=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 292 PROJECTIONPARAMETERS6=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 293 PROJECTIONPARAMETERS7=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 294 PROJECTIONPARAMETERS8=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 295 PROJECTIONPARAMETERS9=0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
## 296 QAPERCENTINTERPOLATEDDATA=0.0
## 297 QAPERCENTMISSINGDATA=0.000288301613181829
## 298 QAPERCENTOUTOFBOUNDSDATA=0.000288301613181829
## 299 QUADRANTCLOUDCOVERAGE=0, 0, 6, 6
## 300 RADIOMETRICDBVERSION=04.08, 2008-07-05
## 301 RECEIVINGCENTER=EDOS
## 302 RECURRENTCYCLENUMBER=196, 210
## 303 REPROCESSINGACTUAL=not reprocessed
## 304 RESMETHOD1=CC
## 305 RESMETHOD10=CC
## 306 RESMETHOD11=CC
## 307 RESMETHOD12=CC
## 308 RESMETHOD13=CC
## 309 RESMETHOD14=CC
## 310 RESMETHOD2=CC
## 311 RESMETHOD3B=END_OBJECT
## 312 RESMETHOD3N=CC
## 313 RESMETHOD4=END_OBJECT
## 314 RESMETHOD5=END_OBJECT
## 315 RESMETHOD6=END_OBJECT
## 316 RESMETHOD7=END_OBJECT
## 317 RESMETHOD8=END_OBJECT
## 318 RESMETHOD9=END_OBJECT
## 319 SCENECENTER=34.6895299079244, -119.808997
## 320 SCENECENTERMETERS=3.842205699298686e+06, 2.426733743902585e+05
## 321 SCENECLOUDCOVERAGE=3
## 322 SCIENCEREVIEWDATE=20080705
## 323 SENSORNAME.1=VNIR
## 324 SENSORNAME.2=SWIR
## 325 SENSORNAME.3=TIR
## 326 SENSORSHORTNAME=ASTER_VNIR, ASTER_SWIR, ASTER_TIR
## 327 SETTINGTIMEOFPOINTING.1=2008-07-17T18:45:11Z
## 328 SETTINGTIMEOFPOINTING.2=2008-07-17T18:45:10Z
## 329 SETTINGTIMEOFPOINTING.3=2008-07-17T18:45:09Z
## 330 SHORTNAME=AST_L1T
## 331 SIZEMBDATAGRANULE=93.66
## 332 SOLARDIRECTION=125.521348, 68.636272
## 333 SOURCEDATAPRODUCT=ASTL1A 0807171852160902110683, 2015-05-24T23:02:36.000Z, PDS
## 334 SOUTHBOUNDINGCOORDINATE=34.3530214355242
## 335 SPATIALRESOLUTION=15, 30, 90
## 336 SPHEROIDCODE=WGS84
## 337 STANDARDDEVIATIONOFFSET10=0.260870009660721, 0.223257005214691
## 338 STANDARDDEVIATIONOFFSET11=0.260870009660721, 0.223257005214691
## 339 STANDARDDEVIATIONOFFSET12=0.260870009660721, 0.223257005214691
## 340 STANDARDDEVIATIONOFFSET13=0.260870009660721, 0.223257005214691
## 341 STANDARDDEVIATIONOFFSET14=0.260870009660721, 0.223257005214691
## 342 STANDARDDEVIATIONOFFSET4=0.0, 0.0
## 343 STANDARDDEVIATIONOFFSET5=0.0, 0.0
## 344 STANDARDDEVIATIONOFFSET6=0.0, 0.0
## 345 STANDARDDEVIATIONOFFSET7=0.0, 0.0
## 346 STANDARDDEVIATIONOFFSET8=0.0, 0.0
## 347 STANDARDDEVIATIONOFFSET9=0.0, 0.0
## 348 THRESHOLD10=0.699999988079071, 2.0, 3.0, 13.0
## 349 THRESHOLD11=0.699999988079071, 2.0, 3.0, 13.0
## 350 THRESHOLD12=0.699999988079071, 2.0, 3.0, 13.0
## 351 THRESHOLD13=0.699999988079071, 2.0, 3.0, 13.0
## 352 THRESHOLD14=0.699999988079071, 2.0, 3.0, 13.0
## 353 THRESHOLD4=0.0, 0.0, 0.0, 0.0
## 354 THRESHOLD5=0.0, 0.0, 0.0, 0.0
## 355 THRESHOLD6=0.0, 0.0, 0.0, 0.0
## 356 THRESHOLD7=0.0, 0.0, 0.0, 0.0
## 357 THRESHOLD8=0.0, 0.0, 0.0, 0.0
## 358 THRESHOLD9=0.0, 0.0, 0.0, 0.0
## 359 TIMEOFDAY=185216255000Z
## 360 UPPERLEFT=35.0236619068926, -120.291114948846
## 361 UPPERLEFTM=3.8806200e+06, 1.997100e+05
## 362 UPPERRIGHT=35.0455564809338, -119.349965967973
## 363 UPPERRIGHTM=3.8806200e+06, 2.856600e+05
## 364 UTMZONECODE1=11
## 365 UTMZONECODE10=11
## 366 UTMZONECODE11=11
## 367 UTMZONECODE12=11
## 368 UTMZONECODE13=11
## 369 UTMZONECODE14=11
## 370 UTMZONECODE2=11
## 371 UTMZONECODE3B=11
## 372 UTMZONECODE3N=11
## 373 UTMZONECODE4=11
## 374 UTMZONECODE5=11
## 375 UTMZONECODE6=11
## 376 UTMZONECODE7=11
## 377 UTMZONECODE8=11
## 378 UTMZONECODE9=11
## 379 UTMZONENUMBER=11
## 380 WESTBOUNDINGCOORDINATE=-120.263823002482
## 381 Subdatasets:
## 382 SUBDATASET_1_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":VNIR_Swath:ImageData2
## 383 SUBDATASET_1_DESC=[5125x5731] ImageData2 VNIR_Swath (8-bit unsigned integer)
## 384 SUBDATASET_2_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":VNIR_Swath:ImageData1
## 385 SUBDATASET_2_DESC=[5125x5731] ImageData1 VNIR_Swath (8-bit unsigned integer)
## 386 SUBDATASET_3_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":VNIR_Swath:ImageData3N
## 387 SUBDATASET_3_DESC=[5125x5731] ImageData3N VNIR_Swath (8-bit unsigned integer)
## 388 SUBDATASET_4_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":TIR_Swath:ImageData10
## 389 SUBDATASET_4_DESC=[855x956] ImageData10 TIR_Swath (16-bit unsigned integer)
## 390 SUBDATASET_5_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":TIR_Swath:ImageData11
## 391 SUBDATASET_5_DESC=[855x956] ImageData11 TIR_Swath (16-bit unsigned integer)
## 392 SUBDATASET_6_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":TIR_Swath:ImageData12
## 393 SUBDATASET_6_DESC=[855x956] ImageData12 TIR_Swath (16-bit unsigned integer)
## 394 SUBDATASET_7_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":TIR_Swath:ImageData13
## 395 SUBDATASET_7_DESC=[855x956] ImageData13 TIR_Swath (16-bit unsigned integer)
## 396 SUBDATASET_8_NAME=HDF4_EOS:EOS_SWATH:"AST_L1T_00307172008185216_20150524180249_98171.hdf":TIR_Swath:ImageData14
## 397 SUBDATASET_8_DESC=[855x956] ImageData14 TIR_Swath (16-bit unsigned integer)
## 398 SUBDATASET_9_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":0
## 399 SUBDATASET_9_DESC=[11x11] Latitude (64-bit floating-point)
## 400 SUBDATASET_10_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":1
## 401 SUBDATASET_10_DESC=[11x11] Longitude (64-bit floating-point)
## 402 SUBDATASET_11_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":2
## 403 SUBDATASET_11_DESC=[5125x5731] ImageData2 (8-bit unsigned integer)
## 404 SUBDATASET_12_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":3
## 405 SUBDATASET_12_DESC=[5125x5731] ImageData1 (8-bit unsigned integer)
## 406 SUBDATASET_13_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":4
## 407 SUBDATASET_13_DESC=[5125x5731] ImageData3N (8-bit unsigned integer)
## 408 SUBDATASET_14_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":5
## 409 SUBDATASET_14_DESC=[11x11] Latitude (64-bit floating-point)
## 410 SUBDATASET_15_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":6
## 411 SUBDATASET_15_DESC=[11x11] Longitude (64-bit floating-point)
## 412 SUBDATASET_16_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":7
## 413 SUBDATASET_16_DESC=[855x956] ImageData10 (16-bit unsigned integer)
## 414 SUBDATASET_17_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":8
## 415 SUBDATASET_17_DESC=[855x956] ImageData11 (16-bit unsigned integer)
## 416 SUBDATASET_18_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":9
## 417 SUBDATASET_18_DESC=[855x956] ImageData12 (16-bit unsigned integer)
## 418 SUBDATASET_19_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":10
## 419 SUBDATASET_19_DESC=[855x956] ImageData13 (16-bit unsigned integer)
## 420 SUBDATASET_20_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":11
## 421 SUBDATASET_20_DESC=[855x956] ImageData14 (16-bit unsigned integer)
## 422 SUBDATASET_21_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":12
## 423 SUBDATASET_21_DESC=[9602x58] VNIR_Supplement (8-bit unsigned integer)
## 424 SUBDATASET_22_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":13
## 425 SUBDATASET_22_DESC=[71x13] TIR_Supplement_Temp (32-bit unsigned integer)
## 426 SUBDATASET_23_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":14
## 427 SUBDATASET_23_DESC=[71x100x10x8] TIR_Supplement_Chopper (8-bit unsigned integer)
## 428 SUBDATASET_24_NAME=HDF4_SDS:UNKNOWN:"AST_L1T_00307172008185216_20150524180249_98171.hdf":15
## 429 SUBDATASET_24_DESC=[71x935] TIR_Supplement_Encoder (16-bit unsigned integer)
## 430 Corner Coordinates:
## 431 Upper Left ( 0.0, 0.0)
## 432 Lower Left ( 0.0, 512.0)
## 433 Upper Right ( 512.0, 0.0)
## 434 Lower Right ( 512.0, 512.0)
## 435 Center ( 256.0, 256.0)
As you can see by the number of rows contained in the dataframe above, there is a lot of information contained in the file metadata. You can click through the pages in the bottom right hand corner of the table to see more information. But how can we search for the specific information that we are interested in?
Below, we show how to use regular expressions (grep()), an extremely helpful pattern matching function.
Here we search the very large string object to find solar direction, used in the conversion from radiance to reflectance.
# Here we use grep to search for the text string 'SOLARDIRECTION=' inside of the metadata
grep('SOLARDIRECTION=', md)
## [1] 332
So '332' above is telling us that the specific text string we searched for is contained in line 332 of the file metadata.
Here we take it a step further by inserting the grep call into a call to the large string object, which then prints out the info for line 332 of the metadata.
# Need Solar Zenith Angle (SZA)--calculate by searching for solar direction info
sza <- md[grep('SOLARDIRECTION=', md)]
sza
## [1] " SOLARDIRECTION=125.521348, 68.636272"
Ultimately we are only interested in the second value of the solar direction, which states the solar elevation or solar zenith angle (SZA) needed to convert radiance to reflectance. Below we combine a regular expression call with substr() (to clip out the second value in the string) and convert it to a number using as.numeric().
# Parse out the text and only keep the SZA value
sza <- as.numeric((substr(sza, (regexpr(', ', sza) + 2), 10000)))
sza
## [1] 68.63627
It is important to explain regular expressions and subsetting strings, because we use these tools frequently to extract values from the metadata that are needed to perform the conversions from DN to TOA Reflectance. These tools are useful for working with remote sensing data in R in general.
Below we search for the gain designation for each ASTER L1T VNIR band.
# Search for the gain designation for each band
gain_01 <- gsub(' ', '', strsplit(md[grep('GAIN.*=01', md)], ',')[[1]][[2]])
gain_02 <- gsub(' ', '', strsplit(md[grep('GAIN.*=02', md)], ',')[[1]][[2]])
gain_03b <- gsub(' ', '', strsplit(md[grep('GAIN.*=3B', md)], ',')[[1]][[2]])
gain_03n <- gsub(' ', '', strsplit(md[grep('GAIN.*=3N', md)], ',')[[1]][[2]])
gain_01
## [1] "HGH"
The gain designation keywords are:
- HGH = High
- NOR = Normal
- LOW = Low
4. Georeference and Define the CRS
The ASTER L1T data product is derived from ASTER Level 1A data that has been geometrically corrected and reprojected to a north-up Universal Transverse Mercator (UTM) projection. Below we search for the coordinates needed to define the coordinate reference system (CRS) and georeference the raster arrays.
# Search for lower right (LR) and upper left (UL) bounding box values from the metadata
lr <- substr(md[grep('LOWERRIGHTM', md)], 15, 50)
ul <- substr(md[grep('UPPERLEFTM', md)], 14, 50)
# Define LR and UL corners, notice here we offset the pixel center by dividing the spatial resolution (15 m) by two.
ul_y <- as.numeric((substr(ul, 1, (regexpr(', ' , ul) - 1)))) + 7.5
ul_x <- as.numeric((substr(ul, (regexpr(', ' , ul) + 2), 10000))) - 7.5
lr_y <- as.numeric((substr(lr, 1, (regexpr(', ', lr) - 1)))) - 7.5
lr_x <- as.numeric((substr(lr, (regexpr(', ', lr) + 2) , 10000))) + 7.5
# Search for and define Universal Transverse Mercator (UTM) zone from metadata
utm_z <- substr(md[grep('UTMZONECODE', md)[1]], 1, 50)
utm_z <- substr(utm_z, regexpr('=', utm_z) + 1, 50)
# Configure the bounding box (extent) properties
y_min <- min(ul_y, lr_y); y_max <- max(ul_y, lr_y)
x_min <- min(ul_x, lr_x); x_max <- max(ul_x, lr_x)
# Here we define the extent using the bounding box values and the extent() function from the raster package.
raster_dims_15m <- extent(x_min, x_max, y_min, y_max)
raster_dims_15m
## class : Extent
## xmin : 199702.5
## xmax : 285667.5
## ymin : 3803753
## ymax : 3880628
Next, we compile the CRS into a PROJ.4 string to attach projection information to rasters. For more information on PROJ.4 strings, see https://proj4.org//.
crs_string <- paste('+proj=utm +zone=', utm_z, ' +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0', sep = '')
crs_string
## [1] "+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
# Remove unnecessary variables
rm(lr,lr_x,lr_y,md,ul,ul_x,ul_y,utm_z,x_min,x_max,y_min,y_max,doy,date)
Now we can get a list of science dataset (SDS) names using get_subdatasets() from the gdalUtils package.
sds <- get_subdatasets(file_name)
sds
## [1] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:VNIR_Swath:ImageData2"
## [2] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:VNIR_Swath:ImageData1"
## [3] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:VNIR_Swath:ImageData3N"
## [4] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:TIR_Swath:ImageData10"
## [5] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:TIR_Swath:ImageData11"
## [6] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:TIR_Swath:ImageData12"
## [7] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:TIR_Swath:ImageData13"
## [8] "HDF4_EOS:EOS_SWATH:AST_L1T_00307172008185216_20150524180249_98171.hdf:TIR_Swath:ImageData14"
## [9] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:0"
## [10] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:1"
## [11] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:2"
## [12] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:3"
## [13] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:4"
## [14] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:5"
## [15] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:6"
## [16] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:7"
## [17] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:8"
## [18] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:9"
## [19] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:10"
## [20] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:11"
## [21] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:12"
## [22] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:13"
## [23] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:14"
## [24] "HDF4_SDS:UNKNOWN:AST_L1T_00307172008185216_20150524180249_98171.hdf:15"
Using the information above, we can set which SDS coincides with the correct band in the hdf file. The ASTER VNIR bands are:
- 1 (green)
- 2 (red)
- 3N (NIR)
5. Subset the VNIR SDS
First, we set the band designations and generate output filenames for each SDS.
band1 <- sds[2]; band2 <- sds[1]; band3 <- sds[3]
# Clean up the names of the specific SDS
clip2 <- c(max(unlist((gregexpr(':', band1)))), max(unlist((gregexpr(':', band2)))), max(unlist((gregexpr(':', band3)))))
# Generate output filename
new_file_name <- strsplit(file_name, '.hdf')
# Define band names by combining the filename with the band name
b1_name <- paste(new_file_name, substr(band1, (clip2[1] + 1), 10000), sep = '_')
b2_name <- paste(new_file_name, substr(band2, (clip2[2] + 1), 10000), sep = '_')
b3_name <- paste(new_file_name, substr(band3, (clip2[3] + 1), 10000), sep = '_')
# Add output directory to the filenames. We will later use this to export the output files
band1_tif_name <- paste(out_dir, b1_name, '.tif', sep='')
band2_tif_name <- paste(out_dir, b2_name,'.tif', sep='')
band3_tif_name <- paste(out_dir, b3_name,'.tif', sep='')
band1_tif_name
## [1] "c:/ASTER_L1T/ASTER_L1T_DEMO_output/AST_L1T_00307172008185216_20150524180249_98171_ImageData1.tif"
Below we extract the specific datasets using gdal_translate() and specify which SDS to open by setting the 'sd_index' parameter:
# Extract specified SDS and export as GeoTIFF (leaving in DN)
invisible(gdal_translate(file_name, band1_tif_name, sd_index=2))
invisible(gdal_translate(file_name, band2_tif_name, sd_index=1))
invisible(gdal_translate(file_name, band3_tif_name, sd_index=3))
Here we use the raster() command to set the CRS using the PROJ.4 string that we created earlier. We then use the extent() call to define the correct extent within the specified CRS.
# Set CRS for each band
aster_b1 <- raster(band1_tif_name, crs = crs_string)
aster_b2 <- raster(band2_tif_name, crs = crs_string)
aster_b3 <- raster(band3_tif_name, crs = crs_string)
# Define Extent
extent(aster_b1) <- raster_dims_15m
extent(aster_b2) <- raster_dims_15m
extent(aster_b3) <- raster_dims_15m
aster_b1
## class : RasterLayer
## dimensions : 5125, 5731, 29371375 (nrow, ncol, ncell)
## resolution : 15, 15 (x, y)
## extent : 199702.5, 285667.5, 3803753, 3880628 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : c:\ASTER_L1T\ASTER_L1T_DEMO_output\AST_L1T_00307172008185216_20150524180249_98171_ImageData1.tif
## names : AST_L1T_00307172008185216_20150524180249_98171_ImageData1
## values : 0, 255 (min, max)
# Remove unnecessary variables
rm(clip2, sds, b1_name, b2_name, b3_name, band1, band2, band3)
6. Apply Functions on Rasters
Below we determine the ucc and irradiance values specific to each band, and ultimately apply the calc_radiance() and calc_reflectance() functions (defined earlier) to each SDS.
# Create an if-else statement to define High/Normal/Low Gain setting for band 1
if (gain_01 == 'HGH'){
ucc1 <- ucc[1, 2]
}else if(gain_01 == 'NOR'){
ucc1 <- ucc[1, 3]
}else{
ucc1 <- ucc[1, 4]
}
# Set irradiance for band 1
irradiance1 <- irradiance[1]
Here we use the calc() function from the raster package, defining which raster to apply the function on and which function to execute.
# Convert from DN to Radiance for band 1
rad_b1 <- calc(aster_b1, calc_radiance)
# Set zeros to avoid performing calculation on missing data
rad_b1[rad_b1 == calc_radiance(0)] <- 0
# Convert from Radiance to TOA Reflectance for band 1
ref_b1 <- calc(rad_b1, calc_reflectance)
rm(rad_b1)
Now perform the same steps for bands 2 and 3:
# Logic to define High/Normal/Low Gain setting for band 2
if (gain_02 == 'HGH'){
ucc1 <- ucc[2, 2]
}else if(gain_02 == 'NOR'){
ucc1 <- ucc[2, 3]
}else{
ucc1 <- ucc[2, 4]
}
# Set irradiance for band 2
irradiance1 <- irradiance[2]
# Convert from DN to Radiance for band 2
rad_b2 <- calc(aster_b2, calc_radiance)
# Set zeros to avoid performing calculation on missing data
rad_b2[rad_b2 == calc_radiance(0)] <- 0
# Convert from Radiance to TOA Reflectance for band 2
ref_b2 <- calc(rad_b2, calc_reflectance)
rm(rad_b2)
# Logic to define High/Normal/Low Gain setting for band 3
if (gain_03n == 'HGH'){
ucc1 <- ucc[3, 2]
}else if(gain_03n == 'NOR'){
ucc1 <- ucc[3, 3]
}else{
ucc1 <- ucc[3, 4]
}
# Set irradiance for band 3
irradiance1 <- irradiance[3]
# Convert from DN to Radiance for band 3
rad_b3 <- calc(aster_b3, calc_radiance)
# Set zeros to avoid performing calculation on missing data
rad_b3[rad_b3 == calc_radiance(0)] <- 0
# Convert from Radiance to TOA Reflectance for band 3
ref_b3 <- calc(rad_b3, calc_reflectance)
rm(rad_b3)
7. Visualize ASTER L1T Data
Next we provide some examples of visualizing ASTER L1T data in R.
Below we use the basic plot() function, coming from the raster package, to visualize band 3 reflectance.
# Basic Visualization of Data
plot(ref_b3)
Figure 1. ASTER L1T Band 3 TOA Reflectance over Santa Barbara County, CA.
Next, we will crop our ASTER L1T image to a smaller extent for our region of interest (ROI).
NOTE: If you are using this tutorial on your own file/region of interest, you will need to change the extent values specific to your region of interest.
# Set up extent for output rasters using a bounding box (in meters)
y_min <- 3814342.897
y_max <- 3835932.940
x_min <- 222887.657
x_max <-251859.590
# Define the new extent
subextent <- extent(x_min, x_max, y_min, y_max)
# Crop layer to our subset region of interest (ROI) using the 'crop()' function from the raster package.
aster_b1 <- crop(aster_b1, subextent)
aster_b2 <- crop(aster_b2, subextent)
aster_b3 <- crop(aster_b3, subextent)
ref_b1 <- crop(ref_b1, subextent)
ref_b2 <- crop(ref_b2, subextent)
ref_b3 <- crop(ref_b3, subextent)
# Set fill value equal to NA so it is ignored in plotting/statistics
aster_b1[aster_b1 == 0] <- NA
aster_b2[aster_b2 == 0] <- NA
aster_b3[aster_b3 == 0] <- NA
# Basic plot for one of the cropped layers:
plot(ref_b3)
Figure 2. Cropped ASTER L1T Band 3 TOA Reflectance over Santa Barbara and Lake Cachuma.
Now, we can stack the three VNIR bands using the stack() function from the raster package, and plot an RGB 'natural color' composite of our ROI. Below, we pass the RGB raster stack into the function plotRGB(), define which band goes into the Red, Green, and Blue designations, set the stretch using the scale parameter, set the maximum # of pixels to plot, define the NA color as white, and set the stretch type to linear.
# Stack all three VNIR bands to generate RGB composite using the stack() function from the raster package
bands_123 <- stack(aster_b2, aster_b3, aster_b1)
rm(aster_b1, aster_b2, aster_b3)
# Remove the intermediate files
invisible(file.remove(band1_tif_name, band2_tif_name, band3_tif_name))
# Make RGB composite of VNIR bands using the plotRGB() function from the raster package
plotRGB(bands_123, 1,2,3, scale = 800, maxpixels = 2000000, colNA='white', stretch = 'lin')
Figure 3. ASTER L1T Pseudo-true color RGB composite showing the Gap Fire (2008) burn scar.
Above we can see part of Santa Barbara, CA and the Pacific Ocean in the bottom of the image. Moving north, we see the burn scar from the 2008 Gap Fire, and the green 'spine' of the Santa Ynez Mountains. In the upper left hand corner of the image we see Lake Cachuma--note that the reservoir looks full in this image.
Next, we define a function to compute NDVI using band 2 (red) and band 3 (NIR) reflectance, and then plot our results.
# Write a function to calculate NDVI from Red & NIR Reflectance
calc_NDVI <- function(x,y){(x - y)/ (x + y)} # NDVI = (NIR - Red)/(NIR + Red)
# Calculate NDVI by calling NDVI function and designating NIR and Red bands
aster_ndvi <- calc_NDVI(ref_b3, ref_b2)
# Visualize NDVI
plot(aster_ndvi)
Figure 4. ASTER L1T TOA Reflectance-derived NDVI over region of interest from July 17, 2008.
Above, regions with higher NDVI appear green, notably the vegetated Santa Ynez Mountains. If we want to stretch the image, a good way to determine where to set the stretch bounds is by plotting the histogram, which we do below by calling the hist() function, which then computes the frequency of values using the entire raster array. We can customize the aesthetics of our graph by calling parameters within the hist() function.
Here we add a title using the 'main=' parameter, define the color of the histogram using the 'col=' parameter, etc.
# Generate histogram to help stretch the color map:
hist(aster_ndvi, main="Distribution of NDVI Values", col= "steelblue4", plot = TRUE,
maxpixels=ncell(aster_ndvi), breaks = 250, xlab = 'NDVI', ylim = c(0,50000), xlim = c(-0.2, 0.8))
Figure 5. Histogram of NDVI values over the region of interest.
Now that we see the distribution of NDVI values in our ROI, we can plot NDVI and set the stretch bounds using the 'zlim' parameter:
# Plot NDVI for the clipped ROI, setting zlim to reflect the distribution of NDVI values above:
plot(aster_ndvi, zlim=c(-0.1,0.65), xlab = 'UTM Longitude (m)', ylab = 'UTM Latitude (m)', maxpixels = 20000000)
title("ASTER L1T NDVI: July 17, 2008\nSanta Barbara, CA", line = 0.21)
Figure 6. ASTER L1T NDVI over region of interest from July 17, 2008 with title and labels.
Calculate Statistics
Here we use cellStats() from the Raster package to calculate statistics on an entire raster array.
# Mean
print(paste0('Mean NDVI:', cellStats(aster_ndvi, stat='mean', na.rm=TRUE)))
## [1] "Mean NDVI:0.284815870316154"
# Min and Max
print(paste0('Min NDVI:', cellStats(aster_ndvi, stat='min', na.rm=TRUE)))
## [1] "Min NDVI:-0.0899943114904258"
print(paste0('Max NDVI:', cellStats(aster_ndvi, stat='max', na.rm=TRUE)))
## [1] "Max NDVI:0.65735568432493"
# Standard Deviation
print(paste0('Standard Deviation of NDVI:', cellStats(aster_ndvi, stat='sd', na.rm=TRUE)))
## [1] "Standard Deviation of NDVI:0.11893559583146"
8. Export Results to GeoTIFFs
Below we export our results to GeoTIFFs by calling the writeRaster() command.
# Set up output filenames
b1_ref_out_name <- gsub('.tif', '_reflectance.tif', band1_tif_name)
b2_ref_out_name <- gsub('.tif', '_reflectance.tif', band2_tif_name)
b3_ref_out_name <- gsub('.tif', '_reflectance.tif', band3_tif_name)
ndvi_outname <- paste(out_dir, new_file_name, '_NDVI.tif', sep = '')
rgb_outname <- paste(out_dir, new_file_name, '_RGB.tif', sep = '')
# Export the raster files to the output directory
writeRaster(ref_b1, filename = b1_ref_out_name, options = 'INTERLEAVE=BAND',
NAflag = 0, format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(ref_b2, filename = b2_ref_out_name, options = 'INTERLEAVE=BAND',
NAflag = 0, format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(ref_b3, filename = b3_ref_out_name, options = 'INTERLEAVE=BAND',
NAflag = 0, format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(aster_ndvi, filename = ndvi_outname, options = 'INTERLEAVE=BAND',
format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(bands_123, filename = rgb_outname, options = 'INTERLEAVE=BAND',
format = 'GTiff', datatype = 'INT2S', overwrite = TRUE)
Repeat Steps 1-8 for the Second ASTER L1T File
Below is where we look at comparing two ASTER L1T images by creating a difference map between them. But first, we need to re-run all of the previous commands on the second file, and then use the output (NDVI) for comparison. The second file is over the same region, but it is from March 18, 2016.
# Here we have grouped all of the commands for steps 1-8 into a single codeblock.
# For comments and explanations of the lines below, see the code blocks above.
file_name2 <- file_names[1]
date <- paste(substr(file_name2, 16, 19), substr(file_name2, 12, 13), substr(file_name2, 14, 15), sep = '-')
doy <- as.numeric(strftime(date, format = '%j'))
earth_sun_dist <- (1 - 0.01672 * cos(calc_radians(0.9856 * (doy - 4))))
md <- gdalinfo(file_name2)
sza <- md[grep('SOLARDIRECTION=', md)]
sza <- as.numeric((substr(sza, (regexpr(', ', sza) + 2), 10000)))
gain_01 <- gsub(' ', '', strsplit(md[grep('GAIN.*=01', md)], ',')[[1]][[2]])
gain_02 <- gsub(' ', '', strsplit(md[grep('GAIN.*=02', md)], ',')[[1]][[2]])
gain_03b <- gsub(' ', '', strsplit(md[grep('GAIN.*=3B', md)], ',')[[1]][[2]])
gain_03n <- gsub(' ', '', strsplit(md[grep('GAIN.*=3N', md)], ',')[[1]][[2]])
lr <- substr(md[grep('LOWERRIGHTM', md)], 15, 50)
ul <- substr(md[grep('UPPERLEFTM', md)], 14, 50)
ul_y <- as.numeric((substr(ul, 1, (regexpr(', ' , ul) - 1)))) + 7.5
ul_x <- as.numeric((substr(ul, (regexpr(', ' , ul) + 2), 10000))) - 7.5
lr_y <- as.numeric((substr(lr, 1, (regexpr(', ', lr) - 1)))) - 7.5
lr_x <- as.numeric((substr(lr, (regexpr(', ', lr) + 2) , 10000))) + 7.5
utm_z <- substr(md[grep('UTMZONECODE', md)[1]], 1, 50)
utm_z <- substr(utm_z, regexpr('=', utm_z) + 1, 50)
y_min <- min(ul_y, lr_y); y_max <- max(ul_y, lr_y)
x_min <- min(ul_x, lr_x); x_max <- max(ul_x, lr_x)
raster_dims_15m <- extent(x_min, x_max, y_min, y_max)
crs_string <- paste('+proj=utm +zone=', utm_z, ' +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0', sep = '')
rm(lr,lr_x,lr_y,md,ul,ul_x,ul_y,utm_z,x_min,x_max,y_min,y_max,doy,date)
sds <- get_subdatasets(file_name2)
band1 <- sds[2]; band2 <- sds[1]; band3 <- sds[3]
clip2 <- c(max(unlist((gregexpr(':', band1)))), max(unlist((gregexpr(':', band2)))), max(unlist((gregexpr(':', band3)))))
new_file_name <- strsplit(file_name2, '.hdf')
b1_name <- paste(new_file_name, substr(band1, (clip2[1] + 1), 10000), sep = '_')
b2_name <- paste(new_file_name, substr(band2, (clip2[2] + 1), 10000), sep = '_')
b3_name <- paste(new_file_name, substr(band3, (clip2[3] + 1), 10000), sep = '_')
band1_tif_name <- paste(out_dir, b1_name, '.tif', sep='')
band2_tif_name <- paste(out_dir, b2_name,'.tif', sep='')
band3_tif_name <- paste(out_dir, b3_name,'.tif', sep='')
invisible(gdal_translate(file_name2, band1_tif_name, sd_index=2))
invisible(gdal_translate(file_name2, band2_tif_name, sd_index=1))
invisible(gdal_translate(file_name2, band3_tif_name, sd_index=3))
aster_b1 <- raster(band1_tif_name, crs = crs_string)
aster_b2 <- raster(band2_tif_name, crs = crs_string)
aster_b3 <- raster(band3_tif_name, crs = crs_string)
extent(aster_b1) <- raster_dims_15m
extent(aster_b2) <- raster_dims_15m
extent(aster_b3) <- raster_dims_15m
rm(clip2, sds, b1_name, b2_name, b3_name, band1, band2, band3)
if (gain_01 == 'HGH'){
ucc1 <- ucc[1, 2]
}else if(gain_01 == 'NOR'){
ucc1 <- ucc[1, 3]
}else{
ucc1 <- ucc[1, 4]
}
irradiance1 <- irradiance[1]
rad_b1 <- calc(aster_b1, calc_radiance)
rad_b1[rad_b1 == calc_radiance(0)] <- 0
ref_b1 <- calc(rad_b1, calc_reflectance)
rm(rad_b1)
if (gain_02 == 'HGH'){
ucc1 <- ucc[2, 2]
}else if(gain_02 == 'NOR'){
ucc1 <- ucc[2, 3]
}else{
ucc1 <- ucc[2, 4]
}
irradiance1 <- irradiance[2]
rad_b2 <- calc(aster_b2, calc_radiance)
rad_b2[rad_b2 == calc_radiance(0)] <- 0
ref_b2 <- calc(rad_b2, calc_reflectance)
rm(rad_b2)
if (gain_03n == 'HGH'){
ucc1 <- ucc[3, 2]
}else if(gain_03n == 'NOR'){
ucc1 <- ucc[3, 3]
}else{
ucc1 <- ucc[3, 4]
}
irradiance1 <- irradiance[3]
rad_b3 <- calc(aster_b3, calc_radiance)
rad_b3[rad_b3 == calc_radiance(0)] <- 0
ref_b3 <- calc(rad_b3, calc_reflectance)
y_min <- 3814342.897
y_max <- 3835932.940
x_min <- 222887.657
x_max <-251859.590
subextent <- extent(x_min, x_max, y_min, y_max)
aster_b1 <- crop(aster_b1, subextent)
aster_b2 <- crop(aster_b2, subextent)
aster_b3 <- crop(aster_b3, subextent)
ref_b1 <- crop(ref_b1, subextent)
ref_b2 <- crop(ref_b2, subextent)
ref_b3 <- crop(ref_b3, subextent)
aster_b1[aster_b1 == 0] <- NA
aster_b2[aster_b2 == 0] <- NA
aster_b3[aster_b3 == 0] <- NA
bands_123_2 <- stack(aster_b2, aster_b3, aster_b1)
rm(aster_b1, aster_b2, aster_b3)
invisible(file.remove(band1_tif_name, band2_tif_name, band3_tif_name))
aster_ndvi2 <- calc_NDVI(ref_b3, ref_b2)
b1_ref_out_name <- gsub('.tif', '_reflectance.tif', band1_tif_name)
b2_ref_out_name <- gsub('.tif', '_reflectance.tif', band2_tif_name)
b3_ref_out_name <- gsub('.tif', '_reflectance.tif', band3_tif_name)
ndvi_outname <- paste(out_dir, new_file_name, '_NDVI.tif', sep = '')
rgb_outname <- paste(out_dir, new_file_name, '_RGB.tif', sep = '')
writeRaster(ref_b1, filename = b1_ref_out_name, options = 'INTERLEAVE=BAND',
NAflag = 0, format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(ref_b2, filename = b2_ref_out_name, options = 'INTERLEAVE=BAND',
NAflag = 0, format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(ref_b3, filename = b3_ref_out_name, options = 'INTERLEAVE=BAND',
NAflag = 0, format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(aster_ndvi2, filename = ndvi_outname, options = 'INTERLEAVE=BAND',
format = 'GTiff', datatype = 'FLT8S', overwrite = TRUE)
writeRaster(bands_123_2, filename = rgb_outname, options = 'INTERLEAVE=BAND',
format = 'GTiff', datatype = 'INT2S', overwrite = TRUE)
9. Create Difference Maps
Now that both images are georeferenced and converted from DN -> TOA Reflectance -> NDVI, we can compare and contrast the images.
First, we can look at the earlier image:
# Plot NDVI for the clipped ROI
plot(aster_ndvi, zlim=c(-0.1,0.75), xlab = 'UTM Longitude (m)', ylab = 'UTM Latitude (m)', maxpixels = 20000000)
title("ASTER L1T NDVI: July 17, 2008\nSanta Barbara, CA", line = 0.21)
Figure 7. ASTER L1T NDVI over region of interest from July 17, 2008.
Next, the later image:
# Plot NDVI for the clipped ROI
plot(aster_ndvi2, zlim=c(-0.1,0.75), xlab = 'UTM Longitude (m)', ylab = 'UTM Latitude (m)', maxpixels = 20000000)
title("ASTER L1T NDVI: March 18, 2016\nSanta Barbara, CA", line = 0.21)
Figure 8. ASTER L1T NDVI over region of interest from March 18, 2016.
Based on the images above, we can see that the ROI is much greener in spring 2016 than in summer 2008. Also note how much Lake Cachuma has shrunk by 2016, and the regeneration of vegetation over the area where the Gap Fire occurred in 2008.
# Perform band math, subtracting the older image from the newer image to look at changes between the observations.
ndvi_change <- aster_ndvi2 - aster_ndvi
# Generate a basic plot of the changes:
plot(ndvi_change)
Figure 9. NDVI difference map over region of interest.
Finally, we can clean up the difference map, define a more meaningful colormap, and export our plot to an image file.
# Create a custom color ramp
brgr = diverge_hcl(50, h = c(55,160), l = c(17,83),c = 90, p = 1.9)
# Plot NDVI change
plot(ndvi_change, zlim=c(-0.7,0.7), main="NDVI Change: 07/17/2008 to 04/18/2016",xlab = 'UTM Longitude (m)', ylab = 'UTM Latitude (m)', maxpixels = 20000000, col = brgr)
Figure 10. NDVI difference map over region of interest with custom color map.
# Export difference map to a png image file.
invisible(dev.copy(png, paste0(out_dir,'NDVI_difference_map.png')))
invisible(dev.off())
Above, green indicates increases in vegetation greenness and brown indicates decreases in vegetation greenness between the two observations. We see positive NDVI change over the 2008 Gap Fire burn scar, as well as areas where vegetation has grown in where Lake Cachuma has receded. There are a few areas that had higher NDVI in July of 2008--possibly irrigated agricultural fields.
We can visualize both of the RGB composites once more to confirm our findings:
# Make RGB composite of VNIR bands from 2008 image:
plotRGB(bands_123, 1,2,3, scale = 800, maxpixels = 2000000, colNA='white', stretch = 'lin')
Figure 11. ASTER L1T Pseudo-true color RGB composite from July 17, 2008.
# Make RGB composite of VNIR bands from 2016 image:
plotRGB(bands_123_2, 1,2,3, scale = 800, maxpixels = 2000000, colNA='white', stretch = 'lin')
Figure 12. ASTER L1T Pseudo-true color RGB composite from March 18, 2016.
If you are interested in automated batch processing scripts for ASTER L1T files, check out the LP DAAC Data Prep Scripts
References
- ABRAMS, M., HOOK, S., and RAMACHANDRAN, B., 1999, Aster user handbook, Version 2, NASA/Jet Propulsion Laboratory, Pasadena, CA, at https://asterweb.jpl.nasa.gov/content/03_data/04_Documents/aster_user_guide_v2.pdf.
- ARCHARD, F., AND D'SOUZA, G., 1994, Collection and pre-processing of NOAA-AVHRR 1km resolution data for tropical forest resource assessment. Report EUR 16055, European Commission, Luxembourg, at http://bookshop.europa.eu/en/collection-and-pre-processing-of-noaa-avhrr-1-km-resolution-data-for-tropical-forest-resource-assessment-pbCLNA16055/.
- EVA, H., AND LAMBIN, E.F., 1998, Burnt area mapping in Central Africa using ATSR data, International Journal of Remote Sensing, v. 19, no. 18, 3473-3497, at http://dx.doi.org/10.1080/014311698213768.
- THOME, K.J., BIGGAR, S.F., and SLATER, P.N., 2001, Effects of assumed solar spectral irradiance on intercomparisons of earth-observing sensors. In International Symposium on Remote Sensing, International Society for Optics and Photonics, pp. 260-269, at http://dx.doi.org/10.1117/12.450668.
- Landsat 7 Science Data Users Handbook: https://landsat.gsfc.nasa.gov/wp-content/uploads/2016/08/Landsat7_Handbook.pdf.
Contact Information
Material written by Cole Krehbiel1
Contact: LPDAAC@usgs.gov
Voice: +1-605-594-6116
Organization: Land Processes Distributed Active Archive Center (LPDAAC)
Website: https://lpdaac.usgs.gov/
Date last modified: 07-13-2017
1 Innovate!, Inc., contractor to the U.S. Geological Survey, Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, USA. Work performed under USGS contract G15PD00403 for LP DAAC2.
2 LP DAAC Work performed under NASA contract NNG14HH33I.