[GRASS-user] Temporal framework: calculating annual 5-day extremes

GRASS GRASS user
1

Hi,

I have a data set containing multiple years of daily raster layers and would
like to aggregate and output annual raster layers of 5-day extremes
(maxima).

Essentially, for each grid cell, I need to calculate rolling 5-day sums of
each year and then find the annual max of the latter sums, and output as a
series of raster layers of annual 5-day extremes.

However, I'm trying to work out the best way in GRASS of doing this. Overall
t.rast.aggregate comes closest to the type of functionality needed. I've
also looked at t.rast3d.mapcalc, but unsure if calculating a rolling sum is
possible.

I'd be grateful for any suggestions on how I might approach this.

Thanks,

Richard

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014.html
Sent from the Grass - Users mailing list archive at Nabble.com.

2

r.series for the computation and g.mlist for selecting the rasters is how I have used in the past.
Sitansu

···

On Thu, Apr 6, 2017 at 10:01 AM, RichardCooper <richtcooper@hotmail.com> wrote:

Hi,

I have a data set containing multiple years of daily raster layers and would
like to aggregate and output annual raster layers of 5-day extremes
(maxima).

Essentially, for each grid cell, I need to calculate rolling 5-day sums of
each year and then find the annual max of the latter sums, and output as a
series of raster layers of annual 5-day extremes.

However, I’m trying to work out the best way in GRASS of doing this. Overall
t.rast.aggregate comes closest to the type of functionality needed. I’ve
also looked at t.rast3d.mapcalc, but unsure if calculating a rolling sum is
possible.

I’d be grateful for any suggestions on how I might approach this.

Thanks,

Richard

–
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014.html
Sent from the Grass - Users mailing list archive at Nabble.com.

grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

3

Hi,
i am not sure if i understand the problem correctly. However, you can use
t.rast.accumulate [1] to create the rolling sum for an arbitrary interval
(5 days, one year?) and then use t.rast.aggregate
to compute the yearly maximum value time series based on the time series
output of t.rast.accumulate.

Or you can use t.rast.aggregate two time, first to compute the 5 day
aggregation (sum of all values in 5 day interval) and then use
t.rast.aggregate to compute the yearly maximums on the output of the first
operation.

Best regards
Soeren

[1] https://grass.osgeo.org/grass72/manuals/t.rast.accumulate.html

Am 06.04.2017 06:31 schrieb "RichardCooper" <richtcooper@hotmail.com>:

Hi,

I have a data set containing multiple years of daily raster layers and
would
like to aggregate and output annual raster layers of 5-day extremes
(maxima).

Essentially, for each grid cell, I need to calculate rolling 5-day sums of
each year and then find the annual max of the latter sums, and output as a
series of raster layers of annual 5-day extremes.

However, I'm trying to work out the best way in GRASS of doing this.
Overall
t.rast.aggregate comes closest to the type of functionality needed. I've
also looked at t.rast3d.mapcalc, but unsure if calculating a rolling sum is
possible.

I'd be grateful for any suggestions on how I might approach this.

Thanks,

Richard

--
View this message in context: http://osgeo-org.1560.x6.nabbl
e.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014.html
Sent from the Grass - Users mailing list archive at Nabble.com.
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

4

I have a time series of rainfall data, and for each year I want to calculate
the five-day period with maximum rainfall. So I would need to calculate the
sum of day1 to day5, then day2 to day6, then day3 to day7 etc for the whole
year, and then output the maximum grid cell 5-day values for each year into
a single raster.

To do this in t.rast.accumulate, I can see how to set a temporal cycle of 1
year (cycle= "12 months"), but not sure how to specify such a rolling sum
calculation of 5 days as described above. The default method is the 'mean'
as indicated in r.series.accumulation? I'm not too sure how the accumulation
is applied in the module.

Best regards,
Richard

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316076.html
Sent from the Grass - Users mailing list archive at Nabble.com.

5

Hi Richard

t.rast.aggregate might be the function you are looking for. It has the method 'max'.

For the cycling you might need to somehow loop over the data with different starting days.

Best, Mira

