Using MODIS in R

To continue on for Part 2, make sure you have completed all prior steps explained in Part 1.

Section 1: Create Folder and Set Directory

Step 1:

Before we start go to your desktop and create a new folder names MODIS_Tutorial

Step 2: Open a new session of R studio

After opening R studio your window should look similar to this

Step 3: Open a new R script

Once you have opened R studio, open a new empty R script by clicking the highlighted icon in the image below, and selecting R script.

Step 4: Save script

Save empty script in the MODIS_Tutorial folder we created in Step 1.

Step 5: Set working directory

Click on the tab Session on the task bar, then click on Set Working Directory, and then click Choose Directory…. Chose the MODIS_Tutorial folder.

Step 6: Create a new folder directly from R

Create another folder for the MODIS data inside the MODIS_Tutorial folder, and change the working directory to your new folder. Type or Copy/Paste the following lines into your R script. The function dir.create() will create a folder within your working directory, and the function setwd() will set your working directory to the new created folder.

dir.create("MODIS_data")
setwd("MODIS_data")

Section 2: Installing and Loading Packages for MODIStsp

Step 1:

MODIStsp requires R v >= 3.2.1 and GDAL (Geospatial Data Abstraction Library) v >= 1.11.1 with support for HDF4 raster format to be installed in your system. The gWidgets API is intended to be a cross-toolkit API for working with GUI objects.

# rgdal provides bindings to the 'Geospatial' Data Abstraction library ('GDAL')
install.packages("rgdal") 

install.packages("gWidgetsRGtk2")

# MODIStsp performs several preprocessing steps (e.g., download, mosaicing, reprojection and resize) on MODIS data available within a given time period
install.packages("MODIStsp") 

## To install the development version (containing the latest improvements and bug fixes) from github:
install.packages("devtools")

library(rgdal)
library(gWidgetsRGtk2)
library(MODIStsp)
library(devtools)

If you get any errors when loading the libraries look for potential fixes in the common problem section, here. Note, fixes can be as simple as installing a new library.

Section 3: Run MODIS

Step 1: Run MODIStsp()

To deploy the modis application run the following code:

MODIStsp()

Please note:

* MODIStsp takes awhile to load.

* When prompted, you will need to enter the number one to select the option that starts with "Yes"

* Afterwards, you need to wait a moment for the application to check for a GDAL installation.

If MODIStsp() throws and error, look in the common issues page for possible fixes. here.

Step 2:

Step 3:

Step 4:

NOTE: if your GDAL program is not up to date or you have not downloaded it MODIS wont run. Look in the here for how to download GDAL in your computer.

Section 4: Understanding How to Use MODIS

If everything has ran successfully a window like the one below should have popped up in your screen.

Step 1:

Step 2:

Step 3:

Step 4:

Step 5:

Step 6:

Step 7:

Step 8:

Step 9:

Step 10:

You can enter the number range of the tiles you will like to use by viewing the tiles displayed on Google Earth. Follow this link to access tiles overlaid on an Earth image in Google Earth Modis Tile. Below is an example of how the image will appear in Google Earth:

Note, this file will only open if you have Google Earth downloaded in your computer. If you do not have Google Earth downloaded you can download it here.

Step 11:

Step 12:

Step 13:

Change projection to the following string:

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Step 14:

Step 15:

Step 16:

Step 17:

Step 18:

Step 19:

Step 20:

Step 21:

The time it would take for your data to download will depend on the time range and tile range selected.

Part 3: Working with Raster files

Step 1: Download Modis Data

Now, we will provide an example. Data for this example has the following characteristics:

Category: Ecosystem Variables - Vegetation Indices
Product: Vegetation_Index_16Days_250m 
Platform (aka satellite): Terra
Original MODIS layer: 16 day EVI average

We will need to download two different set of date ranges to compare the vegetation loss, this means you’ll have to process the request twice, and save each year in a separate folder.

First date range: 2003-01-01 to 2003-01-31
Second date range: 2019-01-01 to 2019-01-31

And with the MODIS tiles:

Horizontal: 10 to 10
Verical: 9 to 9

Remember to change projection to User Defined and write down the following string:

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

Set the appropriate folders for data and HDF files, and lets get started!

Step 2: R Code

Step 2.1: Load Necessary Packages

Run the following libraries, we will need all of them to make the analysis and visualization of our data. If you have never used this libraries before you might need to install them first by using the function install.packages("name of library").

library(raster)
library(rasterVis)

Step 2.2: Set Working Directory

Now, lets load in our MODIS data into R. The first step to do this is to set your working directory to the location of your files using the function setwd(). For example, it can look something like this:

setwd("C:/Users/Desktop/Modis_Tutorial/MODIS_data/2003")

Step 2.3: Load Data

To load all the files in our working directory we need to use the function list.files(). This function reads in all the files in a given folder with an specified pattern (in this case .tif).

files_2003 <- list.files(getwd(), 
                         pattern = ".tif$",
                         full.names = TRUE, 
                         recursive =  TRUE)

Step 2.4: Stack raster and Calculate Mean

To stack a list of rasters we use the function stack() in the package raster. To calculate the mean we use the function calc() and as the function to perform in raster stack is mean. Note, MODIS does not allow for decimals, so if you want only values from 0-1 we dived need to divided by 10000.

files_2003_stack <- stack(files_2003)
mean_2003 <- calc(files_2003_stack, fun= mean, na.rm=TRUE)/10000

As an additional thing, you can always save your new created raster using the function writeRaster(). For example:

writeRaster(mean_2003, format="ascii", filename= "mean_2003.ascii", overwrite=TRUE)

Step 2.5: Load The Rest of the Data (i.e., 2019)

Repeat steps 2.2 - 2.4, creating new variable names and setting a new working directory, just as follows:

setwd("C:/Users/Desktop/Modis_Tutorial/MODIS_data/2019")

files_2019 <- list.files(getwd(), 
                         pattern = ".tif$",
                         full.names = TRUE, 
                         recursive =  TRUE)
files_2019_stack <- stack(files_2019)
mean_2019 <- calc(files_2019_stack, fun= mean, na.rm=TRUE)/10000

Step 2.6: Calculate the Vegetation Loss

To calculate the vegetation loss within the years of study, we simply subtract the both mean rasters created in the previous steps, and name the resulting raster something appropriate like “diff.evi”.

diff.evi <- mean_2019-mean_2003

Step 2.7: Reclassify the data

Since we only want to know if there was a loss in vegetation, we need to make all positive numbers equal NA. To do this we use the function reclassify(). Inside this function there are two arguments, the first one represents the raster to reclassify (diff.evi), and the second one is the values within the raster that need to be reclassify. In addition, c(0, Inf, NA) read as: to all the values from 0 to Inf give them a value of NA.

diff.evi_reclass <- reclassify(diff.evi, c(0, Inf,NA))

Step 2.8: Visualizing Data

To visualize our new reclassified vegetation loss raster (aka diff.evi_reclass) we can use the function levelplot() in the package rasterVis.

levelplot(diff.evi_reclass)