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

gene martin.laloux at gmail.com
Sun Apr 7 01:13:24 PDT 2013

```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) / gt)
rastery = int((y - gt) / gt)
# only one band here
array([])

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)
[]
[]
[]

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) / gt)
py =int((y - gt) / gt)
for j in range(bands):
band = layer.GetRasterBand(j+1)
col.append(data)
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.
```