pymodis, GDAL and GRASS GIS: A FOSS4G example to build time series of MODIS data

Last updated on May 24, 2020 3 min read

Recently, the National Aeronautics and Space Administration (NASA) has re-processed all the MODIS archive and made available the sixth version for all land products. In this post, I share a script that combines different FOSS4G tools to download, mosaic, re-project, subset, import, apply quality assessment flags and build a time series of MODIS LST daily products. First, we use the pymodis tool  modis_download  to automatically download 4 tiles of the MOD11A1 product for the period 2010-2016.  Note that you need to register yourself at https://urs.earthdata.nasa.gov/users/new  and enable the applications at https://urs.earthdata.nasa.gov/home  to be able to download MODIS land products.

# Download MOD11A1.006 - daily land surface temperature (1km) 
modis_download.py -P your_password -U your_user -s MOLT -p MOD11A1.006 -t h19v04,h19v05,h20v04,h20v05 -f 2010-01-01 -e 2016-12-31 .

Then, we use modis_mosaic to build daily mosaics for LST Day and QC Day bands by means of VRT files.

modis_mosaic.py -v -s "1 1 0 0 0 0 0 0 0 0 0 0" list_of_files.txt

After that, we re-project the mosaics to EPSG 3035 and convert them to GTiff format with modis_convert.

# cut the first part of the filename containing date (year+doy): A2015365 
LIST_DATES=`ls *.vrt | cut -d'_' -f1` 

for DAY in $LIST_DATES ; do 
  # convert to tif and project 
  modis_convert.py -g 1000 -e 3035 -v -s "(1)" -o ${DAY}_LST_Day_1km_mosaic ${DAY}_None_LST_Day_1km.vrt 
done

Once we have the tif files, we perform a spatial subset for our area of interest. For this purpose, we combine the GRASS command g.region and the possibility to use eval command with gdal_translate  utility from GDAL/OGR library.

# set region 
g.region -p region=your_study_area 
eval `g.region -g` 

for DAY in $LIST_DATES ; do 
  # spatial subseting 
  gdal_translate -of GTiff -projwin $w $n $e $s ${DAY}_LST_Day_1km_mosaic.tif ${DAY}_LST_Day_1km_subset.tif 
done

Once the subset is done, we import the resulting maps into our GRASS database with r.in.gdal, and we remove all VRT, HDF and TIF files.

for DAY in $LIST_DATES ; do 
  # import into GRASS 
  r.in.gdal input=${DAY}_LST_Day_1km_subset.tif output=${DAY}_LST_Day_1km 
  
  # remove vrt, tiff and hdf files 
  rm -f ${DAY}_None_LST_Day_1km.vrt ${DAY}_None_LST_Day_1km_mosaic.tif ${DAY}_LST_Day_1km_subset.tif 
done

We then use i.modis.qc to get QC flags for LST and apply them as in Metz et al. (2014) to keep the highest quality data.

for map in `g.list type=raster pattern=*_QC_*` ; do 
  i.modis.qc input=${map} output=${map}_mandatory_qa productname=mod11A1 qcname=mandatory_qa_11A1 
  i.modis.qc input=${map} output=${map}_lst_error_qa productname=mod11A1 qcname=lst_error_11A1 
done

list_lst=`g.list rast pat=*LST_Day_1km` 

for m in ${list_lst} ; do 
  # cut common part of filenames 
  i=`echo $m | cut -c 1-23` 
  r.mapcalc --o expression="${m} = if(${i}_mandatory_qa < 2 && ${i}_lst_error_qa == 0, ${m}, null())" 
done

With all the former steps done, we can build a daily time series for LST data…

t.create type=strds temporaltype=absolute \
  output=LST_Day_daily \
  title="LST Day 1km MOD11A1.006" \
  description="Daily LST Day 1km MOD11A1.006, 2010-2016" 

g.list rast pattern=*LST_Day output=filelist_lst.txt 

t.register -i input=LST_${lst}_greece_daily \
  file=filelist_lst.txt \
  start="2010-01-01" \
  increment="1 day"

and you are ready to play! You may, for example, use the GRASS Add-ons r.hants or r.series.lwr to reconstruct/fill the gaps in the time series and then, you have a whole set of temporal modules at your fingertips to perform different tasks. You can check the Temporal Data Processing wiki for several examples.

Extra: The extra beauty of all this, is that it is possible to run all these steps as a script with the new --exec  interface of GRASS GIS. For example, if you name your script lst_processing.sh (do not forget to write the shebang in the beginning of the file: #!/bin/bash), then, you can run it in a non-interactive GRASS session as:

grass78 /path/to/grassdata/location/mapset --exec sh lst_processing.sh

Enjoy! :) and Happy New Year!

gdal GRASS GIS MODIS pymodis time series
Veronica Andreo
Veronica Andreo
Researcher and Lecturer

My research interests include remote sensing and GIS applied to public health problems.