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! 🙂