On 06/04/17 11:09, RichardCooper wrote:

I have a time series of rainfall data, and for each year I want to calculate
the five-day period with maximum rainfall. So I would need to calculate the
sum of day1 to day5, then day2 to day6, then day3 to day7 etc for the whole
year, and then output the maximum grid cell 5-day values for each year into
a single raster.

To do this in t.rast.accumulate, I can see how to set a temporal cycle of 1
year (cycle= "12 months"), but not sure how to specify such a rolling sum
calculation of 5 days as described above. The default method is the 'mean'
as indicated in r.series.accumulation? I'm not too sure how the accumulation
is applied in the module.

Best regards,
Richard

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316076.html
Sent from the Grass - Users mailing list archive at Nabble.com.
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

--
Dr. Mira Kattwinkel
Quantitative Landscape Ecology
Institute for Environmental Sciences
University of Koblenz-Landau
Fortstraße 7
76829 Landau
Germany
Phone: + 49 6341 280-31553
Office: Building I, Room 2.02

6

Hi,

2017-04-06 11:09 GMT+02:00 RichardCooper <richtcooper@hotmail.com>:

I have a time series of rainfall data, and for each year I want to calculate
the five-day period with maximum rainfall. So I would need to calculate the
sum of day1 to day5, then day2 to day6, then day3 to day7 etc for the whole
year, and then output the maximum grid cell 5-day values for each year into
a single raster.

What you need is a temporal moving window with the size of 5 days to
compute for each day the 5 day aggregate of the future.
You can convert your time series into a voxel dataset (3d raster) and
use r3.mapcalc with the neighbor index operator map[y][z] (if i
remember correctly):

agg_map3d = map3d[0][0][0] + map3d[0][0][1] + ... map3d[0][0][4]

Or you use t.rast.algebra [1] and the temporal neighbor operator strds[t]:

agg_strds = prec_strds[0] + prec_strds[1] + ... prec_strds[4]

Best regards
Soeren

[1] https://grass.osgeo.org/grass73/manuals/t.rast.algebra.html

To do this in t.rast.accumulate, I can see how to set a temporal cycle of 1
year (cycle= "12 months"), but not sure how to specify such a rolling sum
calculation of 5 days as described above. The default method is the 'mean'
as indicated in r.series.accumulation? I'm not too sure how the accumulation
is applied in the module.

Best regards,
Richard

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316076.html
Sent from the Grass - Users mailing list archive at Nabble.com.
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

7

Many thanks for the suggestions.

I'd really like to use the voxel approach but am getting the following error
on running t.rast.to.rast3 to create a voxel.

I've increased the hard and soft limits for open files to 40000 on my
system, but I still am unable to convert a space time raster dataset into a
3D raster map with only 5000 raster layers (note: I have 50K daily raster
layers that I need to analyse i voxel)

$ cat /proc/sys/fs/file-max
800532
$ ulimit -Sn
400000
$ ulimit -Hn
400000

Thanks,
Richard

Creation of STRDS and error on trying to convert to 3D data set:
t.create output=capha_test_5 temporaltype=relative semantictype=max
title=cahpa_test_5 description=cahpa_test_5

t.register --overwrite --verbose input=capha_test_5@cahpa
file=/home/rcooper/grassdatacl/climdata/cahpa/.tmp/rcooper-dell/3528.4
start=1 unit=day increment=1
Gathering map information...
Registering maps in the temporal database...
Registering maps in the space time dataset...
Updating space time dataset...
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_5@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_50@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_25@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_10@cahpa>

t.rast.to.rast3 --overwrite --verbose input=capha_test_5@cahpa
output=cahpa_test_5_3d
Traceback (most recent call last):
  File "/usr/local/grass-7.2.0/scripts/t.rast.to.rast3",
line 194, in <module>
    main()
  File "/usr/local/grass-7.2.0/scripts/t.rast.to.rast3",
line 152, in main
    output=output, overwrite=grass.overwrite())
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 408, in run_command
    ps = start_command(*args, **kwargs)
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 377, in start_command
    return Popen(args, **popts)
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 74, in __init__
    subprocess.Popen.__init__(self, args, **kwargs)
  File "/usr/lib/python2.7/subprocess.py", line 710, in
