# Lesson 3. Work with MODIS remote sensing data in in R.

## Learning Objectives

After completing this tutorial, you will be able to:

• Open MODIS imagery in R
• Create NBR index using MODIS imagery in R
• Calculate total burned area in R

## What you need

You will need a computer with internet access to complete this lesson and the data for week 6 of the course.

First, let’s import MODIS data. Below notice that we have used a slightly different version of the list.files() pattern= argument.

We have used glob2rx("*sur_refl*.tif$") to select all layers that both 1. Contain the word sur_refl in them and 2. Contain the extension .tif Let’s import our MODIS image stack. # open modis bands (layers with sur_refl in the name) all_modis_bands_july7 <-list.files("data/week06/modis/reflectance/07_july_2016/crop", pattern=glob2rx("*sur_refl*.tif$"),
full.names = T)
# create spatial raster stack
all_modis_bands_st_july7 <- stack(all_modis_bands_july7)
## Error in x[[1]]: subscript out of bounds

# view range of values in stack
all_modis_bands_st_july7[[2]]

# view band names
names(all_modis_bands_st_july7)
# clean up the band names for neater plotting
names(all_modis_bands_st_july7) <- gsub("MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b", "Band",
names(all_modis_bands_st_july7))

# view cleaned up band names
names(all_modis_bands_st_july7)


## Reflectance values range 0-1

As we’ve discussed in class, the normal range of reflectance values is 0-1 where 1 is the BRIGHTEST values and 0 is the darkest value. Have a close look at the min and max values in the second raster layer of our stack, above. What do you notice?

The min and max values are widely outside of the expected range of 0-1 - min: -32768, max: 32767 What could be causing this? We need to better understand our data before we can work with it more. Have a look at the table in the MODIS users guide. The data that we are working with is the MOD09GA product. Look closely at the table on page 14 of the guide. Part of the table can be seen below.

Click here to check out the MODIS user guide - check out page 14 for the MOD09GA table.

The column headings for the table below:

