Procedure

Install and/or load necessary packages and set your working directory**

ipak <- function(pkg){
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    if (length(new.pkg)) 
        install.packages(new.pkg, dependencies = TRUE)
    sapply(pkg, require, character.only = TRUE)
}


packages <- c("rgdal","raster")
ipak(packages)

Prepare Workspace

Like for every new session of R, you should set your working directory.

# Set the working directory ("wd")
setwd("YOURDRIVE\\yourpath\\yourfolder")

Load Data

Load all MODIS bands that end in sur_refl or “surface reflectance”.

# open modis bands (layers with sur_refl in the name)
all_modis_bands_july7 <- list.files("reflectance\\07_july_2016\\crop",
                                    pattern = glob2rx("*sur_refl*.tif$"),
                                    full.names = TRUE)

Stack

Stack and brick all bands and change the names to a more meaningful way.

# create spatial raster stack
all_modis_bands_pre_st <- stack(all_modis_bands_july7)
all_modis_bands_pre_br <- brick(all_modis_bands_pre_st)

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

## class      : RasterLayer 
## dimensions : 2400, 2400, 5760000  (nrow, ncol, ncell)
## resolution : 463.3127, 463.3127  (x, y)
## extent     : -10007555, -8895604, 3335852, 4447802  (xmin, xmax, ymin, ymax)
## crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
## source     : memory
## names      : MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b02_1 
## values     : -1e+06, 100390000  (min, max)

# view band names
names(all_modis_bands_pre_br)

## [1] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b01_1"
## [2] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b02_1"
## [3] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b03_1"
## [4] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b04_1"
## [5] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b05_1"
## [6] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b06_1"
## [7] "MOD09GA.A2016189.h09v05.006.2016191073856_sur_refl_b07_1"

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

# view cleaned up band names
names(all_modis_bands_pre_br)

## [1] "Band01_1" "Band02_1" "Band03_1" "Band04_1" "Band05_1" "Band06_1" "Band07_1"

Explore your data

Plot the histogram of pixel distribution of every MODIS band.

# turn off scientific notation
options("scipen" = 100, "digits" = 4)

# plot the histograms of every band
# use x11() for better viewing expierience
hist(all_modis_bands_pre_br,
     col = "springgreen",
     xlab = "Reflectance Value",
     main= "Distribution of MODIS reflectance values for each band\n Data not scaled")

Scale factor

Perform a radiometric calibration and change the format of the rasters from Digital Numbers to Reflectance. Plot the correspondent histograms and compare values range and distribution with those of the previous histograms.

# deal with nodata value --  -28672
all_modis_bands_pre_br <- all_modis_bands_pre_br * .0001

# plot the histograms of every band
hist(all_modis_bands_pre_br,
     xlab = "Reflectance Value",
     col = "springgreen",
     main = "Distribution of MODIS reflectance values for each band\n Scale factor applied")

No data values

Remove No data values. TIP: For better identification of these No data values, check the images in QGIS. Plot the correspondent histograms and compare with previous.

# deal with nodata value --  -100
all_modis_bands_pre_br[all_modis_bands_pre_br < 0 ] <- NA
 
# plot histogram

hist(all_modis_bands_pre_br,
     xlab = "Reflectance Value",
     col = "springgreen")
mtext("Distribution of reflectance values for each band", outer = TRUE, cex = 1.5)

Plot RGB

Plot a true color version of the MODIS brick.

plotRGB(all_modis_bands_pre_br, r=1, g=4, b=3, axes=TRUE, stretch="hist")

Remove clouds

Open cloud mask layer and mask out the values that correspond to clouds. TIP: check these images in QGIS as well in order to identify the clouds values which will be used for the mask.

# open cloud mask layer
cloud_mask <- raster("reflectance\\07_july_2016\\crop\\cloud_mask_july7_500m.tif")

# mask clouds
par(xpd = FALSE, mar = c(0, 0, 1, 5))

# create cloud & cloud shadow mask
cloud_mask[cloud_mask > 0] <- NA

# plot
plot(cloud_mask,
     main = "MODIS - Cloud mask layer")

Plot the masked data

plot(cloud_mask,
     main = "The Raster Mask",
     col = c("green"),
     legend = FALSE,
     axes = FALSE,
     box = FALSE)
# add legend to map
par(xpd = TRUE) # force legend to plot outside of the plot extent
legend(x = cloud_mask@extent@xmax, cloud_mask@extent@ymax,
       c("Not masked", "Masked"),
       fill = c("green", "white"),
       bty = "n")

Mask the clouds out of the raster brick AOI

In this step, you merge both cloud mask and MODIS brick, with the result of a MODIS set of bands without clouds.

all_modis_bands_pre_br_mask <- mask(all_modis_bands_pre_br, mask = cloud_mask)

# check the memory situation
inMemory(all_modis_bands_pre_br_mask)

## [1] TRUE

# check the classes or type of the MODIS brick after cloud masking
class(all_modis_bands_pre_br_mask)

## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"

Plot your cloud corrected image

# first turn all axes to the color white and turn off ticks
par(col.axis = "white", col.lab = "white", tck = 0)
# then plot the data
plotRGB(all_modis_bands_pre_br_mask,
        stretch = "lin",
        r = 1, g = 4, b = 3,
        main = "MODIS without clouds",
        axes = TRUE)

Modis Cloudmasking

Yrneh Ulloa, Severin Herzsprung

September 2020

Overview

Objectives

Data