[GRASS-user] GRASS 6.4 How to extract point-values from multiple rasters

GRASS GRASS user
1

Dear GRASS community,

I have to extract isolated value from several raster.
I've found this script but I do not understand how to run it. Would someone have an other script with more comments (notably what are the input I have to give, what are the grass or python function/module) or explain me how to run this script.

I would also be interested by a bash script if someone has such a script.

Thank you for your support.

Kind regards

Lucien

2

Hi!
I recently had nearly the same problem … and solved it in different GIS Applications … (there were not so many points i needed) but during research I found this link

http://gis.stackexchange.com/questions/3538/extract-raster-values-at-points

Which you might find interesting especially when looking at
isectpntrst(in=“path/to/shapefile”, raster=“path/to/raster”, field=“fieldname”)

when it comes to bash …

or probably
http://gis-techniques.blogspot.co.at/2012/10/extract-raster-values-from-points.html

using R …

It always depends how you’d like to solve it … of course it can be done in GRASS too …

hope this helps

kind regards
Werner

···

On Thu, Apr 4, 2013 at 4:37 PM, BLANDENIER Lucien <lucien.blandenier@unine.ch> wrote:

Dear GRASS community,

I have to extract isolated value from several raster.
I’ve found this script but I do not understand how to run it. Would someone have an other script with more comments (notably what are the input I have to give, what are the grass or python function/module) or explain me how to run this script.

I would also be interested by a bash script if someone has such a script.

Thank you for your support.

Kind regards

Lucien

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

3

ups my fault - did not read the full link by myself … blush seems that ArcGIS is needed for the “lookalike” bash script :wink:

but anyway the other links should give you enough information too …

HTH
Werner

···

On Thu, Apr 4, 2013 at 4:44 PM, Werner Macho <werner.macho@gmail.com> wrote:

Hi!
I recently had nearly the same problem … and solved it in different GIS Applications … (there were not so many points i needed) but during research I found this link

http://gis.stackexchange.com/questions/3538/extract-raster-values-at-points

Which you might find interesting especially when looking at
isectpntrst(in=“path/to/shapefile”, raster=“path/to/raster”, field=“fieldname”)

when it comes to bash …

or probably
http://gis-techniques.blogspot.co.at/2012/10/extract-raster-values-from-points.html

using R …

It always depends how you’d like to solve it … of course it can be done in GRASS too …

hope this helps

kind regards

Werner

On Thu, Apr 4, 2013 at 4:37 PM, BLANDENIER Lucien <lucien.blandenier@unine.ch> wrote:

Dear GRASS community,

I have to extract isolated value from several raster.
I’ve found this script but I do not understand how to run it. Would someone have an other script with more comments (notably what are the input I have to give, what are the grass or python function/module) or explain me how to run this script.

I would also be interested by a bash script if someone has such a script.

Thank you for your support.

Kind regards

Lucien

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

4

On Thu, Apr 4, 2013 at 4:37 PM, BLANDENIER Lucien
<lucien.blandenier@unine.ch> wrote:

Dear GRASS community,

I have to extract isolated value from several raster.
I've found this script but I do not understand how to run it. Would someone have an other script with more comments (notably what are the input I have to give, what are the grass or python function/module) or explain me how to run this script.

The GRASS module r.what [0] extracts raster values at a given
location. As input you need to give the name of the raster map and the
coordinates as east,north

[0] http://grass.osgeo.org/grass64/manuals/r.what.html

Markus M

I would also be interested by a bash script if someone has such a script.

Thank you for your support.

Kind regards

Lucien

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

5

Alternatively you can do it using this script (python-gdal):

http://gis-lab.info/qa/extract-values-rasters-eng.html

Saber Razmjooei
Lutra Consulting
www.lutraconsulting.co.uk

-----Original Message-----
From: grass-user-bounces@lists.osgeo.org [mailto:grass-user-bounces@lists.osgeo.org] On Behalf Of Markus Metz
Sent: 04 April 2013 15:51
To: BLANDENIER Lucien
Cc: grass [grass-user@lists.osgeo.org]
Subject: Re: [GRASS-user] GRASS 6.4 How to extract point-values from multiple rasters

