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 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 - new version of the 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(out_qc_test1_mandatory < 2 && out_qc_test1_lst_error == 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:

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

Enjoy! 🙂

and Happy New Year!

GRASS GIS workshop and talks at FOSS4G 2016 in Bonn – You cannot miss them!

Next Sunday, August 21, and for one week, an historic event will take place in the city of Bonn – Germany … the international FOSS4G Conference!! This conference brings together people from all over the world that love to have FUN and FOSS 🙂

GRASS GIS, being one of the first free and open source GIS, won’t miss the appointment nor the chance to meet friends and share. You shouldn’t either miss all the new goodies that will be presented in the dedicated GRASS GIS workshop and talks. Here are the links:

Workshop:

Unleash the power of GRASS GIS 7

When? Tuesday, August 23, from 9:00 to 13:00

Trainers: Markus Neteler (mundialis GmbH & Co KG), Luca Delucchi (Fondazione Edmund Mach), Martin Landa (Czech Technical University in Prague)

More info: http://2016.foss4g.org/ws27.html


Talks:

#289: A complete toolchain for object-based image analysis with GRASS GIS

Moritz Lennert (http://2016.foss4g.org/talks.html#289)

#290: OpenDEM Generator: combining open access Digital Elevation Models into a homogenized DEM

Luca Delucchi and Markus Neteler (http://2016.foss4g.org/talks.html#290)

#318: Processing Copernicus Sentinel data with GRASS GIS

Markus Neteler and Carmen Tawalika (http://2016.foss4g.org/talks.html#318)

#533: Building applications with FOSS4G bricks: two examples of the use of GRASS GIS modules as a high-level “language”‘ for the analyses of continuous space data in economic geography

Moritz Lennert (http://2016.foss4g.org/talks.html#533)

and so much more…!!!!

Just come… and build bridges!

foss4g-logo

Do you need to process time series??

Then, don’t miss “TGRASS: time series processing in GRASS GIS“, a 4-hours workshop during FOSS4G-ar, the conference on free and open source software for geospatial that will take place next April in Buenos Aires, Argentina.

Screenshot from 2015-06-15 17:47:45

We will learn the basics of time series management, processing and analysis in GRASS GIS by means of a sample dataset of raster maps. We will use a subset of the more than 40 temporal modules available. TGRASS functionalities include tools for management and administration of spatio-temporal data sets and time stamped maps, temporal aggregation and accumulation, obtaining of basic and extended statistics, visualization and animation, import/export, data extraction through vector maps, among others.

Don’t miss this chance to learn such a great tool! Book now! And register yourself for FOSS4G-ar!

See you there!

foss4g-ar-1_black

 

 

How to aggregate daily data into MODIS-like 8 day aggregation pattern?

Several MODIS products come in “8-day” compositions. Suppose we have daily data and we want to aggregate it with this MODIS-like granularity. How can we achieve that?

First idea would be to use the great t.rast.aggregate module in GRASS GIS with granularity=”8 days”. But this presents a problem. It does not re-start every year, but aggregates the whole series with 8-day granularity, not considering the year.

Meanwhile, MODIS 8-day products re-start the aggregation every January 1st. Therefore, second idea could be to use t.rast.aggregate with granularity=”8 days” but looping over years, and then merge the resulting strds with t.merge. However, granularity seems to overrule the where clause. This is evident in the output of  t.rast.list (i.e.: last map granularity covers 3 days from the next year) in the following example:

for YEAR in "2003 2004" "2004 2005" "2005 2006" ; do
  set -- $YEAR ; echo $1 $2
  t.rast.aggregate -s input=daily_temp output=8day_avg_temp_${1} \
  basename=8day_avg_temp method=average granularity="8 days" \
  where="start_time >= '${1}-01-01' and end_time < '${2}-01-01'"
done

# check output maps in time series
t.rast.list 8day_avg_temp_2003
name|mapset|start_time|end_time
8day_avg_temp_2003_01_01|climate|2003-01-01 00:00:00|2003-01-09 00:00:00
8day_avg_temp_2003_01_09|climate|2003-01-09 00:00:00|2003-01-17 00:00:00
...
8day_avg_temp_2003_12_19|climate|2003-12-19 00:00:00|2003-12-27 00:00:00
8day_avg_temp_2003_12_27|climate|2003-12-27 00:00:00|2004-01-04 00:00:00

So, we need a different solution. If we already have a registered time series of MODIS 8-day products, we can just use it to copy its aggregation pattern with t.rast.aggregate.ds. But what if we don’t? How do we create a MODIS-like 8 day granularity to use it as reference to aggregate our daily data??

One solution could be to create a time series (could be raster or vector) with that particular pattern of aggregation and then, use it to aggregate our daily time series. In this example, I’ll show you how to create a raster time series, from now on strds, which stands for “spatio-temporal raster dataset”. But, in order to save space in disk we’ll set a region of only one pixel in the center of the current region

eval `g.region -g`
eval `g.region -cg`
g.region w=$center_easting s=$center_northing \
e=`echo "$center_easting + $ewres" | bc` \
n=`echo "$center_northing + $nsres" | bc` -p

Now, we’ll create daily maps and register them as our daily time series, which will be our input time series (see below).

# create daily maps
for i in `seq -w 1 740` ; do
  r.mapcalc -s expression="daily_ts_${i} = rand(0.0,40.0)"
done
# create daily ts
t.create type=strds temporaltype=absolute \
output=daily_ts title="Daily time series" \
description="Test STRDS with 1 day granularity" --o
# register daily maps
t.register -i type=raster input=daily_ts \
maps=`g.list type=raster pattern=daily_ts_* separator=comma` \
start="2000-01-01" increment="1 days" --o
# check info
t.info input=daily_ts
t.rast.list input=daily_ts

And then we create the 8-day MODIS-like maps and register them as time series. This will be our sample strds, the data set from which we’ll copy the aggregation pattern. We use year and DOY (day of year) to name maps. These will then be converted into calendar dates to register maps as a time series with absolute time.

# create maps every 8 days
for year in `seq 2000 2001` ; do
  for doy in `seq -w 1 8 365` ; do
    r.mapcalc -s expression="8day_${year}_${doy} = rand(0.0,40.0)"
  done
done

After we created maps, we need to assign timestamps to them representing 8-days intervals for all maps, except for the last map in each year, since the aggregation starts all over again, every January 1st.

# mapnames list
g.list type=raster pattern=8day_20??_* > names_list

From de name of each map, we take year and doy, and convert it to a YYYY-MM-DD date for start and end (we need time intervals), considering also if the year is a leap year or not.

# create file with mapnames, start date and end date
for NAME in `cat names_list` ; do
# Parse
YEAR=`echo $NAME | cut -d'_' -f2`
DOY=`echo $NAME | cut -d'_' -f3`
# convert YYYY_DOY to YYYY-MM-DD
# BEWARE: leading zeros make bash assume the number
# is in base 8 system, not base 10!
DOY=`echo "$DOY&quot; | sed 's/^0*//'`
# list with mapname and both start and end time
doy_end=0
if [ $DOY -le "353" ] ; then
doy_end=$(( $DOY + 8 ))
elif [ $DOY -eq "361" ] ; then
if [ $[$YEAR % 4] -eq 0 ] && [ $[$YEAR % 100] -ne 0 ] || [ $[$YEAR % 400] -eq 0 ] ; then
doy_end=$(( $DOY + 6 ))
else
doy_end=$(( $DOY + 5 ))
fi
fi
DATE_START=`date -d "${YEAR}-01-01 +$(( ${DOY} - 1 ))days" +%Y-%m-%d`
DATE_END=`date -d "${YEAR}-01-01 +$(( ${doy_end} -1 ))days" +%Y-%m-%d`
# text file with mapnames, start date and end date
echo "$NAME|$DATE_START|$DATE_END" >> list_map_start_end_time.txt
done

We check the list. Intervals are left open, so end time of one map is the start time of the next. And we also check that the last map of each year ends as expected.

cat list_map_start_end_time.txt
8day_2000_001|2000-01-01|2000-01-09
8day_2000_009|2000-01-09|2000-01-17
...
8day_2000_353|2000-12-18|2000-12-26
8day_2000_361|2000-12-26|2001-01-01
8day_2001_001|2001-01-01|2001-01-09
8day_2001_009|2001-01-09|2001-01-17
...
8day_2001_345|2001-12-11|2001-12-19
8day_2001_353|2001-12-19|2001-12-27
8day_2001_361|2001-12-27|2002-01-01

We are ready now to create our 8-day MODIS-like strds.

# create strds
t.create type=strds temporaltype=absolute \
output=8day_ts title="8 day time series" \
description="STRDS with MODIS like 8 day aggregation" --o
# register maps
t.register type=raster input=8day_ts \
file=list_map_start_end_time.txt --o
# check
t.info input=8day_ts
t.rast.list input=8day_ts

And finally, copy its aggregation to our daily time series.

# copy the aggregation
t.rast.aggregate.ds -s input=daily_ts sample=8day_ts \
output=8day_agg basename=8day_agg \
method=average sampling=contains
# add metadata
t.support input=8day_agg \
title="8 day aggregated ts" \
description="8 day MODIS-like aggregated dataset"
# check
t.info input=8day_agg
+-------------------- Space Time Raster Dataset -----------------------------+
|
+-------------------- Basic information -------------------------------------+
| Id: ........................ 8day_agg@pruebas
| Name: ...................... 8day_agg
| Mapset: .................... pruebas
| Creator: ................... veroandreo
| Temporal type: ............. absolute
| Creation time: ............. 2016-01-09 21:44:54.074235
| Modification time:.......... 2016-01-10 16:39:06.253565
| Semantic type:.............. mean
+-------------------- Absolute time -----------------------------------------+
| Start time:................. 2000-01-01 00:00:00
| End time:................... 2002-01-01 00:00:00
| Granularity:................ 1 day
| Temporal type of maps:...... interval
+-------------------- Spatial extent ----------------------------------------+
| North:...................... -33.0
| South:...................... -33.5
| East:.. .................... -57.0
| West:....................... -57.5
| Top:........................ 0.0
| Bottom:..................... 0.0
+-------------------- Metadata information ----------------------------------+
| Raster register table:...... raster_map_register_115633bc4e0a4195b6cb4fdcab505522
| North-South resolution min:. 0.5
| North-South resolution max:. 0.5
| East-west resolution min:... 0.5
| East-west resolution max:... 0.5
| Minimum value min:.......... 2.935881
| Minimum value max:.......... 28.388835
| Maximum value min:.......... 2.935881
| Maximum value max:.......... 28.388835
| Aggregation type:........... average
| Number of registered maps:.. 92
|
| Title:
| 8 day aggregated ts
| Description:
| 8 day MODIS-like aggregated dataset
| Command history:
| # 2016-01-09 21:44:54
| t.rast.aggregate.ds -s input="daily_ts"
| sample="8day_ts" output="8day_agg" basename="8day_agg"
| method="average" sampling="contains"
| # 2016-01-10 16:39:06
| t.support input="8day_agg"
| title="8 day aggregated ts"
| description="8 day MODIS-like aggregated dataset"
|
+----------------------------------------------------------------------------+
t.rast.list input=8day_agg
name|mapset|start_time|end_time
8day_agg_2000_01_01|pruebas|2000-01-01 00:00:00|2000-01-09 00:00:00
8day_agg_2000_01_09|pruebas|2000-01-09 00:00:00|2000-01-17 00:00:00
8day_agg_2000_01_17|pruebas|2000-01-17 00:00:00|2000-01-25 00:00:00
...
8day_agg_2000_12_18|pruebas|2000-12-18 00:00:00|2000-12-26 00:00:00
8day_agg_2000_12_26|pruebas|2000-12-26 00:00:00|2001-01-01 00:00:00
8day_agg_2001_01_01|pruebas|2001-01-01 00:00:00|2001-01-09 00:00:00
...
8day_agg_2001_12_11|pruebas|2001-12-11 00:00:00|2001-12-19 00:00:00
8day_agg_2001_12_19|pruebas|2001-12-19 00:00:00|2001-12-27 00:00:00
8day_agg_2001_12_27|pruebas|2001-12-27 00:00:00|2002-01-01 00:00:00

Note that some maps from our daily strds remain not aggregated into the 8 day new time series, that’s because we used “sampling=contains”, which assures that only raster maps that are temporally during the time intervals of the strds are considered for computation.

Enjoy! 🙂