A script to create a simple time series using MODIS satellite data with R

MODIS_timeseries.R

library(gdalio)
library(rnaturalearth)
library(dplyr)
library(sf)
library(lubridate)
library(magick)
library(terra)
source(system.file("raster_format/raster_format.codeR", 
                   package = "gdalio", mustWork = TRUE)) # helper functions

modis_time_series <- function(country, res, dates){
  # function to set gdalio grid with an country name.
  set_gdalio_from_country <- function(country, res){
    aoi <- ne_countries(scale = "medium", returnclass = "sf") %>%
      filter(sovereignt==country) %>%
      st_transform(crs=3857)
    
    bounds <- aoi %>%
      st_bbox() 
    grid0 <- list(extent = c(bounds$xmin, bounds$xmax, bounds$ymin, bounds$ymax), 
                  dimension = c((bounds$xmax-bounds$xmin)/res,
                                (bounds$ymax-bounds$ymin)/res), 
                  projection = st_crs(aoi)$wkt)
    gdalio_set_default_grid(grid0)
  }
  
  # function to retrieve a "best" MODIS composite for a given date.
  get_modis <- function(date){
    modis_imagery <- tempfile(fileext = ".xml")
    writeLines(sprintf('<GDAL_WMS>
    <Service name="TiledWMS">
    <ServerUrl>https://gibs.earthdata.nasa.gov/twms/epsg4326/best/twms.cgi?</ServerUrl>
    <TiledGroupName>MODIS Terra CorrectedReflectance TrueColor tileset</TiledGroupName>
    <Change key="${time}">%s</Change>
    </Service>
</GDAL_WMS>', date), modis_imagery)
    
    gdalio_terra(modis_imagery,bands = 1:3, resample = "cubicspline")
  }

  # set gdalio grid
  set_gdalio_from_country(country, res)
  
  dates %>%
    lapply(., get_modis) 
  
}

#map animation function
animate_maps <- function(map_list, dest, quality=1, n_frames=10, fps = 10){
  
  img <- magick::image_graph(dim(map_list[[1]])[2]*quality, dim(map_list[[1]])[1]*quality)
  lapply(map_list, terra::plotRGB)
  dev.off()
  animation <- magick::image_animate(magick::image_morph(img,  frames = n_frames), fps = fps, optimize = TRUE)
  
  magick::image_write(animation, normalizePath(file.path(dest), mustWork = FALSE))
  gc()
  return(animation)
}


ras_list <- modis_time_series(country='Italy', 
                  res=500, 
                  dates=seq(ymd('2021-07-01'), length.out = 7, by ='days'))

animate_maps(ras_list, 
             dest='exports/ItalySmall.gif',
             quality=0.25)