[GRASS-SVN] r61772 - in grass-addons/grass7/vector: . v.centerline
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Aug 29 09:59:07 PDT 2014
Author: mlennert
Date: 2014-08-29 09:59:07 -0700 (Fri, 29 Aug 2014)
New Revision: 61772
Added:
grass-addons/grass7/vector/v.centerline/
grass-addons/grass7/vector/v.centerline/Makefile
grass-addons/grass7/vector/v.centerline/centerline_flightpaths.png
grass-addons/grass7/vector/v.centerline/centerline_river.png
grass-addons/grass7/vector/v.centerline/v.centerline.html
grass-addons/grass7/vector/v.centerline/v.centerline.py
Log:
new module for calculating centerline of a series of lines
Added: grass-addons/grass7/vector/v.centerline/Makefile
===================================================================
--- grass-addons/grass7/vector/v.centerline/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.centerline/Makefile 2014-08-29 16:59:07 UTC (rev 61772)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.centerline
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/vector/v.centerline/centerline_flightpaths.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/vector/v.centerline/centerline_flightpaths.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/vector/v.centerline/centerline_river.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/vector/v.centerline/centerline_river.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/vector/v.centerline/v.centerline.html
===================================================================
--- grass-addons/grass7/vector/v.centerline/v.centerline.html (rev 0)
+++ grass-addons/grass7/vector/v.centerline/v.centerline.html 2014-08-29 16:59:07 UTC (rev 61772)
@@ -0,0 +1,87 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.centerline.closest</em> creates a new map with a line representing an approximation of the central tendency of a series of input lines. This cani for example, be the central line of a river represented by its two sides, or a line representing the general direction of a series of flight paths, etc.
+
+Two algorithms are proposed in the module, both based on the idea of using a reference line, creating a series of reference points along this line and then finding the coordinates of corresponding points on all the input lines. The default algorithm uses closest distance to identify corresponding points, while the second algorithm (<b>t</b> flag) draws perpendicular transversals at the reference points and uses the intersections of these transversals with the other lines as corresponding points.
+
+In detail, the default algorithm goes as follows:
+<ul>
+ <li>choose one of the input lines as reference line</li>
+ <li>create a series of points at regular intervals on this line</li>
+ <li>for each of these points:
+ <ul>
+ <li>find the closest point on all input lines</li>
+ <li>get the coordinates of those points</li>
+ <li>calculate the mean or (mathematical) median of these coordinates</li>
+ </ul>
+ </li>
+ <li>use the calculated means (or medians) as vertices of the new line</li>
+</ul>
+<br>
+
+
+The transversals algorithm goes as follows:
+<ul>
+ <li>choose one of the input lines as reference line</li>
+ <li>create a series of perpendicular (transversal) lines at regular intervals on this line</li>
+ <li>for each of these transversals:
+ <ul>
+ <li>find the intersection points of the transversal with all input lines</li>
+ <li>get the coordinates of those points</li>
+ <li>calculate the mean or (mathematical) median of these coordinates</li>
+ </ul>
+ </li>
+ <li>use the calculated means (or medians) as vertices of the new line</li>
+</ul>
+
+The user can change three parameters in the algorithms: the choice of the reference line (<b>refline</b>), the number of vertices to calculate (<b>vertices</b>) and the search range (<b>range</b>), i.e. the maximum distance of corresponding points for the default algorithm and the length of the transversals on each side of the reference line.
+
+If no reference line is given the module choses the reference line by determining the mean distance of the midpoint of each line to the midpoints of all other lines. The line with the lowest mean distance is then chosen as the reference line. If no range is given, the module uses the mean of the above mean distances as the range for the transversals algorithm, and an unlimited search range for the default algorithm.
+
+If the <b>m</b> flag is set and there are more than 2 lines in the input file, the module calculates the mathematical median of the x and of the y coordinates.
+
+<h2>NOTES</h2>
+
+The median in this module is <b>not</b> the geometric median, but the simple mathematical median respectively of the x and the y coordinates.
+
+The transversals algorithm is very sensitive to the range parameter. The user might want to play around with this parameter to find the best value.
+
+Increasing the number of vertices should have a smoothing effect on the resulting line, but in the case of the transversals algorithm it can possibly lead to more instability.
+
+<h2>EXAMPLE</h2>
+
+<div class="code"><pre>
+v.centerline input=flightpaths output=center_line_mean
+v.centerline -m input=flightpaths output=center_line_median
+v.centerline input=flightpaths output=center_line_mean_5000 range=5000
+v.centerline -t input=flightpaths output=center_line_mean_t
+v.centerline -t input=flightpaths output=center_line_mean_t_8000 range=8000
+</pre></div>
+
+<center>
+ <img src="centerline_flightpaths.png" border="1"><br>
+ Different centerlines resulting from variations in the parameters and flags
+</center>
+
+<div class="code"><pre>
+v.centerline input=river output=center_line
+v.centerline -t input=river output=river_center_t
+</pre></div>
+
+<center>
+ <img src="centerline_river.png" border="1"><br>
+ Mean central line (median only makes sense if number of lines > 2) for distance (red) and transversals (blue) algorithms, the latter with automatically determined range
+</center>
+
+
+<h2>SEE ALSO</h2>
+
+Similar addons:
+<em>
+<a href="v.centerpoint">v.centerpoint</a>
+</em>
+
+<h2>AUTHOR</h2>
+ Moritz Lennert
+
+<p><i>Last changed: $Date$</i>
Added: grass-addons/grass7/vector/v.centerline/v.centerline.py
===================================================================
--- grass-addons/grass7/vector/v.centerline/v.centerline.py (rev 0)
+++ grass-addons/grass7/vector/v.centerline/v.centerline.py 2014-08-29 16:59:07 UTC (rev 61772)
@@ -0,0 +1,319 @@
+#!/usr/bin/env python
+############################################################################
+#
+# MODULE: v.centerline
+# AUTHOR: Moritz Lennert
+# PURPOSE: Takes a map of vector lines and creates a new map containing
+# a central line
+#
+# COPYRIGHT: (c) 2014 Moritz Lennert, and the GRASS Development Team
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+
+#%module
+#% description: Creates a central line of a map of lines
+#% keywords: vector
+#% keywords: lines central
+#%end
+#%option G_OPT_V_INPUT
+#%end
+#%option G_OPT_V_OUTPUT
+#%end
+#%option
+#% key: range
+#% type: double
+#% description: Distance (in map units) of search radius
+#% required: no
+#%end
+#%option
+#% key: refline
+#% type: integer
+#% description: Category of the line to use as initial reference
+#% required: no
+#%end
+#%option
+#% key: vertices
+#% description: Number of vertices for center line
+#% type: integer
+#% answer: 20
+#% required: no
+#%end
+#%flag
+#% key: t
+#% description: Use transversal line algorithm
+#%end
+#%flag
+#% key: m
+#% description: Use approximate median instead of mean
+#%end
+
+import os
+import atexit
+import grass.script as grass
+from grass.pygrass.vector import VectorTopo
+from grass.pygrass.vector import geometry as geo
+
+tmp_points_map = None
+tmp_line_map = None
+tmp_centerpoints_map = None
+
+def cleanup():
+ if grass.find_file(tmp_points_map, element='vector')['name']:
+ grass.run_command('g.remove', vect=tmp_points_map, quiet=True)
+ if grass.find_file(tmp_line_map, element='vector')['name']:
+ grass.run_command('g.remove', vect=tmp_line_map, quiet=True)
+ if grass.find_file(tmp_cleaned_map, element='vector')['name']:
+ grass.run_command('g.remove', vect=tmp_cleaned_map, quiet=True)
+ if grass.find_file(tmp_centerpoints_map, element='vector')['name']:
+ grass.run_command('g.remove', vect=tmp_centerpoints_map, quiet=True)
+ if grass.find_file(tmp_map, element='vector')['name']:
+ grass.run_command('g.remove', vect=tmp_map, quiet=True)
+
+
+def main():
+
+ input = options['input']
+ if options['refline']:
+ refline_cat = int(options['refline'])
+ else:
+ refline_cat = None
+ nb_vertices = int(options['vertices'])
+ if options['range']:
+ search_range = float(options['range'])
+ else:
+ search_range = None
+ output = options['output']
+ transversals = flags['t']
+ median = flags['m']
+
+ global tmp_points_map
+ global tmp_centerpoints_map
+ global tmp_line_map
+ global tmp_cleaned_map
+ global tmp_map
+ tmp_points_map = 'points_map_tmp_%d' % os.getpid()
+ tmp_centerpoints_map = 'centerpoints_map_tmp_%d' % os.getpid()
+ tmp_line_map = 'line_map_tmp_%d' % os.getpid()
+ tmp_cleaned_map = 'cleaned_map_tmp_%d' % os.getpid()
+ tmp_map = 'generaluse_map_tmp_%d' % os.getpid()
+
+ nb_lines=grass.vector_info_topo(input)['lines']
+
+ #Find best reference line and max distance between centerpoints of lines
+ segment_input=''
+ categories=grass.pipe_command('v.category', input=input, option='print',
+ quiet=True)
+ for category in categories.stdout:
+ segment_input += 'P ' + category.strip()
+ segment_input += ' ' + category.strip() + ' 50%\n'
+
+ grass.write_command('v.segment', input=input, output=tmp_centerpoints_map,
+ file='-', stdin=segment_input, quiet=True)
+
+ center_distances=grass.pipe_command('v.distance',
+ _from=tmp_centerpoints_map, to=tmp_centerpoints_map, upload='dist',
+ column='dist', flags='pa', quiet=True)
+
+ cats=[]
+ mean_dists=[]
+ count=0
+ distmax=0
+ for center in center_distances.stdout:
+ if count<2:
+ count+=1
+ continue
+ cat = center.strip().split('|')[0]
+ distsum=0
+ for x in center.strip().split('|')[1:]:
+ distsum+=float(x)
+ mean_dist = distsum/len(center.strip().split('|')[1:])
+ cats.append(cat)
+ mean_dists.append(mean_dist)
+
+ if transversals and not search_range:
+ search_range=sum(mean_dists)/len(mean_dists)
+ grass.message(_("Calculated search range: %.5f." % search_range))
+
+ if not refline_cat:
+ refline_cat = sorted(zip(cats, mean_dists), key = lambda tup: tup[1])[0][0]
+ grass.message(_("Line with category number %s was chosen as reference line." % refline_cat))
+
+ #Use transversals algorithm
+ if transversals:
+
+ #Break any intersections in the original lines so that they do not interfere
+ #further on
+ grass.run_command('v.clean', input=input, output=tmp_cleaned_map,
+ tool='break', quiet=True)
+
+ xmean=[]
+ ymean=[]
+ xmedian=[]
+ ymedian=[]
+ step=100.0/nb_vertices
+
+ os.environ['GRASS_VERBOSE']='-1'
+
+ for vertice in range(0,nb_vertices+1):
+ #v.segment sometimes cannot find points when using 0% or 100% offset
+ length_offset=step*vertice
+ if length_offset<0.00001:
+ length_offset=0.00001
+ if length_offset>99.99999:
+ length_offset=99.9999
+ #Create endpoints of transversal
+ segment_input = 'P 1 %s %.5f%% %f\n' % (refline_cat, length_offset,
+ search_range)
+ segment_input += 'P 2 %s %.5f%% %f\n' % (refline_cat, length_offset,
+ -search_range)
+ grass.write_command('v.segment', input=input, output=tmp_points_map,
+ stdin=segment_input, overwrite=True)
+
+ #Create transversal
+ grass.write_command('v.net', points=tmp_points_map, output=tmp_line_map,
+ operation='arcs', file='-', stdin='99999 1 2', overwrite=True)
+
+ #Patch transversal onto cleaned input lines
+ maps = tmp_cleaned_map + ',' + tmp_line_map
+ grass.run_command('v.patch', input=maps, out=tmp_map,
+ overwrite=True)
+
+ #Find intersections
+ grass.run_command('v.clean', input=tmp_map, out=tmp_line_map,
+ tool='break', error=tmp_points_map, overwrite=True)
+
+ #Add categories to intersection points
+ grass.run_command('v.category', input=tmp_points_map, out=tmp_map,
+ op='add', overwrite=True)
+
+ #Get coordinates of points
+ coords=grass.pipe_command('v.to.db', map=tmp_map,
+ op='coor', flags='p')
+
+ count=0
+ x=[]
+ y=[]
+ for coord in coords.stdout:
+ x.append(float(coord.strip().split('|')[1]))
+ y.append(float(coord.strip().split('|')[2]))
+
+ #Calculate mean and median for this transversal
+ if len(x) > 0:
+ xmean.append(sum(x)/len(x))
+ ymean.append(sum(y)/len(y))
+
+ x.sort()
+ y.sort()
+
+ xmedian.append((x[(len(x)-1)/2]+x[(len(x))/2])/2)
+ ymedian.append((y[(len(y)-1)/2]+y[(len(y))/2])/2)
+
+
+ del os.environ['GRASS_VERBOSE']
+
+
+
+ #Use closest point algorithm
+ else:
+
+ #Get reference line calculate its length
+ grass.run_command('v.extract', input=input, output=tmp_line_map,
+ cats=refline_cat, quiet=True)
+
+ os.environ['GRASS_VERBOSE']='0'
+ lpipe=grass.pipe_command('v.to.db', map=tmp_line_map, op='length',
+ flags='p')
+ del os.environ['GRASS_VERBOSE']
+
+ for l in lpipe.stdout:
+ linelength=float(l.strip().split('|')[1])
+
+ step=linelength/nb_vertices
+
+ #Create reference points for vertice calculation
+ grass.run_command('v.to.points', input=tmp_line_map,
+ output=tmp_points_map, dmax=step, quiet=True)
+
+ nb_points=grass.vector_info_topo(tmp_points_map)['points']
+
+ cat=[]
+ x=[]
+ y=[]
+
+ #Get coordinates of closest points on all input lines
+ if search_range:
+ points=grass.pipe_command('v.distance', _from=tmp_points_map,
+ from_layer=2, to=input, upload='to_x,to_y',
+ col='x,y', dmax=search_range, flags='pa', quiet=True)
+ else:
+ points=grass.pipe_command('v.distance', _from=tmp_points_map,
+ from_layer=2, to=input, upload='to_x,to_y',
+ col='x,y', flags='pa', quiet=True)
+
+ firstline=True
+ for point in points.stdout:
+ if firstline:
+ firstline=False
+ continue
+ cat.append((int(point.strip().split('|')[0])))
+ x.append(float(point.strip().split('|')[2]))
+ y.append(float(point.strip().split('|')[3]))
+
+ #Calculate mean coordinates
+ xsum=[0]*nb_points
+ ysum=[0]*nb_points
+ linecount=[0]*nb_points
+
+ for i in range(len(cat)):
+ index=cat[i]-1
+ linecount[index]+=1
+ xsum[index]=xsum[index]+x[i]
+ ysum[index]=ysum[index]+y[i]
+
+ xmean=[0]*nb_points
+ ymean=[0]*nb_points
+
+ for c in range(0,nb_points):
+ xmean[c]=xsum[c]/linecount[c]
+ ymean[c]=ysum[c]/linecount[c]
+
+ #Calculate the median
+
+ xmedian=[0]*nb_points
+ ymedian=[0]*nb_points
+
+ for c in range(0,nb_points):
+ xtemp=[]
+ ytemp=[]
+ for i in range(len(cat)):
+ if cat[i] == c+1:
+ xtemp.append(x[i])
+ ytemp.append(y[i])
+ xtemp.sort()
+ ytemp.sort()
+ xmedian[c]=(xtemp[(len(xtemp)-1)/2]+xtemp[(len(xtemp))/2])/2
+ ymedian[c]=(ytemp[(len(ytemp)-1)/2]+ytemp[(len(ytemp))/2])/2
+
+
+ #Create new line and write to file
+ if median and nb_lines > 2 :
+ line = geo.Line(zip(xmedian,ymedian))
+ else:
+ if median and nb_lines <= 2:
+ grass.message(_("More than 2 lines necesary for median, using mean."))
+ line = geo.Line(zip(xmean,ymean))
+
+ new = VectorTopo(output)
+ new.open('w')
+
+ new.write(line)
+ new.close()
+
+if __name__=="__main__":
+ options,flags=grass.parser()
+ atexit.register(cleanup)
+ main()
+
Property changes on: grass-addons/grass7/vector/v.centerline/v.centerline.py
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list