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

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 and enable the applications at to be able to download MODIS Land products.

# Download MOD11A1.006 - new version of the daily land surface temperature (1km) -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. -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 -g 1000 -e 3035 -v -s "(1)" -o ${DAY}_LST_Day_1km_mosaic ${DAY}_None_LST_Day_1km.vrt

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

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

for DAY in $LIST_DATES ; do
# import into GRASS 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

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

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(out_qc_test1_mandatory < 2 && out_qc_test1_lst_error == 0, ${m}, null())"

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 (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:

grass72 /path/to/grassdata/location/mapset --exec sh

Enjoy! 🙂

and Happy New Year!