Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question

1 Answer 1

up vote 2 down vote accepted

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...
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.