[GRASS-SVN] r64982 - grass-addons/grass7/vector/v.isochrones
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Apr 4 06:39:00 PDT 2015
Author: mlennert
Date: 2015-04-04 06:39:00 -0700 (Sat, 04 Apr 2015)
New Revision: 64982
Added:
grass-addons/grass7/vector/v.isochrones/v_isochrones_cost.png
grass-addons/grass7/vector/v.isochrones/v_isochrones_net.png
Removed:
grass-addons/grass7/vector/v.isochrones/v_isochrones.png
Modified:
grass-addons/grass7/vector/v.isochrones/v.isochrones.html
grass-addons/grass7/vector/v.isochrones/v.isochrones.py
Log:
Added a second algorithm using v.net.iso
Modified: grass-addons/grass7/vector/v.isochrones/v.isochrones.html
===================================================================
--- grass-addons/grass7/vector/v.isochrones/v.isochrones.html 2015-04-04 04:09:28 UTC (rev 64981)
+++ grass-addons/grass7/vector/v.isochrones/v.isochrones.html 2015-04-04 13:39:00 UTC (rev 64982)
@@ -2,49 +2,80 @@
<em>v.isochrones</em> creates a vector polygon map of isochrones
(<b>isochrones</b>) based on a roads map (<b>map</b>) with speed
-attribute (<b>speed_column</b>), one or several starting points
+or cost attribute (<b>cost_column</b>), one or several starting points
(<b>start_points</b>) and time steps for the isochrones
-(<b>time_steps</b>). Optionally, the module can create a raster map
-with continuous time from the starting points (<b>timemap</b>).
+(<b>time_steps</b>). The module is actually a front-end to different
+GRASS GIS modules, and the user can chose the approach with the <b>method</b> option
+(see the notes section for details).
-<p>
-The user can define the speed to use in offroad areas
-(<b>offroad_speed</b>) and the time unit of the output and the time
-steps (<b>unit</b>).
+
<h2>NOTES</h2>
-The module is a front-end to <em><a href="r.cost.html">r.cost</a></em>.
-It transforms the roads to raster, calculates time cost from the
-starting points and then transforms the result into discrete vector
-polygons based on the time steps chosen.
+Two approaches are currently implemented in the module:
<p>
-Current region settings are used to define the maximal extension and
-the resolution at which the cost map is calculated.
+The first (default) approach is based on <em><a href="v.net.iso">v.net.iso</a></em>. The output
+of that module is then used to assign isochrone values to all pixels based on the nearest
+road segment that is within a given distance (<b>max_distance</b>). The <b>-i</b> flag allows
+to calculate a separate isochrone map for each starting point.
+For this approach, the input road map currently has to be prepared by adding all nodes (using
+<em><a href="v.net">v.net</a></em> with 'operation=nodes' and the '-c' flag. In addition
+and a <b>cost_column</b> has to be present in the attribute table containing for each line segment
+the time of traversal in minutes (i.e. length/speed).
+<p>
+The second approach is based on <em><a href="r.cost.html">r.cost</a></em>.
+It transforms the roads to raster (here <b>cost_column</b> has to point to an attribute
+column containing speed in km/h), assigning a user chosen <b>offroad_speed</b> to all offroad
+pixels, calculates time cost from the starting points and then transforms the result into
+discrete vector polygons based on the time steps chosen. Optionally, the user can chose to
+keep the a raster output of <em><a href="r.cost.html">r.cost</a></em> (<b>timemap</b>).
+Current region settings are used to define the maximal extension and the resolution
+at which the cost map is calculated. The <b>memory</b> option allows to decide how much
+memory the <em><a href="r.cost.html">r.cost</a></em> module can use (see the r.cost man page
+for more details).
+
<h2>EXAMPLE</h2>
<div class="code"><pre>
-g.copy vector=roadsmajor,roads
-v.db.addcolumn roads col="speed int"
-v.db.update roads col=speed value=120 where="ROAD_NAME='US-1'"
-v.db.update roads col=speed value=90 where="(ROAD_NAME like 'US%' AND ROAD_NAME <> 'US-1') OR ROAD_NAME like 'I-%'"
-v.db.update roads col=speed value=50 where="speed is null"
+v.net -c input=roadsmajor operation=nodes output=myroads
+v.db.addcolumn myroads col="speed int"
+v.db.addcolumn myroads col="length double precision"
+v.db.addcolumn myroads col="cost double precision"
+v.db.update myroads col=speed value=120 where="ROAD_NAME='US-1'"
+v.db.update myroads col=speed value=90 where="(ROAD_NAME like 'US%' AND ROAD_NAME <> 'US-1') OR ROAD_NAME like 'I-%'"
+v.db.update myroads col=speed value=50 where="speed is null"
+v.to.db myroads op=length col=length
+v.db.update myroads col=cost qcol="(length/(speed*1000))*60"
-g.region vector=roads res=100 -a -p
+g.region vector=myroads res=50 -a -p
echo "634637|224663" | v.in.ascii input=- output=start x=1 y=2
-v.isochrones map=roads speed_col=speed start_points=start isochrones=isochrones30_60_120 time_steps=30,60,120
+v.isochrones map=myroads cost_column=cost start_points=start isochrones=isochrones_vnetiso time_steps=15,30,60 \
+ max_distance=1000 method=v.net.iso
+
</pre></div>
<center>
- <img src="v_isochrones.png" border="1"><br>
- 30, 60 and 120 minute isochrones from the start point
+ <img src="v_isochrones_net.png" border="1"><br>
+ 15, 30 and 60 minute isochrones from the start point
+ using method v.net.iso
</center>
+v.isochrones map=myroads cost_column=speed start_points=start isochrones=isochrones_rcost time_steps=15,30,60 \
+ method=r.cost memory=1000
+</pre></div>
+
+<center>
+ <img src="v_isochrones_cost.png" border="1"><br>
+ 15, 30 and 60 minute isochrones from the start point
+ using method r.cost
+</center>
+
+
<h2>SEE ALSO</h2>
<em>
Modified: grass-addons/grass7/vector/v.isochrones/v.isochrones.py
===================================================================
--- grass-addons/grass7/vector/v.isochrones/v.isochrones.py 2015-04-04 04:09:28 UTC (rev 64981)
+++ grass-addons/grass7/vector/v.isochrones/v.isochrones.py 2015-04-04 13:39:00 UTC (rev 64982)
@@ -23,11 +23,20 @@
#% required: yes
#%end
#%option G_OPT_V_FIELD
+#% key: roads_layer
+#% label: Layer number of the roads with relevant attributes
#% required: yes
#%end
+#%option G_OPT_V_FIELD
+#% key: node_layer
+#% label: Layer number of the nodes on the network (for method=v.net.iso)
+#% required: no
+#% answer: 2
+#% guisection: v.net.iso
+#%end
#%option G_OPT_DB_COLUMN
-#% key: speed_column
-#% description: Name of speed attribute column (in km/h)
+#% key: cost_column
+#% description: Name of speed attribute column (in km/h) or cost column (in minutes)
#% required: yes
#%end
#%option G_OPT_V_INPUT
@@ -37,174 +46,348 @@
#%end
#%option G_OPT_V_OUTPUT
#% key: isochrones
-#% label: Output vector map with isochrone polygons
+#% label: Output vector map with isochrone polygons (Output prefix with flag -i)
#% required: yes
#%end
+#%option
+#% key: time_steps
+#% type: double
+#% description: Time steps of isochrones (in units defined by unit parameter)
+#% multiple: yes
+#% required: yes
+#%end
#%option G_OPT_R_OUTPUT
#% key: timemap
#% label: Optional output raster map with continuous time from starting points
#% required: no
+#% guisection: r.cost
#%end
#%option
-#% key: time_steps
+#% key: offroad_speed
#% type: double
-#% description: Time steps of isochrones (in units defined by unit parameter)
-#% multiple: yes
-#% required: yes
+#% description: Speed for off-road areas (in km/h > 0)
+#% required: no
+#% answer: 5.0
+#% guisection: r.cost
#%end
#%option
-#% key: offroad_speed
+#% key: memory
+#% description: Amount of memory (in MB) use
#% type: integer
-#% description: Speed for off-road areas (in km/h)
#% required: no
-#% answer: 5
+#% answer: 300
+#% guisection: r.cost
#%end
#%option
-#% key: unit
-#% description: Time unit to use (seconds, minutes, hours)
-#% type: string
+#% key: max_distance
+#% type: double
+#% description: Maximum distance (m) from network to include into isochrones
#% required: no
-#% options: s,m,h
-#% answer: m
+#% guisection: v.net.iso
#%end
+#%option
+#% key: method
+#% description: Method to use for isochrone calculation
+#% type: string
+#% required: yes
+#% options: v.net.iso,r.cost
+#% answer: v.net.iso
+#%end
+#%flag
+#% key: i
+#% description: Create individual isochrone map for each starting point
+#% guisection: v.net.iso
+#%end
+
+
import os
import atexit
import math
import grass.script as grass
-tmp_cost_map = None
-tmp_map = None
-tmp_time_map = None
-tmp_region_map = None
+global isoraw
+global isos_extract
+global isos_extract_rast
+global isos_grow_cat
+global isos_grow_cat_int
+global isos_grow_distance
+global isos_grow_distance_recode
+global isos_poly_all
+global isos_final
+isoraw = None
+isos_extract = None
+isos_extract_rast = None
+isos_grow_cat = None
+isos_grow_cat_int = None
+isos_grow_distance = None
+isos_grow_distance_recode = None
+isos_poly_all = None
+isos_final = None
+
def cleanup():
- # remove temporary cost column from road map
- grass.run_command('v.db.dropcolumn',
- map=roads,
- column=cost_column,
- quiet=True)
+ if method == 'r.cost':
+ # remove temporary cost column from road map
+ if grass.vector_db(roads)[int(layer)]['driver'] == 'dbf':
+ grass.run_command('v.db.dropcolumn',
+ layer=layer,
+ map=roads,
+ column=tmp_cost_column,
+ quiet=True)
- if grass.find_file(tmp_cost_map, element='raster')['name']:
- grass.run_command('g.remove', flags='f', type='raster', name=tmp_cost_map, quiet=True)
- if grass.find_file(tmp_map, element='raster')['name']:
- grass.run_command('g.remove', flags='f', type='raster', name=tmp_map, quiet=True)
- if grass.find_file(tmp_time_map, element='raster')['name']:
- grass.run_command('g.remove', flags='f', type='raster', name=tmp_time_map, quiet=True)
- if grass.find_file(tmp_region_map, element='raster')['name']:
- grass.run_command('g.remove', flags='f', type='raster', name=tmp_region_map, quiet=True)
+ if grass.find_file(tmp_cost_map, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=tmp_cost_map, quiet=True)
+ if grass.find_file(tmp_time_map, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=tmp_time_map, quiet=True)
+ if grass.find_file(tmp_region_map, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=tmp_region_map, quiet=True)
+ elif method == 'v.net.iso':
+ if grass.find_file(isoraw, element='vector')['name']:
+ grass.run_command('g.remove', flags='f', type='vector', name=isoraw, quiet=True)
+ if grass.find_file(isos_extract, element='vector')['name']:
+ grass.run_command('g.remove', flags='f', type='vector', name=isos_extract, quiet=True)
+ if grass.find_file(isos_extract_rast, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=isos_extract_rast, quiet=True)
+ if grass.find_file(isos_grow_cat, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=isos_grow_cat, quiet=True)
+ if grass.find_file(isos_grow_cat_int, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster',
+ name=isos_grow_cat_int, quiet=True)
+ if grass.find_file(isos_grow_distance, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=isos_grow_distance, quiet=True)
+ if grass.find_file(isos_grow_distance_recode, element='cell')['name']:
+ grass.run_command('g.remove', flags='f', type='raster', name=isos_grow_distance_recode, quiet=True)
+ if grass.find_file(isos_poly_all, element='vector')['name']:
+ grass.run_command('g.remove', flags='f', type='vector', name=isos_poly_all, quiet=True)
+ if grass.find_file(isos_final, element='vector')['name']:
+ grass.run_command('g.remove', flags='f', type='vector', name=isos_final, quiet=True)
+
+def isocalc(isoraw):
+
+ global isos_extract
+ global isos_extract_rast
+ global isos_grow_cat
+ global isos_grow_cat_int
+ global isos_grow_distance
+ global isos_grow_distance_recode
+ global isos_poly_all
+
+ isos_extract = 'isos_extract_%d' % os.getpid()
+ isos_extract_rast = 'isos_extract_rast_%d' % os.getpid()
+ isos_grow_cat = 'isos_grow_cat_%d' % os.getpid()
+ isos_grow_cat_int = 'isos_grow_cat_int_%d' % os.getpid()
+ isos_grow_distance = 'isos_grow_distance_%d' % os.getpid()
+ isos_grow_distance_recode = 'isos_grow_distance_recode_%d' % os.getpid()
+ isos_poly_all = 'isos_poly_all_%d' % os.getpid()
+
+ grass.use_temp_region()
+
+ grass.run_command('v.extract', input_=isoraw,
+ cat=output_cats[0:-1], output=isos_extract,
+ overwrite=True)
+ grass.run_command('g.region', vect=isos_extract,
+ flags='a')
+ grass.run_command('v.to.rast', input_=isos_extract,
+ use='cat', out=isos_extract_rast, overwrite=True)
+ grass.run_command('r.grow.distance', input=isos_extract_rast,
+ value=isos_grow_cat, distance=isos_grow_distance,
+ flags='m', overwrite=True)
+ grass.mapcalc(isos_grow_cat_int + ' = int(' + isos_grow_cat + ')')
+ if max_distance:
+ recode_str = '0:%f:1\n%f:*:0' % (max_distance, max_distance)
+ grass.write_command('r.recode', input_=isos_grow_distance,
+ output=isos_grow_distance_recode, rules='-',
+ stdin=recode_str, overwrite=True)
+ grass.run_command('r.mask', raster=isos_grow_distance_recode,
+ maskcats=1, overwrite=True)
+ grass.run_command('r.to.vect', input_=isos_grow_cat_int,
+ output=isos_poly_all, type_='area', flags='sv',
+ overwrite=True)
+ grass.run_command('v.extract', input_=isos_poly_all,
+ output=isos_final, cats=output_cats[0:-1])
+ if max_distance:
+ grass.run_command('r.mask', flags='r')
+
+
+
def main():
- # Input
+ # Input for all methods
global roads
roads = options['map']
- speed_column = options['speed_column']
+ global layer
+ layer = options['roads_layer']
+ node_layer = options['node_layer']
+ cost_column = options['cost_column']
start_points = options['start_points']
time_steps = options['time_steps'].split(',')
- offroad_speed = int(options['offroad_speed'])
- unit = options['unit']
- # Output
isochrones = options['isochrones']
- if options['timemap']:
- timemap = options['timemap']
- else:
- timemap = None
+ global method
+ method = options['method']
- global tmp_cost_map
- global tmp_map
- global tmp_time_map
- global tmp_region_map
- global cost_column
+ if method == 'r.cost':
+ offroad_speed = float(options['offroad_speed'])
+ if offroad_speed == 0:
+ grass.message(_('Offroad speed has to be > 0. Set to 0.00001.'))
+ offroad_speed = 0.00001
- tmp_cost_map = 'cost_map_tmp_%d' % os.getpid()
- tmp_map = 'map_tmp_%d' % os.getpid()
- tmp_time_map = 'time_map_tmp_%d' % os.getpid()
- tmp_region_map = 'region_map_tmp_%d' % os.getpid()
+ memory = int(options['memory'])
+ # Output
+ if options['timemap']:
+ timemap = options['timemap']
+ else:
+ timemap = None
- # get current resolution
- region = grass.region()
- resolution = math.sqrt(float(region['nsres']) * float(region['ewres']))
+ global tmp_cost_map
+ global tmp_time_map
+ global tmp_region_map
+ global tmp_cost_column
- # add cost column to road vector
- cost_column = 'tmp%d' % os.getpid()
- def_cost_column = cost_column + ' DOUBLE PRECISION'
- grass.run_command('v.db.addcolumn',
- map=roads,
- column=def_cost_column,
- quiet=True)
+ tmp_cost_map = 'cost_map_tmp_%d' % os.getpid()
+ tmp_time_map = 'time_map_tmp_%d' % os.getpid()
+ tmp_region_map = 'region_map_tmp_%d' % os.getpid()
- # calculate cost (in seconds) depending on speed
- # (resolution/(speed (in km/h) * 1000 / 3600))
- query_value = "%s / (%s * 1000 / 3600)" % (resolution, speed_column)
- grass.run_command('v.db.update',
- map=roads,
- column=cost_column,
- qcolumn=query_value)
+ # get current resolution
+ region = grass.region()
+ resolution = math.sqrt(float(region['nsres']) * float(region['ewres']))
- # transform to raster
- grass.run_command('v.to.rast',
- input=roads,
- out=tmp_cost_map,
- use='attr',
- attrcolumn=cost_column,
- type='line')
+ if grass.vector_db(roads)[int(layer)]['driver'] == 'dbf':
+ # add cost column to road vector
+ tmp_cost_column = 'tmp%d' % os.getpid()
+ def_cost_column = tmp_cost_column + ' DOUBLE PRECISION'
+ grass.run_command('v.db.addcolumn',
+ map=roads,
+ layer=layer,
+ column=def_cost_column,
+ quiet=True)
- # replace null values with cost for off-road areas
- # (resolution/(off-road speed * 1000 / 3600))
- null_value = resolution / (float(offroad_speed) * 1000 / 3600)
- grass.run_command('r.null', map=tmp_cost_map, null=null_value)
+ # calculate cost (in minutes) depending on speed
+ # (resolution/(speed (in km/h) * 1000 / 60))
+ query_value = "%s / (%s * 1000 / 60)" % (resolution, cost_column)
+ grass.run_command('v.db.update',
+ map=roads,
+ column=tmp_cost_column,
+ qcolumn=query_value)
+ else:
+ tmp_cost_column = "%s / (%s * 1000 / 60)" % (resolution, cost_column)
- # calculate time distance from starting points
- grass.run_command('r.cost',
- input=tmp_cost_map,
- start_points=start_points,
- output=tmp_map)
+ # transform to raster
+ grass.run_command('v.to.rast',
+ input=roads,
+ out=tmp_cost_map,
+ use='attr',
+ attrcolumn=tmp_cost_column,
+ type='line',
+ memory=memory)
- # if necessary, transform time to desired time units
- if unit == 'm':
- expression = tmp_time_map + ' = ' + tmp_map + ' / 60'
- grass.mapcalc(expression)
- elif unit == 'h':
- expression = tmp_time_map + ' = ' + tmp_map + ' / 3600'
- grass.mapcalc(expression)
- else:
- grass.run_command('g.rename', raster = (tmp_map, tmp_time_map))
+ # replace null values with cost for off-road areas
+ # (resolution/(off-road speed * 1000 / 60))
+ null_value = resolution / (offroad_speed * 1000 / 60)
+ grass.run_command('r.null', map=tmp_cost_map, null=null_value)
- if timemap:
- grass.run_command('g.copy', raster=(tmp_time_map, timemap))
- grass.run_command('r.colors', map=timemap, color='grey', flags='ne')
+ # limit the cumulated cost surface calculation to the max time distance
+ # requested
+ max_cost = time_steps[-1]
- # recode time distance to time steps
- recode_rules = '0:%s:%s\n' % (time_steps[0], time_steps[0])
- for count in range(1, len(time_steps)):
- recode_rules += time_steps[count-1] + ':'
- recode_rules += time_steps[count] + ':'
- recode_rules += time_steps[count] + '\n'
+ # calculate time distance from starting points
+ grass.run_command('r.cost',
+ input=tmp_cost_map,
+ start_points=start_points,
+ output=tmp_time_map,
+ max_cost=max_cost,
+ memory=memory)
- grass.write_command('r.recode',
- input=tmp_time_map,
- output=tmp_region_map,
- rules='-',
- stdin=recode_rules)
+ if timemap:
+ grass.run_command('g.copy', raster=(tmp_time_map, timemap))
+ grass.run_command('r.colors', map=timemap, color='grey', flags='ne')
- # transform to vector areas
- grass.run_command('r.to.vect',
- input=tmp_region_map,
- output=isochrones,
- type='area',
- column='traveltime')
+ # recode time distance to time steps
+ recode_rules = '0:%s:%s\n' % (time_steps[0], time_steps[0])
+ for count in range(1, len(time_steps)):
+ recode_rules += time_steps[count-1] + ':'
+ recode_rules += time_steps[count] + ':'
+ recode_rules += time_steps[count] + '\n'
- # give the polygons a default color table
- grass.run_command('v.colors',
- map=isochrones,
- use='attr',
- column='traveltime',
- color='grey')
+ grass.write_command('r.recode',
+ input=tmp_time_map,
+ output=tmp_region_map,
+ rules='-',
+ stdin=recode_rules)
+ # transform to vector areas
+ grass.run_command('r.to.vect',
+ input=tmp_region_map,
+ output=isochrones,
+ type='area',
+ column='traveltime',
+ flags='s')
+
+ # give the polygons a default color table
+ grass.run_command('v.colors',
+ map=isochrones,
+ use='attr',
+ column='traveltime',
+ color='grey')
+
+ elif method == 'v.net.iso':
+ global max_distance
+ if(options['max_distance']):
+ max_distance = float(options['max_distance'])
+ else:
+ max_distance = None
+ global output_cats
+ output_cats = []
+ for i in range(1,len(time_steps)+2):
+ output_cats.append(i)
+ startpoints=grass.read_command('v.distance', from_=start_points,
+ to=roads, to_type='point', to_layer=node_layer,
+ upload='cat', flags='p', quiet=True).split('\n')[1:-1]
+
+ global isoraw
+ isoraw = 'isoraw_temp_%d' % os.getpid()
+ global isos_final
+ isos_final = 'isos_final_%d' % os.getpid()
+
+ if flags['i']:
+ for point in startpoints:
+ startpoint_cat = point.split('|')[0]
+ startnode_cat = point.split('|')[1]
+ grass.run_command('v.net.iso', input_=roads, output=isoraw,
+ center_cats=startnode_cat, costs=time_steps,
+ arc_column=cost_column, overwrite=True)
+ isocalc(isoraw)
+ outname = isochrones + '_' + startpoint_cat
+ grass.run_command('g.rename', vect=isos_final+','+outname)
+ # give the polygons a default color table
+ grass.run_command('v.colors',
+ map=outname,
+ use='cat',
+ color='grey')
+
+ else:
+ startnodes = []
+ for point in startpoints:
+ startnodes.append(point.split('|')[1])
+ grass.run_command('v.net.iso', input_=roads, output=isoraw,
+ center_cats=startnodes, costs=time_steps,
+ arc_column=cost_column, overwrite=True)
+ isocalc(isoraw)
+ grass.run_command('g.rename', vect=isos_final+','+isochrones)
+ # give the polygons a default color table
+ grass.run_command('v.colors',
+ map=isochrones,
+ use='cat',
+ color='grey')
+
+ else:
+ grass.fatal(_("You need to chose at least one of the methods"))
+
+
if __name__ == "__main__":
options, flags = grass.parser()
atexit.register(cleanup)
Deleted: grass-addons/grass7/vector/v.isochrones/v_isochrones.png
===================================================================
(Binary files differ)
Added: grass-addons/grass7/vector/v.isochrones/v_isochrones_cost.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/vector/v.isochrones/v_isochrones_cost.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/vector/v.isochrones/v_isochrones_net.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/vector/v.isochrones/v_isochrones_net.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
More information about the grass-commit
mailing list