On Thu, Apr 4, 2013 at 4:37 PM, BLANDENIER Lucien <lucien.blandenier@unine.ch> wrote:

Dear GRASS community,

I have to extract isolated value from several raster.

I’ve found this script but I do not understand how to run it. Would someone have an other script with more comments (notably what are the input I have to give, what are the grass or python function/module) or explain me how to run this script.

The GRASS module r.what [0] extracts raster values at a given location. As input you need to give the name of the raster map and the coordinates as east,north

[0] http://grass.osgeo.org/grass64/manuals/r.what.html

Markus M

I would also be interested by a bash script if someone has such a script.

Thank you for your support.

Kind regards

Lucien

grass-user mailing list

grass-user@lists.osgeo.org

http://lists.osgeo.org/mailman/listinfo/grass-user

grass-user mailing list

grass-user@lists.osgeo.org

http://lists.osgeo.org/mailman/listinfo/grass-user

This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed. If you have received this email in error please notify the system manager. This message contains confidential information and is intended only for the individual named. If you are not the named addressee you should not disseminate, distribute or copy this e-mail. Please notify the sender immediately by e-mail if you have received this e-mail by mistake and delete this e-mail from your system. If you are not the intended recipient you are notified that disclosing, copying, distributing or taking any action in reliance on the contents of this information is strictly prohibited.

Whilst reasonable care has been taken to avoid virus transmission, no responsibility for viruses is taken and it is your responsibility to carry out such checks as you feel appropriate.

Saber Razmjooei and Peter Wells trading as Lutra Consulting.

6

Dear Lucien,

here is a little bash script, which updates the attribut-table from
your vector file and
add all extracted values from multiple rasters in a mapset.

for A_Raster in $( eval "g.mlist type=rast mapset="..." )
do
v.db.addcol map="..." columns="$A_Raster double precision" --overwrite
v.what.rast vector="..." raster=$A_Raster column=$A_Raster --overwrite
echo "Spatial query of $A_Raster finished at $(date)"
done

Kind regards
Marco

7

Lucien wrote:

> I have to extract isolated value from several raster.

Markus Metz wrote:

The GRASS module r.what [0] extracts raster values at a
given location. As input you need to give the name of the
raster map and the coordinates as east,north

[0] http://grass.osgeo.org/grass64/manuals/r.what.html

note that r.what can take multiple rasters and multiple point-
values as input, all in one pass, to make a table.

Hamish

8

I have to extract isolated value from several raster.

Interestingly.
Yesterday I had the same task to solve:
have a vector with some sampling points and wanting to insert values
from 2 raster layers to the table.

The GRASS module r.what [0] extracts raster values at a
given location. As input you need to give the name of the
raster map and the coordinates as east,north

[0] http://grass.osgeo.org/grass64/manuals/r.what.html

note that r.what can take multiple rasters and multiple point-
values as input, all in one pass, to make a table.

I do not understand why r.what should be of help here.

r.what [1] just prints out.
Isn't the right tool v.sample [2]?
But this tool works in a way the user would not expect it:

Now, the problem is with v.sample:
1) it just accepts 1 input raster only.
2) it creates a new vector with a totally new table instead of updating
the attribute table of the input vector with the raster values

A note to 2): Such a procedure is implemented in v.rast.stats.

What do the GRASS devs think about the items 1) & 2)?
Is it a legitimate enhancement?

Kind regards,
Timmie

[1] http://grass.osgeo.org/grass64/manuals/r.what.html
[2] http://grass.osgeo.org/grass64/manuals/v.sample.html
[3] http://grass.osgeo.org/grass64/manuals/v.rast.stats.html

9

On Sat, Apr 6, 2013 at 1:11 PM, Tim Michelsen
<timmichelsen@gmx-topmail.de> wrote:

