1

Here´s the problem: I have about 40 shapefiles of contour lines (type line) and I need to interpolate them to create DEM's (tiff, or better some grid format). To do so, some batch scripting would be nice but I just can't figure out a good way.

I could do it with the batch tools of QGIS (GRASS v.bspline would do the job), but I would need to set the parameters for EACH single run. Not very slick at all. So I tried to do it with SAGA CMD in a bash script, but the problem is that the interpolation function needs the input of the area extension and I don't really know how to code that again.

It´s my first real contact with coding for GIS or Python but I've programmed before in Matlab and know some Linux bash stuff so that should not be too much of a problem.

So can anyone help me with a simple way to batch script basically what the Interpolation Plugin of QGIS does?

nmtoken
  • 13,355
  • 5
  • 38
  • 87
martin
  • 321
  • 1
  • 9
  • 1
    Possible duplicate of http://gis.stackexchange.com/questions/17432/how-to-smooth-interpolate-a-raster-in-python-using-gdal – nickves Apr 28 '14 at 18:33
  • 1
    @nickves The duplicate that you cite makes no mention of how to batch process which I think is the focus of this Question. Consequently, I am voting to leave it Open. – PolyGeo Apr 28 '14 at 21:38
  • Check this http://gis.stackexchange.com/questions/18354/batch-create-dem-from-contour-maps-open-source – jdeltoro1973 Apr 29 '14 at 01:18

1 Answers1

2

Here's a (very raw!) script for use in GRASS to process a shapefile and create a DEM tile. Four requirements:

  1. You must be in a GRASS session to run this
  2. You must have a GRASS location setup that matches the CRS of the shapefiles
  3. You must have the resolution set to a reasonable number in GRASS's region settings

  4. The contours shapefile needs an elevation column named "level" here

    #!/usr/bin/env python

    """
Interpolate a dem from a line shapefile
The shapefile name is passed on the command line
It must have a column named "level" with elevation values 
"""

import os, sys
import grass.script as grass

if "GISBASE" not in os.environ:
    print "You must be in GRASS GIS to run this program."
    sys.exit(1)

if len(sys.argv) < 2:
    print 'No shapefile specified'
    sys.exit(1)

# Get shapefile name on command line
shp_path    = sys.argv[1]
shp         = os.path.basename(shp_path)
ctour_orig  = os.path.splitext(shp)[0].lower()
ctour       = ctour_orig.replace(' ','_')
ctour_rast  = ctour + '_rast'
dem_tile    = ctour + '_dem'

print "Processing shapefile: %s" % ctour

grass.run_command('v.in.ogr', dsn=shp_path, output=ctour, type='line', overwrite=True)
grass.run_command('g.region', flags='p', vect=ctour, overwrite=True)
grass.run_command('v.to.rast', input=ctour, output=ctour, type='line', use='attr', attrcol='level', overwrite=True)
grass.run_command('r.surf.contour', input=ctour, output=dem_tile, overwrite=True)

print "Completed creating raster tile: %s" % dem_tile

If you save the script as "shp_to_dem.py" then you could loop thru the the shapefiles in a directory and pass each to this script with something like (on Linux):

for s in *.shp; do python shp_to_dem.py $s; done

Maybe this will get you started.

Micha
  • 15,555
  • 23
  • 29