Using a shape file to download MODIS product data for country in R

Question

Is there any way that can be used to parse a shapefile of a country and download MODIS product data within that country using R?

I tried different approaches using the MODIStsp package (https://docs.ropensci.org/MODIStsp/) as well as the MODISTools package (https://docs.ropensci.org/MODISTools/articles/modistools-vignette.html) and they both only allow me to download MODIS product data for a defined site, but not a country.

Answer 1

Here's an example of how you might achieve this.

Firstly, download the MODIS data that you require, in this example I'm using MCD12Q1.006

begin_year and end_year are in the format: Year.Month.Days. shape_file is the shapefile you're using, presumably the extent of the shapefile is the country you're after. Though, I'm only going off by the minimal information you have provided.

library(MODIS)

tifs <- runGdal(product = "MCD12Q1", collection = "006", SDSstring = "01", 
                extent = shape_file %>% st_buffer(dist = 10000), 
                begin = begin_year, end = end_year, 
                outDirPath = "data", job = "modis",
                MODISserverOrder = "LPDAAC") %>% 
  pluck("MCD12Q1.006") %>% 
  unlist()

# rename tifs to have more descriptive names
new_names <- format(as.Date(names(tifs)), "%Y") %>% 
  sprintf("modis_mcd12q1_umd_%s.tif", .) %>% 
  file.path(dirname(tifs), .)
file.rename(tifs, new_names)

landcover <- list.files("data/modis", "^modis_mcd12q1_umd", 
                        full.names = TRUE) %>% 
  stack()
# label layers with year
landcover <- names(landcover) %>% 
  str_extract("(?<=modis_mcd12q1_umd_)[0-9]{4}") %>% 
  paste0("y", .) %>% 
  setNames(landcover, .)

Also, if you require a particular cell size, then you could follow this procedure to get a 5x5 modis cell size.

neighborhood_radius <- 5 * ceiling(max(res(landcover))) / 2

agg_factor <- round(2 * neighborhood_radius / res(landcover))
r <- raster(landcover) %>% 
  aggregate(agg_factor) 
r <- shape_file %>% 
  st_transform(crs = projection(r)) %>% 
  rasterize(r, field = 1) %>% 
  # remove any empty cells at edges
  trim()

Answer 2

Here's an example using MODISTools to automate downloading the correct tiles for the country.

First let's generate a polygon of a country to demonstrate (using Luxembourg as an example):

library(maptools)
library(sf)
data(wrld_simpl)
world = st_as_sf(wrld_simpl)
lux = world[world$NAME=='Luxembourg',]

Now we find the location (centroid) and size of the country:

#find centroid of polygon in long-lat decimal degrees
lux.cent = st_centroid(lux)

#find width and height of country in km
lux.proj = st_transform(lux, 
           "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=km +no_defs")
lux.km_lr = diff(st_bbox(lux.proj)[c(1,3)])
lux.km_ab = diff(st_bbox(lux.proj)[c(2,4)])

Using this info, we can download the correct Modis data (using leaf-area index, lai, as an example):

#download the MODIS tiles for the area we defined
library(MODISTools)
lux_lai <- mt_subset(product = "MOD15A2H",
                          lat = lux.cent$LAT, lon =  lux.cent$LON,
                          band = "Lai_500m",
                          start = "2004-01-01", end = "2004-01-01",
                          km_lr = lux.km_lr, km_ab = lux.km_ab,
                          site_name = "Luxembourg",
                          internal = TRUE, progress = TRUE)

# convert to a spatial raster
lux.rast = mt_to_raster(df = lux_lai, reproject = TRUE)
lux.rast = raster::mask(lux.rast, lux)
plot(lux.rast)
plot(st_geometry(lux),add=T)