I have to extract isolated value from several raster.

Interestingly.
Yesterday I had the same task to solve:
have a vector with some sampling points and wanting to insert values
from 2 raster layers to the table.

The GRASS module r.what [0] extracts raster values at a
given location. As input you need to give the name of the
raster map and the coordinates as east,north

[0] http://grass.osgeo.org/grass64/manuals/r.what.html

note that r.what can take multiple rasters and multiple point-
values as input, all in one pass, to make a table.

I do not understand why r.what should be of help here.

r.what would not help here, but for the original question because
there a vector map was not mentioned.

r.what [1] just prints out.
Isn't the right tool v.sample [2]?

r.what extracts values from one or more raster maps at locations given
by coordinates.

v.what.rast uploads raster values at positions of vector points to the table.

v.rast.stats calculates univariate statistics from a raster map based
on vector areas (polygons) and uploads statistics to new attribute
columns.

HTH,

Markus M

But this tool works in a way the user would not expect it:

Now, the problem is with v.sample:
1) it just accepts 1 input raster only.
2) it creates a new vector with a totally new table instead of updating
the attribute table of the input vector with the raster values

A note to 2): Such a procedure is implemented in v.rast.stats.

What do the GRASS devs think about the items 1) & 2)?
Is it a legitimate enhancement?

Kind regards,
Timmie

[1] http://grass.osgeo.org/grass64/manuals/r.what.html
[2] http://grass.osgeo.org/grass64/manuals/v.sample.html
[3] http://grass.osgeo.org/grass64/manuals/v.rast.stats.html

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

10

I use a Python script with gdal to do that

If you open a grass raster with gdal:

    from osgeo import gdal
    # grass raster dem10m
    file = '/Users/grassdata/geol/MNT/cellhd/dem10m'
    layer = gdal.Open(file)
    gt =layer.GetGeoTransform()
    bands = layer.RasterCount
    print bands
    1
    print gt
   (263104.72544800001, 10.002079999999999, 0.0, 155223.647811, 0.0,
-10.002079999999999)

As a result, you have the number of bands and the geotransform parameters.
If you want to extract the value of the raster under a xy point:

   x,y = (263220.5,155110.6)
   # transform to raster point coordinates
   rasterx = int((x - gt[0]) / gt[1])
   rastery = int((y - gt[3]) / gt[5])
   # only one band here
   print layer.GetRasterBand(1).ReadAsArray(rasterx,rasterx, 1, 1)
   array([[222]])

With 3 raster bands with the same xy point:

   # recalculate the new raster coordinates px2,py2 (geological map here)
   .... (same as precedent)
   # and
   band1 = layer.GetRasterBand(1)
   band2 = layer.GetRasterBand(2)
   band3 = layer.GetRasterBand(3)
   print band1.ReadAsArray(px2,py2, 1, 1)
   [[253]]
   print band2.ReadAsArray(px2,py2, 1, 1)
   [[215]]
   print band3.ReadAsArray(px2,py2, 1, 1)
   [[118]]

So you could make a function that allows to get the values of multiple
rasters under a xy point:

    def Val_raster(x,y,layer,bands,gt):
        col=
        px = int((x - gt[0]) / gt[1])
        py =int((y - gt[3]) / gt[5])
        for j in range(bands):
            band = layer.GetRasterBand(j+1)
            data = band.ReadAsArray(px,py, 1, 1)
            col.append(data[0][0])
      return col

and

    file1 = '/Users/grassdata/geol/MNT/cellhd/dem10m'
    layer1 = gdal.Open(file1)
    gt1 =layer1.GetGeoTransform()
    bands1 = layer1.RasterCount
    print Val_raster(x,y,layer1, bands1,gt1)
    [222.0]
    file2 = '/Users/grassdata/geol/MNT/cellhd/geolmap'
    layer2 = gdal.Open(file2)
    gt2 =layer2.GetGeoTransform()
    bands2 = layer2.RasterCount
    print Val_raster(x,y,layer2, bands2,gt2)
    [253, 215, 118]

