Temporal data processing/seasonal aggregation

Seasonal aggregation example

Introduction

As "seasons" or "astronomical seasons" are not yet implemented as a predefined granularity for aggregation, the workaround to get seasonal aggregation of a certain time series, is to create another time series with the desired granularity (i.e.: astronomical seasons) and apply or copy this granularity to aggregate data in our time series of interest.

Preparation of daily dataset

This is an example showing the use of t.rast.aggregate.ds to aggregate a daily time series using a source time series with seasonal granularity. We first create daily maps for more than one year (400 maps). We then create the corresponding spatio-temporal raster data set (strds) and register the maps into it.

# set the computational region and resolution
g.region -p n=-30 s=-50 e=-50 w=-70 res=1

# we do more than one year
for map in `seq 1 400` ; do
 # generate synthetic maps as a seeding
 r.mapcalc expression="daily_prec_${map} = 1"
 echo daily_prec_${map} >> map_list.txt
done

t.create type=strds temporaltype=absolute \
         output=precipitation_daily \
         title="Daily precipitation" \
         description="Test dataset with daily precipitation"

t.register -i type=raster input=precipitation_daily \
           file=map_list.txt start="2014-03-07" increment="1 days"

# check general information of the daily strds
t.info type=strds input=precipitation_daily

+-------------------- Space Time Raster Dataset -----------------------------+
|                                                                            |
+-------------------- Basic information -------------------------------------+
| Id: ........................ precipitation_daily@meteo
| Name: ...................... precipitation_daily
| Mapset: .................... meteo
| Creator: ................... veroandreo
| Temporal type: ............. absolute
| Creation time: ............. 2015-08-23 21:08:00.744761
| Modification time:.......... 2015-08-23 21:08:10.790248
| Semantic type:.............. mean
+-------------------- Absolute time -----------------------------------------+
| Start time:................. 2014-03-07 00:00:00
| End time:................... 2015-04-11 00:00:00
| Granularity:................ 1 day
| Temporal type of maps:...... interval
+-------------------- Spatial extent ----------------------------------------+
| North:...................... -30.0
| South:...................... -50.0
| East:.. .................... -50.0
| West:....................... -70.0
| Top:........................ 0.0
| Bottom:..................... 0.0
+-------------------- Metadata information ----------------------------------+
| Raster register table:...... raster_map_register_d14f30f1438840ac81cbfc0a554c1e69
| North-South resolution min:. 1.0
| North-South resolution max:. 1.0
| East-west resolution min:... 1.0
| East-west resolution max:... 1.0
| Minimum value min:.......... 1.0
| Minimum value max:.......... 1.0
| Maximum value min:.......... 1.0
| Maximum value max:.......... 1.0
| Aggregation type:........... None
| Number of registered maps:.. 400
|
| Title:
| Daily precipitation
| Description:
| Test dataset with daily precipitation
| Command history:
| # 2015-08-23 21:08:00 
| t.create type="strds" temporaltype="absolute"
|     output="precipitation_daily" title="Daily precipitation"
|     description="Test dataset with daily precipitation"
| # 2015-08-23 21:08:10 
| t.register -i type="raster"
|     input="precipitation_daily" file="map_list.txt" start="2014-03-07"
|     increment="1 days"
| 
+----------------------------------------------------------------------------+

Preparation of seasonal granularity dataset

Now, we'll create a second time series (might be raster, vector or raster3d) of seasonal granularity. We'll then take this granularity to aggregate our daily strds. Just as an example, we'll create a spatio-temporal vector data set (stvds).

# we create 4 maps of points (one per season)
for i in `seq 1 4` ; do 
 v.random output=points_${i} n=20
done

# create the stvds
t.create type=stvds temporaltype=absolute \
         output=points \
         title="Points" \
         description="Points for aggregation"

Now, we set the seasonal granularity for the vector maps. For this, we need to create a file with the respective map names and start and end date corresponding to astronomical seasons. We'll use this file as input for t.register.

points_1|2014-03-20|2014-06-21
points_2|2014-06-21|2014-09-23
points_3|2014-09-23|2014-12-21
points_4|2014-12-21|2015-03-20

t.register type=vector input=points file=seasonal_points.txt
t.info type=stvds input=points

