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'"

# check output maps in time series
t.rast.list 8day_avg_temp_2003
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)"
# 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 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)"

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
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 ))
doy_end=$(( $DOY + 5 ))
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

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

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 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 input=8day_agg \
title="8 day aggregated ts" \
description="8 day MODIS-like aggregated dataset"
# check 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
| input="8day_agg"
| title="8 day aggregated ts"
| description="8 day MODIS-like aggregated dataset"
t.rast.list input=8day_agg
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! 🙂


3 thoughts on “How to aggregate daily data into MODIS-like 8 day aggregation pattern?

  1. Thank you for this great workflow. This is really sophisticated.

    Maybe this can be performed using easier steps:
    t.rast.aggregate will use the granularity to create the start and end time of time intervals for the output STRDS. However, the where statement assures, that only maps of a specific date will be used for aggregation. In your first example the last time stamp ranges from 2003-12-27 00:00:00|2004-01-04 00:00:00 because of the granularity size.
    However, t.rast.aggregate assures that only maps from the current year are used for aggregation, if you use this where statement:

    where=”start_time >= ‘${1}-01-01’ and start_time < '${2}-01-01'"

    Hence, the only thing to do is to adjust the end time of the last interval for each year.

    You should also try the latest QGIS Grass-Data-Explorer plugin to visualize and analyse all your spatio-temporal Grass data. have a look at:

    Best regards


    • Hola huhabla!!

      Thanks much for your comment and especially, thanks for claryfication. Granularity combined with where clause is tricky. I thought aggregation of the last map was wrong given that the timestamp is wrong. So, we would just need to register those maps again setting start and end time properly, right? In this case, for just one of those maps we would use:

      t.register type=raster input=8day_avg_temp_2003 map=8day_avg_temp_2003_12_27 start="2003-12-27" end="2004-01-01" --o

      In the case we had already merged the three time series of the first example into one, to correct timestamps for last maps of every year, we would better use a file with mapname|start_time|end_time formatting. Then, the command would be:

      t.register type=raster input=8day_avg_temp_merged file=list_last_maps.txt --o

      Thanks again 🙂



  2. Pingback: How to aggregate daily data into MODIS-like 8 day aggregation pattern? –

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s