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

Last updated on May 24, 2020 6 min read

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"

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

# 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" | 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. 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"

# register maps 
t.register type=raster input=8day_ts \ 
file=list_map_start_end_time.txt

# check 
t.info input=8day_ts 
t.rast.list input=8day_ts

Finally, we 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! :)

Veronica Andreo

Researcher and Lecturer

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