[GRASS-SVN] r49931 - in grass-addons/grass6/vector: . v.transects
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Dec 26 21:04:58 EST 2011
Author: helena
Date: 2011-12-26 18:04:58 -0800 (Mon, 26 Dec 2011)
New Revision: 49931
Added:
grass-addons/grass6/vector/v.transects/
grass-addons/grass6/vector/v.transects/v.transects.html
grass-addons/grass6/vector/v.transects/v.transects.py
Log:
script for generating transects perpendicular to line features
Added: grass-addons/grass6/vector/v.transects/v.transects.html
===================================================================
--- grass-addons/grass6/vector/v.transects/v.transects.html (rev 0)
+++ grass-addons/grass6/vector/v.transects/v.transects.html 2011-12-27 02:04:58 UTC (rev 49931)
@@ -0,0 +1,115 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+<head>
+<title>GRASS GIS manual: v.transects.py</title>
+<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+<link rel="stylesheet" href="grassdocs.css" type="text/css">
+</head>
+<body bgcolor="white">
+
+<img src="grass_logo.png" alt="GRASS logo"><hr align=center size=6 noshade>
+
+<h2>NAME</h2>
+<em><b>v.transects.py</b></em> - Python script to create vector line or area transects along a line
+<h2>KEYWORDS</h2>
+vector, transects, shoreline
+
+<h2>SYNOPSIS</h2>
+<b>v.transects.py</b> <b>input</b>=<em>name</em> <b>output</b>=<em>name</em> <b>transect_spacing</b>=<em>float</em> [<b>dleft</b>=<em>float</em>] [<b>dright</b>=<em>float</em>] [<b>type</b>="line" | "area"]
+
+
+
+<h3>Parameters:</h3>
+<DL>
+<DT><b>input</b>=<em>name</em></DT>
+<DD>Name for input line vector map</DD>
+
+<DT><b>output</b>=<em>name</em></DT>
+<DD>Name for output transects vector map</DD>
+
+<DT><b>transect_spacing</b>=<em>float</em></DT>
+<DD>Distance between transects</DD>
+
+<DT><b>dleft</b>=<em>float</em></DT>
+<DD>Distance transects extend on left side of line</DD>
+<DD>Default: <em>transect_spacing</em></DD>
+
+<DT><b>dright</b>=<em>float</em></DT>
+<DD>Distance transects extend on right side of line</DD>
+<DD>Default: <em>transect_spacing</em></DD>
+
+<DT><b>type</b>="line" | "area"</DT>
+<DD>Specifies the geometry of the transect</DD>
+<DD>Default: "line"</DD>
+
+</DL>
+
+
+<H2>DESCRIPTION</H2><EM>v.transects</EM> creates equally spaced geometries along
+input lines. The geometries can be lines or quadrilateral areas. Lines
+and areas are generated to be perpendicular to the input line.
+
+
+<H2>NOTES</H2>Input vector lines that are shorter than <B>transect_spacing</B> are ignored.
+<p>
+<B>transect_spacing</B>, <B>dleft</B>, and <B>dright</B> are interpreted to be in horizontal map units (e.g., degrees in the WGS84 projection).
+<p>
+<EM><B>v.transects</B></EM> may fail for a network of lines in Windows.
+
+<H2>EXAMPLES</H2>In these examples, the
+<a href="http://courses.ncsu.edu/mea582/common/media/01/NagsHead_series.zip"><font color="0000FF">Nags Head (19MB)</font></a>
+mapset is used to generate a shoreline and shore-perpendicular geometries. To use the mapset, unpack it into the
+<a href="http://courses.ncsu.edu/mea582/common/media/01/nc_spm_08.zip"><font color="0000FF">nc_spm_08 (150MB)</font></a>
+Location or the
+<a href="http://courses.ncsu.edu/mea792/common/media/gisdemo.zip"><font color="0000FF">gisdemo (47MB)</font></a>
+Location.
+
+
+<h3>Example 1) - Generate line transects along shoreline</h3>
+
+<P>Generate 20 cross-shore transects along 2008 shoreline (1m contour)<BR>
+<DIV class=code><PRE>g.region rast=NH_2008_1m
+r.contour input=NH_2008_1m output=NH_2008_1m level=1 cut=100
+v.report map=NH_2008_1m option=length
+# cat|level|length
+# 1|1|1037.86684790028
+# 1038m / 20transects = 52m per transect (value for transect_spacing)
+v.transects input=NH_2008_1m output=NH_2008_transects transect_spacing=52
+v.info NH_2008_transects
+</PRE></DIV>
+
+<h3>Example 2) - Generate line transects specifying the left and right length</h3>
+
+<P>Generate longer, more parallel transects by specifying dleft and dright and
+smoothing the input line<BR>
+<DIV class=code><PRE>g.region rast=NH_2008_1m
+r.contour input=NH_2008_1m output=NH_2008_1m level=1 cut=100
+v.generalize input=NH_2008_1m output=NH_2008_1m_smoothed \
+ method=sliding_averaging look_ahead=201
+v.transects input=NH_2008_1m_smoothed \
+ output=NH_2008_transects_long_smoothed transect_spacing=52 \
+ dleft=20 dright=300
+</PRE></DIV>
+
+<h3>Example 3) - Generate area transects along shoreline</h3>
+
+<P>Generate longer, more parallel transects by specifying dleft and dright and
+smoothing the input line<BR>
+<DIV class=code><PRE>g.region rast=NH_2008_1m
+r.contour input=NH_2008_1m output=NH_2008_1m level=1 cut=100
+v.transects input=NH_2008_1m output=NH_2008_areas \
+ transect_spacing=52 dleft=20 dright=300 type=area
+v.db.addtable NH_2008_areas
+v.db.addcolumn map=NH_2008_areas columns='demStats DOUBLE PRECISION'
+v.rast.stats vector=NH_2008_areas raster=NH_2008_1m column_prefix=NH2008
+v.db.select NH_2008_areas
+</PRE></DIV>
+
+
+<H2>SEE ALSO</H2><EM><A
+href="http://www4.ncsu.edu/~ejhardi2/v.generalize.html">v.generalize</A>, <A
+href="http://www4.ncsu.edu/~ejhardi2/r.transect.html">r.transect</A> </EM>
+
+
+<H2>AUTHOR</H2>Eric Hardin, Helena Mitasova, Updates by John Lloyd
+<P><I>Last changed: $Date$</I> </P></BODY></HTML>
Added: grass-addons/grass6/vector/v.transects/v.transects.py
===================================================================
--- grass-addons/grass6/vector/v.transects/v.transects.py (rev 0)
+++ grass-addons/grass6/vector/v.transects/v.transects.py 2011-12-27 02:04:58 UTC (rev 49931)
@@ -0,0 +1,251 @@
+#!/usr/bin/env python
+#
+#############################################################################
+#
+# MODULE: v.transects.py
+# AUTHOR(S): Eric Hardin
+# PURPOSE: Creates lines or quadrilateral areas perpendicular to a polyline
+# COPYRIGHT: (C) 2011
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+# UPDATES: John Lloyd November 2011
+#
+#############################################################################
+#%Module
+#% description: Transect Generator
+#%End
+#%option
+#% key: input
+#% type: string
+#% description: Name of input vector map
+#% required : yes
+#%end
+#%option
+#% key: output
+#% type: string
+#% description: Name of output vector map
+#% required : yes
+#%end
+#%option
+#% key: transect_spacing
+#% type: double
+#% description: Transect spacing
+#% required : yes
+#%end
+#%option
+#% key: dleft
+#% type: double
+#% description: distance transect extends to the left of input line.
+#% required : no
+#%end
+#%option
+#% key: dright
+#% type: double
+#% description: distance transect extends to the right of input line.
+#% required : no
+#%end
+#%option
+#% key: type
+#% type: string
+#% description: Feature type
+#% required : no
+#% options: line,area
+#% answer : line
+#%end
+
+from subprocess import Popen, PIPE, STDOUT
+from numpy import array
+from math import sqrt
+import grass.script as grass
+import tempfile
+####################################
+# load vector lines into python list
+# returns v
+# len(v) = number of lines in vector map
+# len(v[i]) = number of vertices in ith line
+# v[i][j] = [ xij, yij ] ,i.e., jth vertex in ith line
+def loadVector( vector ):
+ expVecCmmd = 'v.out.ascii format=standard input='+vector
+# JL p = Popen(expVecCmmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=True)
+ p = Popen(expVecCmmd, shell=True, stdin=PIPE, stdout=PIPE, stderr=STDOUT, close_fds=False)
+ vectorAscii = p.stdout.read().strip('\n').split('\n')
+ l = 0
+ while 'ORGANIZATION' not in vectorAscii[l]:
+ l += 1
+ while ':' in vectorAscii[l]:
+ l += 1
+ v = []
+ while l < len(vectorAscii):
+ line = vectorAscii[l].split()
+ if line[0] in ['L','B','A']:
+ skip = len(line)-2
+ vertices = int(line[1])
+ l += 1
+ v.append([])
+ for i in range(vertices):
+ v[-1].append( map(float,vectorAscii[l].split()[:2]) )
+ l += 1
+ l += skip
+ elif line[0] in ['P','C','F','K']:
+ skip = len(line)-2
+ vertices = int(line[1])
+ l += 1
+ for i in range(vertices):
+ l += 1
+ l += skip
+ else:
+ grass.fatal(_("Problem with line: <%s>") % vectorAscii[l])
+ if len(v) < 1:
+ grass.fatal(_("Zero lines found in vector map <%s>") % vector)
+ return v
+
+####################################
+# get transects from input vector
+def getTransects( vector, dleft, dright, transect_spacing ):
+ transect_locs = [] # holds locations where transects should intersect input vector lines
+ for line in vector:
+ transect_locs.append([line[0]])
+ # i starts at the beginning of the line.
+ # j walks along the line until j and i are separated by a distance of transect_spacing.
+ # then, a transect is placed at j, i is moved to j, and this is iterated until the end of the line is reached
+ i = 0
+ j = i+1
+ while j < len(line):
+ d = dist( line[i],line[j] )
+ if d > transect_spacing:
+ r = (transect_spacing - dist(line[i],line[j-1]) )/( dist(line[i],line[j]) - dist(line[i],line[j-1]) )
+ transect_locs[-1].append( (r*array( line[j] ) + (1-r)*array( line[j-1] )).tolist() )
+ i = j-1
+ line[i] = transect_locs[-1][-1]
+ else:
+ j += 1
+ transect_ends = []
+ for transect in transect_locs:
+ if len(transect) < 2: # if a line in input vec was shorter than transect_spacing
+ continue # then don't put a transect on it
+ transect_ends.append([])
+ transect = array( transect )
+ v = NR( transect[0], transect[1] ) # vector pointing parallel to transect
+ transect_ends[-1].append( [ transect[0]+dleft*v, transect[0]-dright*v ] )
+ for i in range(1,len( transect )-1,1):
+ v = NR( transect[i-1], transect[i+1] )
+ transect_ends[-1].append( [ transect[i]+dleft*v, transect[i]-dright*v ] )
+ v = NR( transect[-2], transect[-1] )
+ transect_ends[-1].append( [ transect[-1]+dleft*v, transect[-1]-dright*v ] )
+ return transect_ends
+
+####################################
+# calculate scalar distance between two points
+def dist( ip, fp ):
+ return sqrt( (ip[0]-fp[0])**2 + (ip[1]-fp[1])**2 )
+
+####################################
+# take a vector, normalize and rotate it 90 degrees
+def NR( ip, fp ):
+ x = fp[0] - ip[0]
+ y = fp[1] - ip[1]
+ r = sqrt( x**2 + y**2 )
+ return array([ -y/r, x/r ])
+
+####################################
+# write transects
+def writeTransects( transects, output ):
+ transects_str = ''
+ for transect in transects:
+ transects_str += '\n'.join( [ 'L 2\n'+' '.join(map(str,end_points[0]))+'\n'+' '.join(map(str,end_points[1]))+'\n' for end_points in transect ] )
+ # JL Rewrote Temporary File Logic for Windows
+ _, temp_path = tempfile.mkstemp()
+ a = open(temp_path, 'w')
+ a.write( transects_str )
+ a.seek(0)
+ a.close()
+ grass.run_command('v.in.ascii', flags='n', input=temp_path, output=output, format='standard')
+
+
+####################################
+# writes areas
+def writeQuads( transects, output ):
+ quad_str = ''
+ cnt = 1
+ for line in transects:
+ for tran in range( len(line)-1 ):
+ pt1 = ' '.join( map(str,line[tran][0]) )
+ pt2 = ' '.join( map(str,line[tran][1]) )
+ pt3 = ' '.join( map(str,line[tran+1][1]) )
+ pt4 = ' '.join( map(str,line[tran+1][0]) )
+ pt5 = pt1
+ # centroid is the average of the four corners
+ c = ' '.join( map(str,[ 0.25*(line[tran][0][0]+line[tran][1][0]+line[tran+1][0][0]+line[tran+1][1][0]), 0.25*(line[tran][0][1]+line[tran][1][1]+line[tran+1][0][1]+line[tran+1][1][1]) ]) )
+ quad_str += 'B 5\n' + '\n'.join([pt1,pt2,pt3,pt4,pt5]) + '\n'
+ quad_str += 'C 1 1\n' + c + '\n1 ' + str(cnt) + '\n'
+ cnt += 1
+ # JL Rewrote Temporary File Logic for Windows
+ _, temp_path = tempfile.mkstemp()
+ a = open(temp_path, 'w')
+ a.write( quad_str )
+ a.seek(0)
+ a.close()
+ grass.run_command('v.in.ascii', flags='n', input=a.name, output=output, format='standard')
+ a.close()
+
+####################################
+# Main method
+def main():
+ vector = options['input']
+ output = options['output']
+ # JL Handling Invalid transect_spacing parameter
+ try:
+ transect_spacing = float(options['transect_spacing'])
+ except:
+ grass.fatal(_("Invalid transect_spacing value."))
+ if transect_spacing == 0.0:
+ grass.fatal(_("Zero invalid transect_spacing value."))
+ dleft = options['dleft']
+ dright = options['dright']
+ shape = options['type']
+ print dleft, transect_spacing
+ if not dleft:
+ dleft = transect_spacing
+ else:
+ # JL Handling Invalid dleft parameter
+ try:
+ dleft = float(dleft)
+ except:
+ grass.fatal(_("Invalid dleft value."))
+ if not dright:
+ dright = transect_spacing
+ else:
+ # JL Handling Invalid dright parameter
+ try:
+ dright = float(dright)
+ except:
+ grass.fatal(_("Invalid dright value."))
+ # check if input file does not exists
+ if not grass.find_file(vector, element='vector')['file']:
+ grass.fatal(_("<%s> does not exist.") % vector)
+ # check if output file exists
+ if grass.find_file(output, element='vector')['mapset'] == grass.gisenv()['MAPSET']:
+ if not grass.overwrite():
+ grass.fatal(_("output map <%s> exists") % output)
+
+ #JL Is the vector a line and does if have at least one feature?
+ info = grass.parse_command('v.info', flags = 't', map = vector)
+ if info['lines'] == '0':
+ grass.fatal(_("vector <%s> does not contain lines") % vector)
+
+ #################################
+ v = loadVector( vector )
+ transect_ends = getTransects( v, dleft, dright, transect_spacing )
+ if shape == 'line' or not shape:
+ writeTransects( transect_ends, output )
+ else:
+ writeQuads( transect_ends, output )
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
+
+
Property changes on: grass-addons/grass6/vector/v.transects/v.transects.py
___________________________________________________________________
Added: svn:executable
+ *
More information about the grass-commit
mailing list