__init__
    errread, errwrite)
  File "/usr/lib/python2.7/subprocess.py", line 1327, in
_execute_child
    raise child_exception
OSError: [Errno 7] Argument list too long
(Fri Apr 7 17:30:33 2017) Command finished (2 sec)

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316291.html
Sent from the Grass - Users mailing list archive at Nabble.com.

8

I should add that t.rast.to.rast3 works for a STRDS containing 1000 layers,
but if more than 5000 layers I get the abovementioned error after increasing
the open file limits (to 400000) on my notebook.

The number of open files on the system:
lsof | wc -l
156036

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316295.html
Sent from the Grass - Users mailing list archive at Nabble.com.

9

On 07/04/17 12:45, RichardCooper wrote:

Many thanks for the suggestions.

I'd really like to use the voxel approach but am getting the following error
on running t.rast.to.rast3 to create a voxel.

I've increased the hard and soft limits for open files to 40000 on my
system, but I still am unable to convert a space time raster dataset into a
3D raster map with only 5000 raster layers (note: I have 50K daily raster
layers that I need to analyse i voxel)

$ cat /proc/sys/fs/file-max
800532
$ ulimit -Sn
400000
$ ulimit -Hn
400000

AFAICT, the issue is not with the number of open files, but with a list of map names (arguments) that is too long for the run_command call that calls r.to.rast3. A solution to this could be to allow as input to r.to.rast3 a file with map names (such as for r.series).

You should create a bug report for this.

Moritz

Thanks,
Richard

Creation of STRDS and error on trying to convert to 3D data set:
t.create output=capha_test_5 temporaltype=relative semantictype=max
title=cahpa_test_5 description=cahpa_test_5

t.register --overwrite --verbose input=capha_test_5@cahpa
file=/home/rcooper/grassdatacl/climdata/cahpa/.tmp/rcooper-dell/3528.4
start=1 unit=day increment=1
Gathering map information...
Registering maps in the temporal database...
Registering maps in the space time dataset...
Updating space time dataset...
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_5@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_50@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_25@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_10@cahpa>

t.rast.to.rast3 --overwrite --verbose input=capha_test_5@cahpa
output=cahpa_test_5_3d
Traceback (most recent call last):
  File "/usr/local/grass-7.2.0/scripts/t.rast.to.rast3",
line 194, in <module>
    main()
  File "/usr/local/grass-7.2.0/scripts/t.rast.to.rast3",
line 152, in main
    output=output, overwrite=grass.overwrite())
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 408, in run_command
    ps = start_command(*args, **kwargs)
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 377, in start_command
    return Popen(args, **popts)
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 74, in __init__
    subprocess.Popen.__init__(self, args, **kwargs)
  File "/usr/lib/python2.7/subprocess.py", line 710, in
__init__
    errread, errwrite)
  File "/usr/lib/python2.7/subprocess.py", line 1327, in
_execute_child
    raise child_exception
OSError: [Errno 7] Argument list too long
(Fri Apr 7 17:30:33 2017) Command finished (2 sec)

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316291.html
Sent from the Grass - Users mailing list archive at Nabble.com.
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

10

Patching t.rast.to.rast3 to support file input requires the
modification of r.to.rast3 to implement file support. This may not
happen in the near future.
I would suggest to use t.rast.algebra or to compute the 5day rainfall
extrema in 5 year chunks using t.rast.to.rast3 + r3.mapcalc. Use
t.rast.extract and its "where" option to compute the 5year strds
chunks.

Best regards
Soeren

2017-04-07 13:37 GMT+02:00 Moritz Lennert <mlennert@club.worldonline.be>:

On 07/04/17 12:45, RichardCooper wrote:

Many thanks for the suggestions.

I'd really like to use the voxel approach but am getting the following
error
on running t.rast.to.rast3 to create a voxel.

