Skip to main content

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

  1. Setting up your working environment.
  2. Defining functions for converting ASTER L1T data from Digital Number (DN) <- Radiance <- Top of Atmosphere (TOA) Reflectance.
  3. Reading ASTER L1T HDF-EOS file metadata.
  4. Georeferencing raster arrays and defining the Coordinate Reference System (CRS).
  5. Subsetting the Visible-Near Infrared (VNIR) datasets from the .hdf file.
  6. Applying functions on rasters (convert the data stored as DN to Radiance and convert Radiance to TOA Reflectance).
  7. Visualizing and Plotting ASTER L1T data.
    • Generating an RGB composite image.
    • Calculating Normalized Difference Vegetation Index (NDVI).
    • Visualizing Histograms, Calculating Statistics.
  8. Exporting Results to Georeferenced Tagged Image File Format (GeoTIFF) files.
  9. Creating difference maps between two ASTER L1T observations.
  10. 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

Access 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.

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):

Mathematical equation to convert radiance to Top of Atmosphere TOA reflectance.

Radiance = (DN-1) * ucc

Radiance to TOA Reflectance (from Landsat 7 Users Handbook):

Mathematical equation to convert radiance to Top of Atmosphere TOA reflectance.

where:

Legend of mathematical symbols.

For more information on the equations above, see the references at the end of this document.

UCC Values:

Band No.High gainNormal gainLow gain 1
10.6761.6882.25
20.7081.4151.89
3N0.4230.8621.15
3B0.4230.8621.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

## 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. Access more information on PROJ.4 strings.

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)
ASTER L1T Band 3 TOA Reflectance over Santa Barbara County, CA.

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)
Cropped ASTER L1T Band 3 TOA Reflectance over Santa Barbara and Lake Cachuma.

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')
ASTER L1T Pseudo-true color RGB composite showing the Gap Fire (2008) burn scar.

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)
ASTER L1T TOA Reflectance-derived 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))
Histogram of NDVI values over the region of interest.

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)
ASTER L1T NDVI with labels

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)
ASTER L1T NDVI over region of interest from July 17, 2008.

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)
ASTER L1T NDVI over region of interest from March 18, 2016.

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)
NDVI difference map over region of interest.

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)
NDVI difference map over region of interest with custom color map.

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')
ASTER L1T Pseudo-true color RGB composite showing the Gap Fire (2008) burn scar.

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')
ASTER L1T Pseudo-true color RGB composite from March 18, 2016.

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

Details

Last Updated

June 18, 2025

Published on

Aug. 4, 2017

Author

Cole Krehbiel

Data Centers

Land Processes DAAC (LP DAAC)