You can even get values from a raster in another MAPSET in the same
LOCATION:

    file3 = '/Users/grassdata/geol/topo_maps/cellhd/HCmap'
    layer3 = gdal.Open(file3)
    gt3 =layer3.GetGeoTransform()
    bands3 = layer3.RasterCount
    print Val_raster(x,y,layer3, bands3,gt3)
    [103, 115, 112]

--
View this message in context: http://osgeo-org.1560.n6.nabble.com/GRASS-6-4-How-to-extract-point-values-from-multiple-rasters-tp5044592p5045037.html
Sent from the Grass - Users mailing list archive at Nabble.com.

11

I use a Python script with gdal to do that

maybe a nice example for the wiki?
(http://grasswiki.osgeo.org/wiki/Main_Page)

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.n6.nabble.com/GRASS-6-4-How-to-extract-point-values-from-multiple-rasters-tp5044592p5045051.html
Sent from the Grass - Users mailing list archive at Nabble.com.

12

r.what [1] just prints out.
> Isn't the right tool v.sample [2]?

r.what extracts values from one or more raster maps at locations given
by coordinates.

v.what.rast uploads raster values at positions of vector points to the table.

But why is there v.sample if all other commands have similar functionality?

13

On Sun, Apr 7, 2013 at 1:23 PM, Tim Michelsen
<timmichelsen@gmx-topmail.de> wrote:

r.what [1] just prints out.
> Isn't the right tool v.sample [2]?

r.what extracts values from one or more raster maps at locations given
by coordinates.

v.what.rast uploads raster values at positions of vector points to the table.

But why is there v.sample if all other commands have similar functionality?

v.sample is useful for cross-validation after surface interpolation:
it samples a raster map with nearest neighbor, bilinear or bicubic,
compares the sampled raster value to a point attribute and writes out
the difference in a new vector map. No other module can do that

Markus M

14

Dear Marco,

Thank you for your answer and your script. I tried it but I meet some problem.

I used the script :

for A_Raster in $( eval "g.mlist type=rast mapset=RFE_test" )
do
v.db.addcol extracting_point_ETo_RFE_test columns="$A_Raster double precision" --overwrite
v.what.rast vector="extracting_point_ETo_RFE_test" raster=$A_Raster column=$A_Raster --overwrite
echo "Spatial query of $A_Raster finished at $(date)"
done

and I receive the following error message for each raster of my mapset :

extract_PointValue_from_RFE_raster.sh: line 13: v.db.addcol: command not found
ERREUR :Colonne <RFE_2001001> non trouvÚe
Spatial query of RFE_2001001 finished at Wed Apr 10 13:54:36 GMT 2013
extract_PointValue_from_RFE_raster.sh: line 13: v.db.addcol: command not found
ERREUR :Colonne <RFE_2001002> non trouvÚe
Spatial query of RFE_2001002 finished at Wed Apr 10 13:54:36 GMT 2013
...

I separated commands (v.db.addcol extracting_point_ETo_RFE_test columns="test double precision") and it worked

Do someone has an idea of the problem?

thanks

Lucien

________________________________________
De : grass-user-bounces@lists.osgeo.org [grass-user-bounces@lists.osgeo.org] de la part de Marco Sciaini [sciaini.marco@gmail.com]
Date d'envoi : jeudi 4 avril 2013 16:58
À : grass-user@lists.osgeo.org
Objet : [GRASS-user] Fwd: GRASS 6.4 How to extract point-values from multiple rasters

Dear Lucien,

here is a little bash script, which updates the attribut-table from
your vector file and
add all extracted values from multiple rasters in a mapset.

for A_Raster in $( eval "g.mlist type=rast mapset="..." )
do
v.db.addcol map="..." columns="$A_Raster double precision" --overwrite
v.what.rast vector="..." raster=$A_Raster column=$A_Raster --overwrite
echo "Spatial query of $A_Raster finished at $(date)"
done

Kind regards
Marco
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user