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?

3
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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