Programmatic raster-vector calculation

Question

I currently use GRASS' r.mapcalc to find the difference between vector and raster data using a loop and temporary files and layers.

Is it possible to do the core calculation between raster with a vector data more directly either within GRASS, in some external tool, or some python script?

#!/bin/bash 
TOPO=mytopo
g.region -p rast=mytopo res=10
for file in INPUTS/* ; do
  imbase=${file#INPUTS/}   # strip off the directory prefix
  v.in.ogr input=$file output=temp_a --overwrite
  v.to.rast in=temp_a out=temp_b use=attr column=elev --overwrite
  r.mapcalc "temp_c" = "(if( "temp_b" <= "$TOPO" ,null(),if( "$TOPO" <0,0, "temp_b" - "$TOPO" )))"
   r.out.gdal input=temp_c -c  format=GTiff type=Float32 output=OUTPUTS/$imbase
end
g.remove temp_a
g.remove temp_b
g.remove temp_c

I like GRASS' the ability to do the mapcalc on resolutions and extents other than my source data. but although this works, the process seems seems awkward and slow for numbers of files in the hundreds. Ideally, I'd like something like the below:

 gdalogr_diff.py mytopo.tif PG:dbname=warmerda -sql "SELECT elev[${timeslice}] from temp_output" -f GTiff out_${timeslice}.tif

Answer 1

You could rasterise your vector layer in python, and then apply you condition and save the result:

import gdal
import ogr
import numpy as np
def rasterize_vector ( vector_fname, nx, ny, geo_T, srs, the_attribute ):
    """vector fname is a OGR-friendy database/shapefile filename. nx ad ny
    are the target dataset rows and columns. geoT and srs
    are the geotransform and projection from the target dataset. the_attribute
    is the attribute to rasterize in the vector dataset."""
    # Open the data source
    vector_source = ogr.Open( vector_fname )
    features = vector_source.GetLayer ( 0 )
    source_layer = source_ds.GetLayer(0)
    target_ds = gdal.GetDriverByName( 'MEM' ).Create( "", nx, ny, 1,\
        gdal.GDT_Int32)
    target_ds.SetGeoTransform( geo_T )
    target_ds.SetProjection( srs )
    # Rasterise!
    err = gdal.RasterizeLayer(target_ds, [1], source_layer,
        options=["ATTRIBUTE=%s" % the_attribute ])
    if err != 0:
        raise Exception("error rasterizing layer: %s" % err)
    data = target_ds.ReadAsArray()
    return data

if __name__ == "__main__":
    topo = gdal.Open ( "mytopo.tif" )
    geo_t = topo.GetGeoTransform()
    srs = topo.GetSpatialRef()
    nx = topo.RasterXSize
    ny = topo.RasterYSize
    elev = rasterize_vector ( "my_ogr.shp", nx, ny, geo_T, srs, "elev" )
    # Do your processing here
    # Save output...