Published: Aug. 4, 2017
Visualizing and Plotting ASTER L1T data.
Exporting Results to Georeferenced Tagged Image File Format (GeoTIFF) files.
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
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)
# 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"
'L1T Short Name''Collection Version''Start Date-Time-Group''Production Date-Time-Group'_'Processing Random Number'.hdf.
# 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"
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 |
# 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
# 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)
# 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))}
# 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
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)
grep()
), an extremely helpful pattern matching function.# Here we use grep to search for the text string 'SOLARDIRECTION=' inside of the metadata
grep('SOLARDIRECTION=', md)
## [1] 332
# Need Solar Zenith Angle (SZA)--calculate by searching for solar direction info
sza <- md[grep('SOLARDIRECTION=', md)]
sza
## [1] " SOLARDIRECTION=125.521348, 68.636272"
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.
# 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"
# 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
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)
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"
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"
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))
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)
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]
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)
# 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)
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.
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.
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.
# 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.
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.# 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.
# 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.
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"
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)
# 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)
# 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.
# 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.
# 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.
# 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())
# 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.
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.