Temporal framework: calculating annual 5-day extremes

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

Temporal framework: calculating annual 5-day extremes

RichardCooper
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
SBP
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

SBP
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 <[hidden email]> 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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user


_______________________________________________
grass-user mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

Sören Gebbert
In reply to this post by RichardCooper
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


Am 06.04.2017 06:31 schrieb "RichardCooper" <[hidden email]>:
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
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user

_______________________________________________
grass-user mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

RichardCooper
This post was updated on .
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'. I'm not too sure how the accumulation is applied in the module.

Best regards,
Richard
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

Mira Kattwinkel
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
> [hidden email]
> 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

_______________________________________________
grass-user mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

Sören Gebbert
In reply to this post by RichardCooper
Hi,

2017-04-06 11:09 GMT+02:00 RichardCooper <[hidden email]>:
> 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[x][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
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/grass-user
_______________________________________________
grass-user mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

RichardCooper
This post was updated on .
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)                            


System Info                                                                    
GRASS version: 7.2.0                                                            
GRASS SVN revision: exported                                                    
Build date: 2017-04-07                                                          
Build platform: x86_64-pc-linux-gnu                                            
GDAL: 1.11.2                                                                    
PROJ.4: 4.9.0                                                                  
GEOS: 3.4.2                                                                    
SQLite: 3.8.2                                                                  
Python: 2.7.6                                                                  
wxPython: 2.8.12.1                                                              
Platform: Linux-3.13.0-95-generic-x86_64-with-LinuxMint-17.1-rebecca            
                       








Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

RichardCooper
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
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

Moritz Lennert
In reply to this post by RichardCooper
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
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/grass-user
>


_______________________________________________
grass-user mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

Sören Gebbert
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 <[hidden email]>:

> 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
>> [hidden email]
>> https://lists.osgeo.org/mailman/listinfo/grass-user
>>
>
>
> _______________________________________________
> grass-user mailing list
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/grass-user
_______________________________________________
grass-user mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/grass-user
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

RichardCooper
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:

t.rast.algebra -g expression=agg4_cahpa = capha_rel_2[0] + capha_rel_2[1] + capha_rel_2[2] + capha_rel_2[3] + capha_rel_2[4] basename=agg4

NOTE: I added the -g flag thinking that this may be a solution, but I still get an error.

ERROR: Temporal type of space time dataset <agg4_cahpa@cahpa> and map <agg4_0@cahpa> are different
Traceback (most recent call last):
  File "/usr/local/grass-7.2.0/scripts/t.rast.algebra", line
113, in <module>
    sys.exit(main())
  File "/usr/local/grass-7.2.0/scripts/t.rast.algebra", line
105, in main
    pc = p.parse(expression, basename,
grass.script.overwrite())
  File "/usr/local/grass-7.2.0/etc/python/grass/temporal/tem
poral_raster_algebra.py", line 96, in parse
    self.parser.parse(expression)
  File "/usr/lib/python2.7/dist-packages/ply/yacc.py", line
269, in parse
    return
self.parseopt_notrack(input,lexer,debug,tracking,tokenfunc)
  File "/usr/lib/python2.7/dist-packages/ply/yacc.py", line
975, in parseopt_notrack
    p.callable(pslice)
  File "/usr/local/grass-7.2.0/etc/python/grass/temporal/tem
poral_raster_algebra.py", line 107, in p_statement_assign
    TemporalRasterBaseAlgebraParser.p_statement_assign(self,
t)
  File "/usr/local/grass-7.2.0/etc/python/grass/temporal/tem
poral_raster_base_algebra.py", line 666, in
p_statement_assign
    success = resultstds.register_map(map_i, dbif)
  File "/usr/local/grass-7.2.0/etc/python/grass/temporal/abs
tract_space_time_dataset.py", line 2110, in register_map
    'map': map.get_map_id()})
  File "/usr/local/grass-7.2.0/etc/python/grass/pygrass/mess
ages/__init__.py", line 269, in fatal
    raise FatalError(message)
grass.exceptions.FatalError: Temporal type of space time
dataset <agg4_cahpa@cahpa> and map <agg4_0@cahpa> are
different

Error 2:
Attempting to register raster layers in an absolute SRTDS:
t.register -i --overwrite input=cahpa_abs_2@cahpa file=/home/rcooper/grassdatacl/climdata/cahpa/.tmp/rcooper-dell/32525.2 start=1951-01-01 increment=1 days type=raster
Gathering map information...
ERROR: invalid literal for int() with base 10: '1951-01-01'
Reply | Threaded
Open this post in threaded view
|

Re: Temporal framework: calculating annual 5-day extremes

RichardCooper
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