[GRASS-dev] Sampling STRDS at point location(s)

Dear all,

On a regularly basis I am asked by colleagues to extract raster based time series data (with > 20k maps) for a (set of) point(s).

After both t.vect.observe.strds, t.sample and t.vect.what.strds seem to have a different purposes, I do this using a little shell script which mainly uses r.what (please find it attached).

My first question is, did I overlook an existing module which does the job?

The data behind g.gui.tplot would actually be what I was interested in (also for more than just one point)…

I could further-develop the attached shell-script into a new Python-AddOn if that would be of interest, but some improvements in r.what would not make it necessary to have an new module. The main problem is the 400 maps limit in r.what (which is named as a “to do” in the r.what manual). Are there any plans on working on this limitation?

Another nice feature in r.what would be to be able to name an attribute column which is written to the site_name column of the r.what output when coordinates are taken from a vector map (like it is done when coordinates with a site name are used). The latter can be achieved for the “cat” column when “v.to.db”-output is piped to r.what. However it would be nice to be able to use clear text names here as well (which at the moment – as far as I can see – requires to write the coordinates to the attribute table (v.to.db) and then fetch them along with the full-text site name (using v.db.select) in order to pipe them to r.what).

As you can see in the attached script I also transpose the output of r.what because for time series with thousands of maps my colleagues will likely get problems with the number of columns in the software they are using. Therfore a switch in r.what between “horizontal” and “vertical” output (like in v.db.select) would be useful when dealing with thousands of maps in r.what…

Any thought? Would you prefer an enhancement request for r.what or a new addon or none of both?

Cheers

Stefan

That is because I forgot to attach it... :wink:
Here it comes...

-----Original Message-----
From: Sören Gebbert [mailto:soerengebbert@googlemail.com]
Sent: 16. januar 2015 10:38
To: Blumentrath, Stefan
Subject: Re: [GRASS-dev] Sampling STRDS at point location(s)

Hi Stefan,
can you please attach your referenced script? It did not reached me.

Best regards
Sören

(attachments)

tmp.sh (2.26 KB)

Hi Stefan,
from my understanding of your requirements and the study of your
script i would say, that t.vect.observe.strds should be the right tool
for you.

This module has two inputs. The first input is a vector map layer with
vector points. The second input is one or several space time raster
datasets that should be sampled over time at the vector point
positions. The space time raster dataset will be sampled over its
whole temporal extent (from start to end). You can adjust the range
using the temporal where condition t_where. You need to specify a
column name for each input space time raster dataset.

The result is a new space time vector dataset that contains a single
new vector map and as many time stamped attribute tables linked to the
new vector map as raster map layer are present in the input space time
raster dataset.
Hence, for each time step in the space time raster dataset a new
attribute table is created. The GRASS GIS Temporal Framework allows to
time stamp attribute tables that can be linked to a single vector map
layer.
The module v.what.rast is used for sampling the time stamped raster
map layers. All sampled values of a single time stamped raster map
layer are written into a new time stamped attribute table.

You can use t.vect.db.select to print attribute values of the space
time vector dataset to stdout.

Example:
Vector map layer sample_points contains the vector points at which the
STRDS temperature_daily should be sampled over time. Output is a new
STVDS sampled_stvds and a new vector map sample_vect that links to the
time stamped attribute tables.

t.vect.observe.strds input=sample_points strds=temperature_daily
output=sample_stvds vect=sample_vect columns=temp

The output is generated with t.vect.db.select and should look like this:

t.vect.db.select input=sample_stvds columns=cat,temp

start_time | end_time | cat | temp
2001-01-01 | 2001-01-02 | 1| -6
2001-01-01 | 2001-01-02 | 2| -4
2001-01-01 | 2001-01-02 | 3| -8
....
2001-01-31 | 2001-02-01 | 1| -3
2001-01-31 | 2001-02-01 | 2| -2
2001-01-31 | 2001-02-01 | 3| -5

I hope this is what you are searching for
Best regards
Soeren

2015-01-16 10:08 GMT+01:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

Dear all,

On a regularly basis I am asked by colleagues to extract raster based time
series data (with > 20k maps) for a (set of) point(s).

After both t.vect.observe.strds, t.sample and t.vect.what.strds seem to have
a different purposes, I do this using a little shell script which mainly
uses r.what (please find it attached).

My first question is, did I overlook an existing module which does the job?