GroupScience Data Sets (HDF Layers (21))UnitsData TypeFill ValueValid RangeScale Factor
surf_Refl_b01: 500m Surface Reflectance Band 1 (620-670 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001
surf_Refl_b02: 500m Surface Reflectance Band 2 (841-876 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001
surf_Refl_b03: 500m Surface Reflectance Band 3 (459-479 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001
surf_Refl_b04: 500m Surface Reflectance Band 4 (545-565 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001
surf_Refl_b05: 500m Surface Reflectance Band 5 (1230-1250 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001
surf_Refl_b06: 500m Surface Reflectance Band 6 (1628-1652 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001
surf_Refl_b07: 500m Surface Reflectance Band 7 (2105-2155 nm)Reflectance16-bit signed integer-28672-100 - 160000.0001

Looking at the table, answer the following questions

1. What is valid range of values for our data?
2. What is the scale factor associated with our data?

## Explore our data

Looking at histograms of our data, we can see that the range of values is not what we’d expect. We’d expect values between -100 to 10000 yet instead we have much larger numbers.

# turn off scientific notation
options("scipen"=100, "digits"=4)
# bottom, left, top and right
#par(mfrow=c(4, 2))
hist(all_modis_bands_st_july7,
col = "springgreen",
xlab = "Reflectance Value")
## Error in hist(all_modis_bands_st_july7, col = "springgreen", xlab = "Reflectance Value"): object 'all_modis_bands_st_july7' not found
mtext("Distribution of MODIS reflectance values for each band\n Data not scaled",
outer = TRUE, cex = 1.5)
## Error in mtext("Distribution of MODIS reflectance values for each band\n Data not scaled", : plot.new has not been called yet


## Scale Factor

Looking at the metadata, we can see that our data have a scale factor. Let’s deal with that before doing anything else. The scale factor is .0001. This means we should multiple each layer by that value to get the actual reflectance values of the data.

We can apply this math to all of the layers in our stack using a simple calculation shown below:

# deal with nodata value --  -28672
all_modis_bands_st_july7 <- all_modis_bands_st_july7 * .0001
# view histogram of each layer in our stack
# par(mfrow=c(4, 2))
hist(all_modis_bands_st_july7,
xlab = "Reflectance Value",
col = "springgreen")
## Error in hist(all_modis_bands_st_july7, xlab = "Reflectance Value", col = "springgreen"): object 'all_modis_bands_st_july7' not found
mtext("Distribution of MODIS reflectance values for each band\n Scale factor applied", outer = TRUE, cex = 1.5)
## Error in mtext("Distribution of MODIS reflectance values for each band\n Scale factor applied", : plot.new has not been called yet


Great - now the range of values in our data appear more reasonable. Next, let’s get rid of data that are outside of the valid data range.

## NoData Values

Next, let’s deal with no data values. We can see that our data have a “fill” value of -28672 which we can presume to be missing data. But also we see that valid range values begin at -100. Let’s set all values less than -100 to NA to remove the extreme negative values that may impact out analysis.

# deal with nodata value --  -28672
all_modis_bands_st_july7[all_modis_bands_st_july7 < -100 ] <- NA
## Error in all_modis_bands_st_july7[all_modis_bands_st_july7 < -100] <- NA: object 'all_modis_bands_st_july7' not found
#par(mfrow=c(4,2))
# plot histogram
hist(all_modis_bands_st_july7,
xlab = "Reflectance Value",
col = "springgreen")
## Error in hist(all_modis_bands_st_july7, xlab = "Reflectance Value", col = "springgreen"): object 'all_modis_bands_st_july7' not found
mtext("Distribution of reflectance values for each band", outer = TRUE, cex = 1.5)
## Error in mtext("Distribution of reflectance values for each band", outer = TRUE, : plot.new has not been called yet


Next we plot MODIS layers. Use the MODIS band chart to figure out what bands you need to plot to create a RGB (true color) image.

BandWavelength range (nm)Spatial Resolution (m)Spectral Width (nm)
Band 1 - red620 - 6702502.0
Band 2 - near infrared841 - 8762506.0
Band 3 - blue/green459 - 4795006.0
Band 4 - green545 - 5655003.0
Band 5 - near infrared1230 – 12505008.0
Band 6 - mid-infrared1628 – 165250018
Band 7 - mid-infrared2105 - 215550018

In the plot below, i’ve called attention to the AOI boundary with a yellow color. Why is it so hard to figure out where the study area is in this MODIS image?

## Error in ogrListLayers(dsn = dsn): Cannot open data source

## Error in plotRGB(all_modis_bands_st_july7, r = 1, g = 4, b = 3, stretch = "lin", : object 'all_modis_bands_st_july7' not found
## Error in box(col = "white"): plot.new has not been called yet
## Error in plot(fire_boundary_sin, add = T, border = "yellow", lwd = 50): object 'fire_boundary_sin' not found


Next, we can deal with clouds in the same way that we dealt with them using Landsat data. However, our cloud mask in this case is slightly different with slightly different cloud cover values as follows:

StateTranslated ValueCloud Condition
000clear
011cloudy
102mixed
113not set, assumed clear

The metadata for the MODIS data are a bit trickier to figure out. If you are interested, the link to the MODIS user guide is below.

The MODIS data are also stored natively in a H4 format which we will not be discussing in this class. For the purposes of this assignment, use the table above to assign cloud cover “values” and to create a mask.

Use the cloud cover layer data/week06/modis/reflectance/07_july_2016/crop/cloud_mask_july7_500m to create your mask.

Set all values >0 in the cloud cover layer to NA.

## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)
## Error in strwidth(legend, units = "user", cex = cex, font = text.font): plot.new has not been called yet

## Error in mask(all_modis_bands_st_july7, cloud_mask_7July): object 'all_modis_bands_st_july7' not found


Plot the masked data. Notice that now the clouds are gone as they have been assigned the value NA.

## Error in plotRGB(all_modis_bands_st_mask, r = 1, g = 4, b = 3, stretch = "lin", : object 'all_modis_bands_st_mask' not found
## Error in box(col = "white"): plot.new has not been called yet
## Error in plot(fire_boundary_sin, add = T, col = "yellow", lwd = 1): object 'fire_boundary_sin' not found


Finally crop the data to see just the pixels that overlay our study area.

## Error in crop(all_modis_bands_st_mask, fire_boundary_sin): object 'all_modis_bands_st_mask' not found
## Error in box(col = "white"): plot.new has not been called yet

SEVERITY LEVELNBR RANGE
Enhanced Regrowth-700 to -100
Unburned-100 to +100
Low Severity+100 to +270
Moderate Severity+270 to +660
High Severity+660 to +1300
## Error in overlay(all_modis_bands_st_mask[[2]], all_modis_bands_st_mask[[7]], : object 'all_modis_bands_st_mask' not found
## Error in plot(modis_nbr_cl, main = "MODIS NBR for the Cold Springs site", : object 'modis_nbr_cl' not found
## Error in strwidth(legend, units = "user", cex = cex, font = text.font): plot.new has not been called yet


After we’ve calculated NBR, we may want to calculate total burn AREA. We can do this using the freq() function in R. This function gives us the total number of pixels associated with each value in our classified raster.

1. Calculate frequency ignoring NA values: freq(modis_nbr_cl, useNA='no')
2. Calculate frequency, ignore NA & only include values that equal 5: freq(modis_nbr_cl, useNA='no', value=5)

Let’s use the MODIS data from 7 July 2016 to calculate the total area of land classified as:

1. Burn: moderate severity
2. Burn: high severity
# get summary counts of each class in raster
freq(modis_nbr_cl, useNA='no')

final_burn_area_high_sev <- freq(modis_nbr_cl, useNA='no', value=5)
## Error in freq(modis_nbr_cl, useNA = "no", value = 5): object 'modis_nbr_cl' not found
final_burn_area_moderate_sev <- freq(modis_nbr_cl, useNA='no', value=4)
## Error in freq(modis_nbr_cl, useNA = "no", value = 4): object 'modis_nbr_cl' not found

## Error in x[[1]]: subscript out of bounds
## Error in all_modis_bands_st_july17[all_modis_bands_st_july17 < -100] <- NA: object 'all_modis_bands_st_july17' not found
## Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. (file does not exist)


We can perform the steps that we performed above, on the MODIS post-fire data too. Below is a plot of the July 17 data.

par(col.axis="white", col.lab="white", tck=0)
# clouds removed
1,4,3,
stretch="lin",
main = "Final data plotted with mask\n Post Fire - 17 July 2016",
axes=T)
box(col = "white")
## Error in box(col = "white"): plot.new has not been called yet


Next we calculate NBR on our post fire data. Then we can crop and finally plot the final results!

## Error in overlay(all_modis_bands_st_mask_july17[[2]], all_modis_bands_st_mask_july17[[7]], : object 'all_modis_bands_st_mask_july17' not found


## Post fire NBR results

the_colors = c("palevioletred4","palevioletred1","ivory1")
barplot(modis_nbr_july17_cl,
main = "Distribution of burn values - Post Fire",
col=rev(the_colors),
names.arg=c("Low Severity","Moderate Severity","High Severity"))
## Error in barplot(modis_nbr_july17_cl, main = "Distribution of burn values - Post Fire", : object 'modis_nbr_july17_cl' not found


Finally, plot the reclassified data. Note that we only have 3 classes: 2, 3 and 4.

## Error in plot(modis_nbr_july17_cl, main = "MODIS NBR for the Cold Springs site \n Post fire", : object 'modis_nbr_july17_cl' not found