I've increased the hard and soft limits for open files to 40000 on my
system, but I still am unable to convert a space time raster dataset into
a
3D raster map with only 5000 raster layers (note: I have 50K daily raster
layers that I need to analyse i voxel)

$ cat /proc/sys/fs/file-max
800532
$ ulimit -Sn
400000
$ ulimit -Hn
400000

AFAICT, the issue is not with the number of open files, but with a list of
map names (arguments) that is too long for the run_command call that calls
r.to.rast3. A solution to this could be to allow as input to r.to.rast3 a
file with map names (such as for r.series).

You should create a bug report for this.

Moritz

Thanks,
Richard

Creation of STRDS and error on trying to convert to 3D data set:
t.create output=capha_test_5 temporaltype=relative semantictype=max
title=cahpa_test_5 description=cahpa_test_5

t.register --overwrite --verbose input=capha_test_5@cahpa
file=/home/rcooper/grassdatacl/climdata/cahpa/.tmp/rcooper-dell/3528.4
start=1 unit=day increment=1
Gathering map information...
Registering maps in the temporal database...
Registering maps in the space time dataset...
Updating space time dataset...
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_5@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_50@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_25@cahpa>
Update metadata, spatial and temporal extent from all registered maps of
<capha_test_10@cahpa>

t.rast.to.rast3 --overwrite --verbose input=capha_test_5@cahpa
output=cahpa_test_5_3d
Traceback (most recent call last):
  File "/usr/local/grass-7.2.0/scripts/t.rast.to.rast3",
line 194, in <module>
    main()
  File "/usr/local/grass-7.2.0/scripts/t.rast.to.rast3",
line 152, in main
    output=output, overwrite=grass.overwrite())
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 408, in run_command
    ps = start_command(*args, **kwargs)
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 377, in start_command
    return Popen(args, **popts)
  File
"/usr/local/grass-7.2.0/etc/python/grass/script/core.py",
line 74, in __init__
    subprocess.Popen.__init__(self, args, **kwargs)
  File "/usr/lib/python2.7/subprocess.py", line 710, in
__init__
    errread, errwrite)
  File "/usr/lib/python2.7/subprocess.py", line 1327, in
_execute_child
    raise child_exception
OSError: [Errno 7] Argument list too long
(Fri Apr 7 17:30:33 2017) Command finished (2 sec)

--
View this message in context:
http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316291.html
Sent from the Grass - Users mailing list archive at Nabble.com.
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-user

11

I'm looking at t.rast.algebra, but I'm getting the following error which I
haven't been able to resolve. I'm not sure if I need to create an absolute
SRTDS as shown in the first error below. I then try and register a set of
raster layers with an absolute SRTDS but get the second error. I'd be
grateful for any suggestions.

Error 1:
On using t.rast.algebra with a relative SRTDS I get the following error:

Error 2:
Attempting to register raster layers in an absolute SRTDS:

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316354.html
Sent from the Grass - Users mailing list archive at Nabble.com.

12

With a manual work around I was able to use t.rast.algebra, but needed to
both manually create a relative STRDS and then register the output rasters
from t.rast.algebra into the latter.

The error encountered in using t.rast.algebra (I'm not sure if me or a bug)
appears to be that an absolute STRDS is being created by t.rast.algebra but
the generated raster layers have a relative timestamp (e.g. 'Timestamp: 1
day / 2 days').

This expression output the required summed layers (and I then manually
created the relative STRDS and registered the output layers with
t.register):
t.rast.algebra -s --verbose expression=ag_rel200=(cahpa_rel200[0] +
cahpa_rel200[1] + cahpa_rel200[2] + cahpa_rel200[3] + cahpa_rel200[4])
basename=summed

Regarding the second error with t.register as noted above: the input raster
layers have a relative timestamp and I understand that if I wish to create
an absolute STRDS, then I need to modify each layer's timestamp from
relative to absolute accordingly.

Please let me know if the above makes sense.

Best regards,
Richard

--
View this message in context: http://osgeo-org.1560.x6.nabble.com/Temporal-framework-calculating-annual-5-day-extremes-tp5316014p5316391.html
Sent from the Grass - Users mailing list archive at Nabble.com.