Skip to content

Instantly share code, notes, and snippets.

@joshdoe
Created September 1, 2011 15:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save joshdoe/1186390 to your computer and use it in GitHub Desktop.
Save joshdoe/1186390 to your computer and use it in GitHub Desktop.
GRASS script to convert cost points generated by pgRouting's driving_distance() to polygons of 1/4 mile increments
#!/usr/bin/env python
import os
import atexit
import sys
import subprocess
gisbase = os.environ['GISBASE'] = '/usr/lib/grass64'
gisdbase = os.path.join(os.environ['HOME'], 'grassdata')
location = 'fairfax'
mapset = 'PERMANENT'
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
from grass.script import core as grass
import grass.script.setup as gsetup
gsetup.init(gisbase, gisdbase, location, mapset)
print grass.gisenv()
def main():
cost_points_file = '/home/user/points_2283.shp'
out_file = '/home/user/polys.shp'
# generate contour polygons
grass.run_command('v.in.ogr', dsn=cost_points_file, output='cost_points',
overwrite=True)
grass.run_command('g.region', vect='cost_points', rows=512, cols=512)
grass.run_command('v.surf.rst', input='cost_points', elev='interp',
zcolumn='cost', overwrite=True)
grass.run_command('r.mapcalc', result='int(interp/0.4)')
grass.run_command('r.to.vect', flags='s', input='result',
output='polys', feature='area', overwrite=True)
grass.run_command('v.out.ogr', input='polys', dsn=out_file, type='area',
flags='c')
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment