In ptaconet/opendapr: Fast and Efficient Access to MODIS Earth Observation Data

NOT_CRAN <- !identical(tolower(Sys.getenv("NOT_CRAN")), "true") # vignette will not be executed when tested on the cran
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  purl = NOT_CRAN
  )

library(modisfast)
library(sf)
library(terra)
library(mapview)

In this example we want to download, import and plot over the 3500 km^2^ wide region of interest (mapped below) :

• a 1 year-long time series of monthly NDVI from MODIS Terra Vegetation indices (spatial resolution : 1 km ; temporal resolution : 1 month),
• the same 1 year-long times series of monthly precipitations from Global Precipitation Measurement (GPM) (spatial resolution : 1° ; temporal resolution : 1 month)

roi <- st_as_sf(data.frame(geom="POLYGON ((-5.82 9.54, -5.42 9.55, -5.41 8.84, -5.81 8.84, -5.82 9.54))"), wkt="geom", crs = 4326)
mapview::mapview(roi, legend = FALSE)

Check which collections are available

Before starting, let's identify our collections of interest among the collections that are available for download, using the function mf_list_collections() :

collections <- mf_list_collections()
print(str(collections))

So our collections of interest are :

• MOD13A3.061 : MODIS/Terra Vegetation Indices Monthly L3 Global 1 km SIN Grid
• GPM_3IMERGM.07 : GPM IMERG Final Precipitation L3 1 month 0.1 degree x 0.1 degree V07

Setup the region and time range of interest

First we prepare the script : define the ROI, the time frame, and the folder were the data will be downloaded :

# Set ROI and time range of interest
roi <- st_as_sf(data.frame(id="Korhogo", geom="POLYGON ((-5.82 9.54, -5.42 9.55, -5.41 8.84, -5.81 8.84, -5.82 9.54))"), wkt="geom", crs = 4326)
time_range <- as.Date(c("2017-01-01","2017-12-31"))

Login to EOSDIS Earthdata

And we login to EOSDIS Earthdata with our credentials. To create an account go to : https://urs.earthdata.nasa.gov/.

username <- Sys.getenv("earthdata_un")
password <- Sys.getenv("earthdata_pw")
log <- mf_login(credentials = c(username,password),verbose = "quiet")

# Login to Earthdata servers with username and password. 
username <- "earthdata_un"
password <- "earthdata_pw"
log <- mf_login(credentials = c(username,password))

username <- Sys.getenv("earthdata_un")
password <- Sys.getenv("earthdata_pw")
log <- mf_login(credentials = c(username,password),verbose = "quiet")

Get the URLs of the data to download {#get-url}

With the function mf_get_url(), we get the https URLs for our collections of interest (MOD11A1.061 and GPM_3IMERGDF.07) given our ROI and time frame.

:arrow_right: For each collection, we can select a subset only of variables (sometimes called "dimensions", or "bands") that are of interest for us, through the parameter variables. By default, this parameter is set to NULL, which means that all the available variables for the specified collection are included. To get the variables available for a given collection along with information for each of them (description, etc.), use the function mf_list_variables().

For more information on the variables in the context of OPeNDAP, see the OPeNDAP Data model documentation.

## Get the URLs of MODIS Vegetation indexes
urls_mod11a1 <- mf_get_url(
  collection = "MOD13A3.061",
  variables = c("_1_km_monthly_NDVI","_1_km_monthly_EVI","_1_km_monthly_VI_Quality"),  # get the available variables with : mf_list_variables("MOD13A3.061"), or set to NULL (default) to download all the variables
  roi = roi,
  time_range = time_range )

## Get the URLs of GPM daily
urls_gpm <- mf_get_url(
  collection = "GPM_3IMERGM.07",
  variables = c("precipitation"),  # get the available variables with : mf_list_variables("GPM_3IMERGM.07")
  roi = roi,
  time_range = time_range )

nrow(urls_mod11a1)
urls_mod11a1

nrow(urls_gpm)
head(urls_gpm,3)

Download the data {#download}

Now we download the data with the function mf_download_data().

The destination file for each dataset is specified in the column destfile of the dataframes urls_mod11a1 and urls_gpm. A destination file is specified by default but it can of course be modified.

Setting the argument parallel to TRUE will parallelize - therefore fasten - the download in case numerous datasets need to be downloaded.

# Hint : to fasten the download in case of multiple collections, bind them in a single data.frame (as shown below) before executing the function 'mf_download_data()' :
df_to_dl <- rbind(urls_mod11a1, urls_gpm)

# Then download the data (by default, in temporary directory, but the target folder can be set with the argument 'path'
res_dl <- mf_download_data(df_to_dl = df_to_dl, 
                           parallel = TRUE)

print(str(res_dl))

Import the data in R {#import}

Various packages and related classes can be used to read the data downloaded through OPeNDAP. If terra (successor of raster) is surely the most famous class for raster objects, many packages facilitate the use of spatiotemporal data cubes in formats such as those proposed through modisfast (e.g. NetCDF). For instance, MODIS or VIIRS products can be imported as a stars object from the excellent stars package for data cubes manipulation. All the data can also be imported as ncdf4 objects using e.g. the ncdf4 package.

:warning: :warning: :warning: :warning:

Important warning regarding the use of the data after download through modisfast {#warning-import}

Care must be taken when importing data that was downloaded through the OPeNDAP data providers servers. Depending on the collection, some "issues" were raised. These issues are independent from modisfast : they result most of time of a kind of lack of full implementation of the OPeNDAP framework by the data providers. These issues are :

• for MODIS and VIIRS collections : CRS has to be provided
• for GPM collections : CRS has to be provided + data have to be flipped

End warning

These issues can easily be dealt at the import phase in R using the function mf_import_data(). This function enables to open the data in R as a SpatRaster object (from the package terra) :

# import MOD11A1.061
# argument 'path' is the folder where the data are stored (here, in a temporary directory)
modis_ts <- mf_import_data(path = dirname(list.files(path = tempdir(), pattern = "MOD13A3.061", recursive = TRUE, full.names = TRUE)),
                           collection = "MOD13A3.061",
                           proj_epsg = 4326)

# import GPM
gpm_ts <- mf_import_data(path = unique(dirname(list.files(path = tempdir(), pattern = "3IMERG", recursive = TRUE, full.names = TRUE))),
                         collection = "GPM_3IMERGM.07",
                         proj_epsg = 4326)

modis_ts
gpm_ts

Plot the data

Let's finally plot the data !

(Note that only the first 12 dates are plotted here)

# Land surface temperature
terra::plot(modis_ts["_1_km_monthly_NDVI"])
# Precipitation
terra::plot(gpm_ts["precipitation"])

ptaconet/opendapr documentation built on Nov. 20, 2024, 10:04 p.m.