+-------------------- Space Time Vector Dataset -----------------------------+
|                                                                            |
+-------------------- Basic information -------------------------------------+
| Id: ........................ points@meteo
| Name: ...................... points
| Mapset: .................... meteo
| Creator: ................... veroandreo
| Temporal type: ............. absolute
| Creation time: ............. 2015-08-23 21:10:11.462742
| Modification time:.......... 2015-08-23 21:49:59.554593
| Semantic type:.............. mean
+-------------------- Absolute time -----------------------------------------+
| Start time:................. 2014-03-20 00:00:00
| End time:................... 2015-03-20 00:00:00
| Granularity:................ 1 day
| Temporal type of maps:...... interval
+-------------------- Spatial extent ----------------------------------------+
| North:...................... -30.260367
| South:...................... -49.990657
| East:.. .................... -50.000035
| West:....................... -69.181185
| Top:........................ 0.0
| Bottom:..................... 0.0
+-------------------- Metadata information ----------------------------------+
| Vector register table:...... vector_map_register_727ad0a4157448ee8662d4345c248267
| Number of points ........... 80
| Number of lines ............ 0
| Number of boundaries ....... 0
| Number of centroids ........ 0
| Number of faces ............ 0
| Number of kernels .......... 0
| Number of primitives ....... 80
| Number of nodes ............ 0
| Number of areas ............ 0
| Number of islands .......... 0
| Number of holes ............ 0
| Number of volumes .......... 0
| Number of registered maps:.. 4
|
| Title:
| Points
| Description:
| Points for aggregation
| Command history:
| # 2015-08-23 21:10:11 
| t.create type="stvds" temporaltype="absolute"
|     output="points" title="Points" description="Points for aggregation"
| # 2015-08-23 21:49:59 
| t.register type="vector" input="points"
|     file="seasonal_points.txt"
| 
+----------------------------------------------------------------------------+

Aggregation of daily to seasonal maps

Finally, we aggregate our daily time series "precipitation daily" into seasonally accumulated precipitation with t.rast.aggregate.ds using method=sum.

# use everything between start date and end date
t.rast.aggregate.ds input=precipitation_daily \
                    output=precipitation_seasonal \
                    sample=points type=stvds \
                    base=prec_seasonal \
                    method=sum sampling=contains

# update information for the newly created seasonal time series
t.support input=precipitation_seasonal \
          title="Aggregated precipitation" \
          description="Aggregated seasonal precipitation dataset"

# check info                 
t.info type=strds input=precipitation_seasonal

+-------------------- Space Time Raster Dataset -----------------------------+
|                                                                            |
+-------------------- Basic information -------------------------------------+
| Id: ........................ precipitation_seasonal@meteo
| Name: ...................... precipitation_seasonal
| Mapset: .................... meteo
| Creator: ................... veroandreo
| Temporal type: ............. absolute
| Creation time: ............. 2015-08-23 21:50:26.134279
| Modification time:.......... 2015-08-23 21:50:34.365946
| Semantic type:.............. mean
+-------------------- Absolute time -----------------------------------------+
| Start time:................. 2014-03-20 00:00:00
| End time:................... 2015-03-20 00:00:00
| Granularity:................ 1 day
| Temporal type of maps:...... interval
+-------------------- Spatial extent ----------------------------------------+
| North:...................... -30.0
| South:...................... -50.0
| East:.. .................... -50.0
| West:....................... -70.0
| Top:........................ 0.0
| Bottom:..................... 0.0
+-------------------- Metadata information ----------------------------------+
| Raster register table:...... raster_map_register_be1ca044a59d49509e57a413e3781bc0
| North-South resolution min:. 1.0
| North-South resolution max:. 1.0
| East-west resolution min:... 1.0
| East-west resolution max:... 1.0
| Minimum value min:.......... 89.0
| Minimum value max:.......... 94.0
| Maximum value min:.......... 89.0
| Maximum value max:.......... 94.0
| Aggregation type:........... sum
| Number of registered maps:.. 4
|
| Title:
| Aggregated precipitation
| Description:
| Aggregated seasonal precipitation dataset
| Command history:
| # 2015-08-23 21:50:26 
| t.rast.aggregate.ds input="precipitation_daily"
|     output="precipitation_seasonal" sample="points" type="stvds"
|     base="prec_seasonal" method="sum" sampling="contains"
| # 2015-08-23 21:50:34 
| t.support input="precipitation_seasonal"
|     title="Aggregated precipitation"
|     description="Aggregated seasonal precipitation dataset"
| 
+----------------------------------------------------------------------------+

Verification

At the beginning we created 400 maps, remember??? We'll check now that only the maps "contained" in the period start_date-end_date of the stvds are used for the aggregation (the rest are left aside). That's why we use sampling=contains. Check t.rast.aggregate.ds for other method and sampling options. We use t.rast.series for this verification:

# check that the sum is 365
t.rast.series in=precipitation_seasonal method=sum out=suma_prec_seasonal
d.mon wx0
d.rast suma_prec_seasonal
# use GUI tools to check raster values