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

[gene]

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.