The data behind g.gui.tplot would actually be what I was interested in (also
for more than just one point)...

I could further-develop the attached shell-script into a new Python-AddOn if
that would be of interest, but some improvements in r.what would not make it
necessary to have an new module. The main problem is the 400 maps limit in
r.what (which is named as a “to do” in the r.what manual). Are there any
plans on working on this limitation?

Another nice feature in r.what would be to be able to name an attribute
column which is written to the site_name column of the r.what output when
coordinates are taken from a vector map (like it is done when coordinates
with a site name are used). The latter can be achieved for the “cat” column
when “v.to.db”-output is piped to r.what. However it would be nice to be
able to use clear text names here as well (which at the moment – as far as I
can see – requires to write the coordinates to the attribute table (v.to.db)
and then fetch them along with the full-text site name (using v.db.select)
in order to pipe them to r.what).

As you can see in the attached script I also transpose the output of r.what
because for time series with thousands of maps my colleagues will likely get
problems with the number of columns in the software they are using. Therfore
a switch in r.what between “horizontal” and “vertical” output (like in
v.db.select) would be useful when dealing with thousands of maps in r.what…

Any thought? Would you prefer an enhancement request for r.what or a new
addon or none of both?

Cheers

Stefan

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

Hi Sören,

And thanks for clarification.
I thought that t.vect.observe.strds would do what I am up to, but it is rather time-consuming for >20k maps.
So I thought it`s main purpose must be to generate a new STVDS from an STRDS.
Using r.what is much, much faster (15 min compared to > 12h for t.vect.observe.strds (and here I have not run t.vect.db.select yet because t.vect.observe.strds is still running (since yesterday evening)). When I only want a matrix (text output) as a result, creating a new STVDS seems to be too much overhead...
This was the main reason for me to use r.what...

Cheers
Stefan

-----Original Message-----
From: Sören Gebbert [mailto:soerengebbert@googlemail.com]
Sent: 16. januar 2015 13:33
To: Blumentrath, Stefan
Cc: GRASS developers list (grass-dev@lists.osgeo.org)
Subject: Re: [GRASS-dev] Sampling STRDS at point location(s)

Hi Stefan,
from my understanding of your requirements and the study of your script i would say, that t.vect.observe.strds should be the right tool for you.

This module has two inputs. The first input is a vector map layer with vector points. The second input is one or several space time raster datasets that should be sampled over time at the vector point positions. The space time raster dataset will be sampled over its whole temporal extent (from start to end). You can adjust the range using the temporal where condition t_where. You need to specify a column name for each input space time raster dataset.

The result is a new space time vector dataset that contains a single new vector map and as many time stamped attribute tables linked to the new vector map as raster map layer are present in the input space time raster dataset.
Hence, for each time step in the space time raster dataset a new attribute table is created. The GRASS GIS Temporal Framework allows to time stamp attribute tables that can be linked to a single vector map layer.
The module v.what.rast is used for sampling the time stamped raster map layers. All sampled values of a single time stamped raster map layer are written into a new time stamped attribute table.

You can use t.vect.db.select to print attribute values of the space time vector dataset to stdout.

Example:
Vector map layer sample_points contains the vector points at which the STRDS temperature_daily should be sampled over time. Output is a new STVDS sampled_stvds and a new vector map sample_vect that links to the time stamped attribute tables.

t.vect.observe.strds input=sample_points strds=temperature_daily output=sample_stvds vect=sample_vect columns=temp

The output is generated with t.vect.db.select and should look like this:

t.vect.db.select input=sample_stvds columns=cat,temp

start_time | end_time | cat | temp
2001-01-01 | 2001-01-02 | 1| -6
2001-01-01 | 2001-01-02 | 2| -4
2001-01-01 | 2001-01-02 | 3| -8
....
2001-01-31 | 2001-02-01 | 1| -3
2001-01-31 | 2001-02-01 | 2| -2
2001-01-31 | 2001-02-01 | 3| -5

I hope this is what you are searching for Best regards Soeren

2015-01-16 10:08 GMT+01:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

Dear all,

On a regularly basis I am asked by colleagues to extract raster based
time series data (with > 20k maps) for a (set of) point(s).

After both t.vect.observe.strds, t.sample and t.vect.what.strds seem
to have a different purposes, I do this using a little shell script
which mainly uses r.what (please find it attached).

My first question is, did I overlook an existing module which does the job?

The data behind g.gui.tplot would actually be what I was interested in
(also for more than just one point)...

I could further-develop the attached shell-script into a new
Python-AddOn if that would be of interest, but some improvements in
r.what would not make it necessary to have an new module. The main
problem is the 400 maps limit in r.what (which is named as a “to do”
in the r.what manual). Are there any plans on working on this limitation?

Another nice feature in r.what would be to be able to name an
attribute column which is written to the site_name column of the
r.what output when coordinates are taken from a vector map (like it is
done when coordinates with a site name are used). The latter can be
achieved for the “cat” column when “v.to.db”-output is piped to
r.what. However it would be nice to be able to use clear text names
here as well (which at the moment – as far as I can see – requires to
write the coordinates to the attribute table (v.to.db) and then fetch
them along with the full-text site name (using v.db.select) in order to pipe them to r.what).

As you can see in the attached script I also transpose the output of
r.what because for time series with thousands of maps my colleagues
will likely get problems with the number of columns in the software
they are using. Therfore a switch in r.what between “horizontal” and
“vertical” output (like in
v.db.select) would be useful when dealing with thousands of maps in
r.what…

Any thought? Would you prefer an enhancement request for r.what or a
new addon or none of both?

Cheers

Stefan

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

Hi Stefan,
the most time consuming part of t.vect.observe.strds is the SQL
database handling in v.what.rast. Indeed, r.what is in this regards
much faster.

You can extent t.vect.observe.strds to use r.what instead of
v.what.rast. In this case, an output file with values is generated
rather than a new vector map and a STVDS. This module can be called
t.rast.what then? :slight_smile:

Best regards
Sören

2015-01-16 13:49 GMT+01:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

Hi Sören,

And thanks for clarification.
I thought that t.vect.observe.strds would do what I am up to, but it is rather time-consuming for >20k maps.
So I thought it`s main purpose must be to generate a new STVDS from an STRDS.
Using r.what is much, much faster (15 min compared to > 12h for t.vect.observe.strds (and here I have not run t.vect.db.select yet because t.vect.observe.strds is still running (since yesterday evening)). When I only want a matrix (text output) as a result, creating a new STVDS seems to be too much overhead...
This was the main reason for me to use r.what...

Cheers
Stefan

-----Original Message-----
From: Sören Gebbert [mailto:soerengebbert@googlemail.com]
Sent: 16. januar 2015 13:33
To: Blumentrath, Stefan
Cc: GRASS developers list (grass-dev@lists.osgeo.org)
Subject: Re: [GRASS-dev] Sampling STRDS at point location(s)

Hi Stefan,
from my understanding of your requirements and the study of your script i would say, that t.vect.observe.strds should be the right tool for you.

This module has two inputs. The first input is a vector map layer with vector points. The second input is one or several space time raster datasets that should be sampled over time at the vector point positions. The space time raster dataset will be sampled over its whole temporal extent (from start to end). You can adjust the range using the temporal where condition t_where. You need to specify a column name for each input space time raster dataset.

The result is a new space time vector dataset that contains a single new vector map and as many time stamped attribute tables linked to the new vector map as raster map layer are present in the input space time raster dataset.
Hence, for each time step in the space time raster dataset a new attribute table is created. The GRASS GIS Temporal Framework allows to time stamp attribute tables that can be linked to a single vector map layer.
The module v.what.rast is used for sampling the time stamped raster map layers. All sampled values of a single time stamped raster map layer are written into a new time stamped attribute table.

You can use t.vect.db.select to print attribute values of the space time vector dataset to stdout.

Example:
Vector map layer sample_points contains the vector points at which the STRDS temperature_daily should be sampled over time. Output is a new STVDS sampled_stvds and a new vector map sample_vect that links to the time stamped attribute tables.

t.vect.observe.strds input=sample_points strds=temperature_daily output=sample_stvds vect=sample_vect columns=temp

The output is generated with t.vect.db.select and should look like this:

t.vect.db.select input=sample_stvds columns=cat,temp

start_time | end_time | cat | temp
2001-01-01 | 2001-01-02 | 1| -6
2001-01-01 | 2001-01-02 | 2| -4
2001-01-01 | 2001-01-02 | 3| -8
....
2001-01-31 | 2001-02-01 | 1| -3
2001-01-31 | 2001-02-01 | 2| -2
2001-01-31 | 2001-02-01 | 3| -5

I hope this is what you are searching for Best regards Soeren

2015-01-16 10:08 GMT+01:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

Dear all,

On a regularly basis I am asked by colleagues to extract raster based
time series data (with > 20k maps) for a (set of) point(s).

After both t.vect.observe.strds, t.sample and t.vect.what.strds seem
to have a different purposes, I do this using a little shell script
which mainly uses r.what (please find it attached).

My first question is, did I overlook an existing module which does the job?

The data behind g.gui.tplot would actually be what I was interested in
(also for more than just one point)...

I could further-develop the attached shell-script into a new
Python-AddOn if that would be of interest, but some improvements in
r.what would not make it necessary to have an new module. The main
problem is the 400 maps limit in r.what (which is named as a “to do”
in the r.what manual). Are there any plans on working on this limitation?

Another nice feature in r.what would be to be able to name an
attribute column which is written to the site_name column of the
r.what output when coordinates are taken from a vector map (like it is
done when coordinates with a site name are used). The latter can be
achieved for the “cat” column when “v.to.db”-output is piped to
r.what. However it would be nice to be able to use clear text names
here as well (which at the moment – as far as I can see – requires to
write the coordinates to the attribute table (v.to.db) and then fetch
them along with the full-text site name (using v.db.select) in order to pipe them to r.what).

As you can see in the attached script I also transpose the output of
r.what because for time series with thousands of maps my colleagues
will likely get problems with the number of columns in the software
they are using. Therfore a switch in r.what between “horizontal” and
“vertical” output (like in
v.db.select) would be useful when dealing with thousands of maps in
r.what…

Any thought? Would you prefer an enhancement request for r.what or a
new addon or none of both?

Cheers

Stefan

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

Hi Sören,

Your idea to include this functionality into t.vect.observe.strds is probably the most elegant solution, as long as no C-developer prefers to remove the 400 maps limit in r.what...
If the latter is not the case (or if one would run into memory problems when running r.what with 20k maps anyway) I shall try to come up with a proposal next week.
Have a nice weekend.

Cheers
Stefan

-----Original Message-----
From: Sören Gebbert [mailto:soerengebbert@googlemail.com]
Sent: 16. januar 2015 14:10
To: Blumentrath, Stefan
Cc: GRASS developers list (grass-dev@lists.osgeo.org)
Subject: Re: [GRASS-dev] Sampling STRDS at point location(s)

Hi Stefan,
the most time consuming part of t.vect.observe.strds is the SQL database handling in v.what.rast. Indeed, r.what is in this regards much faster.

You can extent t.vect.observe.strds to use r.what instead of v.what.rast. In this case, an output file with values is generated rather than a new vector map and a STVDS. This module can be called t.rast.what then? :slight_smile:

Best regards
Sören

2015-01-16 13:49 GMT+01:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

Hi Sören,

And thanks for clarification.
I thought that t.vect.observe.strds would do what I am up to, but it is rather time-consuming for >20k maps.
So I thought it`s main purpose must be to generate a new STVDS from an STRDS.
Using r.what is much, much faster (15 min compared to > 12h for t.vect.observe.strds (and here I have not run t.vect.db.select yet because t.vect.observe.strds is still running (since yesterday evening)). When I only want a matrix (text output) as a result, creating a new STVDS seems to be too much overhead...
This was the main reason for me to use r.what...

Cheers
Stefan

-----Original Message-----
From: Sören Gebbert [mailto:soerengebbert@googlemail.com]
Sent: 16. januar 2015 13:33
To: Blumentrath, Stefan
Cc: GRASS developers list (grass-dev@lists.osgeo.org)
Subject: Re: [GRASS-dev] Sampling STRDS at point location(s)

Hi Stefan,
from my understanding of your requirements and the study of your script i would say, that t.vect.observe.strds should be the right tool for you.

This module has two inputs. The first input is a vector map layer with vector points. The second input is one or several space time raster datasets that should be sampled over time at the vector point positions. The space time raster dataset will be sampled over its whole temporal extent (from start to end). You can adjust the range using the temporal where condition t_where. You need to specify a column name for each input space time raster dataset.

The result is a new space time vector dataset that contains a single new vector map and as many time stamped attribute tables linked to the new vector map as raster map layer are present in the input space time raster dataset.
Hence, for each time step in the space time raster dataset a new attribute table is created. The GRASS GIS Temporal Framework allows to time stamp attribute tables that can be linked to a single vector map layer.
The module v.what.rast is used for sampling the time stamped raster map layers. All sampled values of a single time stamped raster map layer are written into a new time stamped attribute table.

You can use t.vect.db.select to print attribute values of the space time vector dataset to stdout.

Example:
Vector map layer sample_points contains the vector points at which the STRDS temperature_daily should be sampled over time. Output is a new STVDS sampled_stvds and a new vector map sample_vect that links to the time stamped attribute tables.

t.vect.observe.strds input=sample_points strds=temperature_daily
output=sample_stvds vect=sample_vect columns=temp

The output is generated with t.vect.db.select and should look like this:

t.vect.db.select input=sample_stvds columns=cat,temp

start_time | end_time | cat | temp
2001-01-01 | 2001-01-02 | 1| -6
2001-01-01 | 2001-01-02 | 2| -4
2001-01-01 | 2001-01-02 | 3| -8
....
2001-01-31 | 2001-02-01 | 1| -3
2001-01-31 | 2001-02-01 | 2| -2
2001-01-31 | 2001-02-01 | 3| -5

I hope this is what you are searching for Best regards Soeren

2015-01-16 10:08 GMT+01:00 Blumentrath, Stefan <Stefan.Blumentrath@nina.no>:

Dear all,

On a regularly basis I am asked by colleagues to extract raster based
time series data (with > 20k maps) for a (set of) point(s).

After both t.vect.observe.strds, t.sample and t.vect.what.strds seem
to have a different purposes, I do this using a little shell script
which mainly uses r.what (please find it attached).

My first question is, did I overlook an existing module which does the job?

The data behind g.gui.tplot would actually be what I was interested
in (also for more than just one point)...

I could further-develop the attached shell-script into a new
Python-AddOn if that would be of interest, but some improvements in
r.what would not make it necessary to have an new module. The main
problem is the 400 maps limit in r.what (which is named as a “to do”
in the r.what manual). Are there any plans on working on this limitation?

Another nice feature in r.what would be to be able to name an
attribute column which is written to the site_name column of the
r.what output when coordinates are taken from a vector map (like it
is done when coordinates with a site name are used). The latter can
be achieved for the “cat” column when “v.to.db”-output is piped to
r.what. However it would be nice to be able to use clear text names
here as well (which at the moment – as far as I can see – requires to
write the coordinates to the attribute table (v.to.db) and then fetch
them along with the full-text site name (using v.db.select) in order to pipe them to r.what).

As you can see in the attached script I also transpose the output of
r.what because for time series with thousands of maps my colleagues
will likely get problems with the number of columns in the software
they are using. Therfore a switch in r.what between “horizontal” and
“vertical” output (like in
v.db.select) would be useful when dealing with thousands of maps in
r.what…

Any thought? Would you prefer an enhancement request for r.what or a
new addon or none of both?

Cheers

Stefan

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

On Fri, Jan 16, 2015 at 4:05 PM, Blumentrath, Stefan
<Stefan.Blumentrath@nina.no> wrote:

Hi Sören,

Your idea to include this functionality into t.vect.observe.strds is probably the most elegant solution, as long as no C-developer prefers to remove the 400 maps limit in r.what...

In r.patch a similar limitation problem was solved in
http://trac.osgeo.org/grass/changeset/23499

Would be great to get rid of that also in r.what.

Markus

BTW: is not this ticket for r.what outdated and can be closed?
https://trac.osgeo.org/grass/ticket/32

-----Original Message-----
From: neteler.osgeo@gmail.com [mailto:neteler.osgeo@gmail.com] On Behalf Of Markus Neteler
Sent: 26. januar 2015 20:06
To: Blumentrath, Stefan
Cc: Sören Gebbert; GRASS developers list (grass-dev@lists.osgeo.org)
Subject: Re: [GRASS-dev] Sampling STRDS at point location(s)

On Fri, Jan 16, 2015 at 4:05 PM, Blumentrath, Stefan <Stefan.Blumentrath@nina.no> wrote:

Hi Sören,

Your idea to include this functionality into t.vect.observe.strds is probably the most elegant solution, as long as no C-developer prefers to remove the 400 maps limit in r.what...

In r.patch a similar limitation problem was solved in
http://trac.osgeo.org/grass/changeset/23499

Would be great to get rid of that also in r.what.

Markus