Context Navigation

source: trunk/anuga_core/anuga/utilities/spatialInputUtil.py @ 9679

Last change on this file since 9679 was 9655, checked in by davies, 10 years ago
Allowing polygon shapefiles used as breaklines to keep final point. Also printing more info on structure inlets
Property svn:executable set to ``*
File size: 56.0 KB

Line
1	"""
2
3	Routines to ease the import of spatial data to ANUGA
4
5	Key routines:
6	readShp_1PolyGeo, readShp_1LineGeo -- read SIMPLE shapefile geometries into ANUGA as a list of points.
7	The supported geometries are quite restricted, see help
8	read_polygon -- Reads either a polygon or line shapefile or a csv as a polygon
9	readShpPointsAndAttributes -- read a multi-point shapefile with its attributes into ANUGA
10
11	ListPts2Wkb -- (Probably for internal use) Convert a list of points to a
12	Wkb geometry, allowing us to use GDALs geometry tools
13	Wbk2ListPts -- reverse of ListPts2Wkb
14
15	addIntersectionPtsToLines -- (Probably for internal use, see add_intersections_to_domain_features)
16	Given 2 intersecting lines, add their intersection points to them, or
17	move existing points if they are within a given distance of the intersection.
18
19	add_intersections_to_domain_features -- Add intersections to bounding_polygon, breaklines, riverwalls
20	This sets them up for being passed to the mesh generator.
21	Points within a given distance of an intersection can be replaced
22	with the intersection point if desired
23
24	rasterValuesAtPoints -- (Probably for internal use, see quantityRasterFun in quantity_setting_functions)
25	Quite efficiently get raster cell values at points in any gdal-compatible raster
26	[gridPointsInPolygon could in future be linked with this to get raster values in a region,
27	if we develop a version of ANUGA with sub-grid topography]
28
29	readRegionPtAreas -- read a shapefile containing regionPtAreas -- xy coordinates + 1 attribute, which is
30	the mesh triangle side length (or area) limit. Can be passed as regionPtAreas in
31	the mesh generation stage.
32
33	readListOfBreakLines -- takes a list of shapefile names, each containing a single line geometry, and reads
34	it in as a dictionary of breaklines.
35
36	readListOfRiverWalls -- takes a list of csv names, each containing a xyz polylines defining riverwalls, with
37	an optional first line defining the non-default riverwall par, and returns a dictionary
38	of the riverwalls and the riverwall_Par
39
40	polygon_from_matching_breaklines -- given a pattern (string) matching exactly 2 breaklines in a directory,
41	convert them to a single polygon. This is useful to e.g. make
42	a polygon defining the extent of a channel that is defined from 2 breaklines
43
44
45	"""
46	import sys
47	import os
48	import os.path
49	import copy
50	import struct
51	import numpy
52	#from matplotlib import path
53	import anuga
54	from anuga.geometry.polygon import inside_polygon
55
56	try:
57	import gdal
58	import ogr
59	gdal_available = True
60	#import osr # Not needed here but important in general
61	except ImportError, err:
62	gdal_available = False
63
64
65	#####################################
66	if gdal_available:
67
68
69	def readShp_1PolyGeo(shapefile, dropLast=True):
70	"""
71	Read a "single polygon" shapefile into an ANUGA_polygon object
72	(a list of lists, each containing a polygon coordinate)
73
74	The attribute table is ignored, and there can only be a single geometry in the shapefile
75
76	INPUTS: shapefile -- full path to shapefile name
77	dropLast -- Logical. If so, cut the last point (which closes
78	the polygon). ANUGA uses this by default for e.g.
79	bounding polygons
80	"""
81
82	# Get the data
83	driver=ogr.GetDriverByName("ESRI Shapefile")
84	dataSrc=driver.Open(shapefile, 0)
85	#dataSrc=ogr.Open(shapefile)
86	layer=dataSrc.GetLayer()
87
88	# Check it is a polygon
89	layerType=ogr.GeometryTypeToName(layer.GetGeomType())
90	if not layerType=='Polygon':
91	msg= shapefile +' is not a polygon shapefile'
92	raise Exception, msg
93
94	# Need a single polygon
95	try:
96	assert(len(layer)==1)
97	except:
98	print shapefile
99
100	boundary_poly=[]
101	for feature in layer:
102	#geom=feature.GetGeometryReg()
103	boundary=feature.GetGeometryRef().Boundary().GetPoints()
104	boundary=[list(pts) for pts in boundary]
105	boundary_poly.extend(boundary)
106
107	if(dropLast):
108	# Return a list of points, without last point [= first point]
109	return boundary_poly[:-1]
110	else:
111	return boundary_poly
112
113	####################
114
115	def readShp_1LineGeo(shapefile):
116	"""
117	Read a "single-line" shapefile into a list of lists (each containing a point),
118	resembling an ANUGA polygon object
119
120	The attribute table is ignored, and there can only be a single geometry in the shapefile
121
122	INPUTS: shapefile -- full path to shapefile name
123	"""
124
125	driver=ogr.GetDriverByName("ESRI Shapefile")
126	dataSrc=driver.Open(shapefile, 0)
127	#dataSrc=ogr.Open(shapefile)
128	layer=dataSrc.GetLayer()
129
130	# Check it is a line
131	layerType=ogr.GeometryTypeToName(layer.GetGeomType())
132	if not layerType=='Line String':
133	msg= shapefile +' is not a line shapefile'
134	raise Exception, msg
135
136	# Need a single line
137	try:
138	assert len(layer)==1
139	except:
140	print shapefile
141
142	line_all=[]
143	for feature in layer:
144	line=feature.GetGeometryRef().GetPoints()
145	line=[list(pts) for pts in line]
146	line_all.extend(line)
147
148	return line_all
149
150	###########################################################################
151	def read_csv_optional_header(filename):
152	"""Read a csv file of numbers, which optionally has a single header
153	row containing alphabetical characters (which is ignored if it
154	exists)
155
156	INPUT:
157	@param filename -- name of appropriate file with ',' delimiter
158
159	OUTPUT:
160	A numpy array with the numeric data
161	"""
162
163	f=open(filename)
164	firstLine=f.readline()
165	f.close()
166	hasLetters=any(c.isalpha() for c in firstLine)
167	outPol = numpy.genfromtxt(
168	filename, delimiter=',',skip_header=int(hasLetters))
169
170	return outPol
171
172	###########################################################################
173	def read_polygon(filename, close_polygon_shapefiles=False):
174	"""
175
176	Read a shapefile (polygon or line) or a csv file as an anuga polygon
177
178	In the case of a csv file, we permit either 1 or 0 header rows, and
179	the file can have > 2 columns [but only the first 2 are used]
180
181	Try to automatically do the correct thing based on the filename
182
183	"""
184	# Check the file exists
185	msg= 'Could not read '+ filename
186	assert os.path.isfile(filename), msg
187
188	# Get the file extension type
189	fileNameNoExtension , fileExtension = os.path.splitext(filename)
190
191	if fileExtension == '.shp':
192	# Read as either a polygon or line shapefile
193	try:
194	outPol=readShp_1PolyGeo(filename, dropLast = (not close_polygon_shapefiles))
195	assert len(outPol)>1
196	except:
197	try:
198	outPol= readShp_1LineGeo(filename)
199	assert len(outPol)>1
200	except:
201	msg= 'Could not read '+ filename +\
202	' as either polygon or line shapefile'
203	raise Exception(msg)
204	else:
205	try:
206	outPol = read_csv_optional_header(filename)
207	# Only take the first 2 columns
208	outPol = outPol[:,0:2].tolist()
209	except:
210	msg = 'Failed reading polygon '+ filename +\
211	' with anuga.utilities.spatialInputUtils.read_polygon'
212	raise Exception(msg)
213
214	return outPol
215
216	####################
217
218	def readShpPtsAndAttributes(shapefile):
219	"""
220	Read a point shapefile with an attribute table into a list
221
222	INPUT: shapefile -- name of shapefile to read
223
224	OUTPUT: List with
225	[ list_of_points, list_of_attributes, names_of_attributes]
226	"""
227
228	driver=ogr.GetDriverByName("ESRI Shapefile")
229	dataSrc=driver.Open(shapefile, 0)
230	#dataSrc=ogr.Open(shapefile)
231	layer=dataSrc.GetLayer()
232
233	pts=[]
234	attribute=[]
235	for i, feature in enumerate(layer):
236	if(i==0):
237	attributeNames=feature.keys()
238	pt=feature.GetGeometryRef().GetPoints()
239	pt=[list(p) for p in pt]
240	pts.extend(pt)
241	att=[feature[i] for i in attributeNames]
242	attribute.extend(att)
243
244	return [pts, attribute, attributeNames]
245
246	########################################
247	def readShpPts(shapefile):
248	"""
249	Wrapper around readShpPtsAndAttributes
250	Only get the points
251	"""
252
253	out=readShpPtsAndAttributes(shapefile)[0]
254	return out
255
256	########################################
257	def read_points(filename):
258	"""
259	Reads x,y geometries from a point shapefile or csv file (anuga polygon type format),
260	and returns as a list of lists
261	"""
262	# Check the file exists
263	msg= 'Could not read '+ filename
264	assert os.path.isfile(filename), msg
265
266	# Get the file extension type
267	fileNameNoExtension , fileExtension = os.path.splitext(filename)
268
269	if fileExtension=='.shp':
270	try:
271	points = readShpPts(filename)
272	except:
273	msg = 'Could not read points from ' + filename
274	raise Exception, msg
275	else:
276	# Assume txt format
277	try:
278	#points=numpy.genfromtxt(filename,delimiter=',',skip_header=1)
279	#points=points[:,0:2].tolist()
280	points=anuga.read_polygon(filename)
281	except:
282	msg = 'Could not read points from ' + filename +\
283	'. Make sure it has a single header row, ' +\
284	'with comma separator, and the first 2 columns are x,y'
285	raise Exception, msg
286	return points
287
288	########################################
289	def ListPts2Wkb( ptsIn, geometry_type='line', appendFirstOnEnd=None):
290	"""
291	Convert list of points to a GDAl Wkb format
292	Can be either points, lines, or polygon
293
294	Purpose is that once data in in Wkb format, we can use GDAL's geometric operations
295	(e.g. to test for intersections, etc)
296
297	INPUT: ptsIn -- list of points in the format [[x0,y0], [x1, y1], ..., [xn, yn]]
298	Actually it is also ok if [x0,y0], .. is a tuple instead
299	geometry_type -- 'point' or 'line' or 'polygon'
300
301	appendFirstOnEnd -- logical. If true, add the first point to the
302	end. Probably wanted for polygons when they are unclosed in ANUGA
303
304	OUTPUT:
305	The points as a Wkb Geometry.
306	geometry_type='point' produces a MULTIPOINT Geometry
307	='line' produces a LINESTRING Geometry
308	='polygon' produces a POLYGON Geometry
309
310	FIXME: This might not sensibly-use the gdal geometry types (although it
311	passes our tests) -- consider revising
312	"""
313	# Avoid modifying ptsIn
314	pts=copy.copy(ptsIn)
315
316	if appendFirstOnEnd is None:
317	if(geometry_type=='polygon'):
318	appendFirstOnEnd=True
319	else:
320	appendFirstOnEnd=False
321
322	if appendFirstOnEnd:
323	pts.append(pts[0])
324
325	if(geometry_type=='point'):
326	data=ogr.Geometry(ogr.wkbMultiPoint)
327	elif(geometry_type=='line'):
328	data=ogr.Geometry(ogr.wkbLineString)
329	elif(geometry_type=='polygon'):
330	data=ogr.Geometry(ogr.wkbLinearRing)
331	else:
332	raise Exception, "Type must be either 'point' or 'line' or 'polygon'"
333
334	for i in range(len(pts)):
335	if(len(pts[i])==2):
336	if(geometry_type=='point'):
337	newPt=ogr.Geometry(ogr.wkbPoint)
338	newPt.AddPoint(pts[i][0], pts[i][1])
339	data.AddGeometryDirectly(newPt)
340	else:
341	data.AddPoint(pts[i][0], pts[i][1])
342	elif(len(pts[i])==3):
343	if(geometry_type=='point'):
344	newPt=ogr.Geometry(ogr.wkbPoint)
345	newPt.AddPoint(pts[i][0], pts[i][1], pts[i][2])
346	data.AddGeometryDirectly(newPt)
347	else:
348	data.AddPoint(pts[i][0], pts[i][1], pts[i][2])
349	else:
350	raise Exception, 'Points must be either 2 or 3 dimensional'
351
352	if(geometry_type=='polygon'):
353	poly = ogr.Geometry(ogr.wkbPolygon)
354	poly.AddGeometry(data)
355	data=poly
356
357	return(data)
358
359	############################################################################
360	def Wkb2ListPts(wkb_geo, removeLast=False, drop_third_dimension=False):
361	"""
362	Reverse of ListPts2Wkb
363	"""
364	if(wkb_geo.GetGeometryName()=='POLYGON'):
365	X=wkb_geo.GetBoundary()
366	new=[ list(X.GetPoints()[i]) for i in range(len(X.GetPoints())) ]
367	elif(wkb_geo.GetGeometryName()=='MULTIPOINT'):
368	new=[ [feature.GetX(), feature.GetY(),feature.GetZ()] for feature in wkb_geo]
369	elif(wkb_geo.GetGeometryName()=='LINESTRING'):
370	new=[ list(wkb_geo.GetPoints()[i]) for i in range(len(wkb_geo.GetPoints())) ]
371	else:
372	raise Exception, 'Geometry type not supported'
373
374	if(removeLast):
375	new=new[:-1]
376	if(drop_third_dimension):
377	new = [new[i][0:2] for i in range(len(new))]
378	return new
379
380	############################################################################
381	def compute_squared_distance_to_segment(pt, line):
382	"""
383	Compute the squared distance between a point and a [finite] line segment
384
385	INPUT: pt -- [x,y]
386	line -- [[x0,y0],[x1,y1]] -- 2 points defining a line segment
387
388	OUTPUT: The distance^2 of pt to the line segment
389
390	"""
391	p0 = line[0]
392	p1 = line[1]
393	#
394	# Get unit vector along segment
395	seg_unitVec_x = float(p1[0]-p0[0])
396	seg_unitVec_y = float(p1[1]-p0[1])
397	segLen = (seg_unitVec_x2+seg_unitVec_y2)**0.5
398	if(segLen == 0.):
399	raise Exception, 'Line has repeated points: Line %s Pt %s' % (str(line),str(pt))
400
401	seg_unitVec_x = seg_unitVec_x/segLen
402	seg_unitVec_y = seg_unitVec_y/segLen
403
404	# Get vector from pt to p0
405	pt_p0_vec_x = float(pt[0]-p0[0])
406	pt_p0_vec_y = float(pt[1]-p0[1])
407	pt_p0_vec_len_squared = (pt_p0_vec_x2 + pt_p0_vec_y2)
408
409	# Get dot product of above vector with unit vector
410	pt_dot_segUnitVec = (pt_p0_vec_x)seg_unitVec_x + (pt_p0_vec_y)seg_unitVec_y
411
412	if( (pt_dot_segUnitVec < segLen) and (pt_dot_segUnitVec > 0.)):
413	# The nearest point on the line is actually between p0 and p1, so we have a 'real' candidate
414	# Get distance^2
415	output = pt_p0_vec_len_squared - pt_dot_segUnitVec**2.
416	else:
417	# Distance is the min distance from p0 and p1.
418	output = min( pt_p0_vec_len_squared, (float(pt[0]-p1[0])2+float(pt[1]-p1[1])2))
419
420	if(output < -1.0e-06):
421	print 'Diagnostic numbers follow: '
422	print output
423	print pt_p0_vec_len_squared
424	print pt_dot_segUnitVec
425	print pt
426	print p1
427	print p0
428	raise Exception, 'round-off in compute_squared_distance_to_segment'
429	if(output < 0.):
430	output=0.
431	return output
432
433	############################################################################
434
435	def find_nearest_segment(pt, segments):
436	"""
437	Given a point and a line, find the line segment nearest to the line
438
439	NOTE: The answer can be ambiguous if one of the segment endpoints is
440	the nearest point. In that case, the behaviour is determined by the behaviour
441	of numpy.argmin. Won't be a problem for this application
442
443	INPUT: pt -- [x,y]
444	segments -- [[x0,y0], [x1,y1], ...]
445	A list of points, consecutive points are interpreted
446	as joined and so defining line segments
447
448	OUTPUT: The squared distance, and the index i of the segment
449	[x_i,y_i],[x_i+1,y_i+1] in segments which is closest to pt
450	"""
451	ll=len(segments)
452	if(ll<=1):
453	raise Exception, 'Segments must have length > 1 in find_nearest_segment'
454
455	ptDist_sq=numpy.zeros(ll-1) # Hold the squared distance from the point to the line segment
456	for i in range(len(segments)-1):
457	# Compute distance from segment
458	ptDist_sq[i]=compute_squared_distance_to_segment(pt, [segments[i],segments[i+1]])
459
460	return [ptDist_sq.min(), ptDist_sq.argmin()]
461
462	######################################################
463
464	def shift_point_on_line(pt, lineIn, nearest_segment_index):
465	"""
466	Support pt is a point, which is near to the 'nearest_segment_index'
467	segment of the line 'lineIn'
468
469	This routine moves the nearest end point of that segment on line to pt.
470
471	INPUTS: pt -- [x,y] point
472	lineIn -- [ [x0, y0], [x1, y1], ..., [xN,yN]]
473	nearest_segment_index = index where the distance of pt to
474	the line from [x_i,y_i] to [x_i+1,y_i+1] is minimum
475
476	OUTPUT: The new line
477	"""
478	# Avoid Changing line
479	line=copy.copy(lineIn)
480
481	# Replace the nearest point on L1 with the intersection point
482	p0 = line[nearest_segment_index]
483	p1 = line[nearest_segment_index+1]
484	d_p0 = ( (pt[0]-p0[0])2 + (pt[1]-p0[1])2)
485	d_p1 = ( (pt[0]-p1[0])2 + (pt[1]-p1[1])2)
486	changeP1=(d_p1<d_p0)
487	line[nearest_segment_index+changeP1][0] = pt[0]
488	line[nearest_segment_index+changeP1][1] = pt[1]
489
490	return line
491
492	#################################################################################
493	def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold, verbose=False):
494	"""
495	Add intersectionPt to line_pts, either by inserting it, or if a point on line_pts is
496	closer than point_movement_threshold, then by moving that point to the intersection point
497
498	INPUTS:
499	intersectionPt -- the intersection point [known to lie on line_pts]
500	line_pts -- ordered list of [x,y] points making a line.
501	point_movement_threshold -- used to decide to move or add intersectionPt
502
503	OUTPUT:
504	new version of lint_pts with the point added
505	"""
506	# Avoid pointer/copy issues
507	L1_pts = copy.copy(line_pts)
508	iP=copy.copy(intersectionPt)
509
510	# Adjust L1
511	tmp = find_nearest_segment(intersectionPt, L1_pts)
512	# Compute the distance from the end points of the segment to the
513	# intersection point. Based on this we decide to add or move the point
514	p0 = L1_pts[tmp[1]]
515	p1 = L1_pts[tmp[1]+1]
516	endPt_Dist_Sq=min( ( (p0[0]-iP[0])2 + (p0[1]-iP[1])2),
517	( (p1[0]-iP[0])2 + (p1[1]-iP[1])2))
518	#
519	if(endPt_Dist_Sq>point_movement_threshold**2):
520	# Insert the intersection point. We do this in a tricky way to
521	# account for the possibility of L1_pts having > 2 coordinates
522	# (riverWalls)
523	if verbose:
524	print ' Inserting new point'
525	dummyPt=copy.copy(L1_pts[tmp[1]])
526	L1_pts.insert(tmp[1]+1,dummyPt)
527	L1_pts[tmp[1]+1][0]=iP[0]
528	L1_pts[tmp[1]+1][1]=iP[1]
529	if(len(L1_pts[tmp[1]+1])==3):
530	# Treat 3rd coordinate
531	# Find distance of inserted point from neighbours, and
532	# Set 3rd coordinate as distance-weighted average of the others
533	d0=((L1_pts[tmp[1]][0]-L1_pts[tmp[1]+1][0])**2.+\
534	(L1_pts[tmp[1]][1]-L1_pts[tmp[1]+1][1])2.)0.5
535	d1=((L1_pts[tmp[1]+2][0]-L1_pts[tmp[1]+1][0])**2.+\
536	(L1_pts[tmp[1]+2][1]-L1_pts[tmp[1]+1][1])2.)0.5
537	L1_pts[tmp[1]+1][2] = (d0L1_pts[tmp[1]+2][2] + d1L1_pts[tmp[1]][2])/(d0+d1)
538
539	else:
540	if verbose:
541	print ' Shifting existing point'
542	# Move a point already on L1
543	L1_pts=shift_point_on_line(iP, L1_pts, tmp[1])
544
545	return L1_pts
546
547	#######################################################################################################
548	def check_polygon_is_small(intersection, buf, tol2=100.):
549	"""
550
551	Elsewhere in the code, we check whether lines intersect by buffering them
552	to polygons with a small buffer = buf, then getting the intersection.
553	[since intersection with polygons is supported by gdal, but apparently
554	not directly with lines].
555
556	The logic of our code only works with point intersections,
557	and it will fails if 2 lines overlap in a line.
558
559	We crudely check for this situation here, by ensuring that the intersection polygon is 'small'
560
561	WARNING: The gdal geometry routines may be a bit rough (?)
562	Intersections not very precise, etc (?)
563
564	INPUT: intersection -- intersection of the 2 lines [gdal geometry]
565	buf -- a length scale giving the size of the intersection extent that we expect for a point
566	tol2 -- Throw an error if the x or y extent is greater than buf*tol2. Seems this needs to be
567	large sometimes -- this might reflect the stated weaknesses of GDALs geometry routines?
568	OUTPUT: True/False
569	False should suggest that the intersection is not describing a point
570	"""
571
572	extent=intersection.GetEnvelope()
573	assert(len(extent)==4) # Make sure this assumption is valid
574	if( (abs(extent[0]-extent[1])>buftol2) or (abs(extent[2]-extent[3]) > buftol2)):
575	return False
576	else:
577	return True
578
579	#######################################################################################################
580
581	def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0, buf=1.0e-05, tol2 = 100,
582	verbose=True, nameFlag=''):
583	"""
584	Add intersection points to lines L1 and L2 if they intersect each other
585
586	This is useful e.g. so that intersections can be exact (important for
587	mesh generation in ANUGA)
588
589	It currently only supports point intersection of 2 lines.
590	Line intersections should fail gracefully
591
592	INPUTS: L1, L2 = Wkb LineString geometries
593
594	point_movement_threshold -- if the distance from the nearest
595	point on L1 or L2 to the intersection is < point_movement_threshold, then the
596	nearest point has its coordinates replaced with the intersection point. This is
597	to prevent points being too close to each other
598
599	buf = tolerence that is used to buffer lines to find
600	intersections. Probably doesn't need modification
601
602	tol2 = [see check_polygon_is_small] Probably doesn't need to change
603
604	nameFlag = This will be printed if intersection occurs.
605	Convenient way to display intersecting filenames
606
607	OUTPUTS: L1,L2 with intersection points added in the right places
608	"""
609
610	if(L1.Intersects(L2)):
611	# Get points on the lines
612	L1_pts=Wkb2ListPts(L1)
613	L2_pts=Wkb2ListPts(L2)
614
615	# Buffer lines by a small amount
616	L1_buf=L1.Buffer(buf)
617	L2_buf=L2.Buffer(buf)
618
619	# Get intersection point[s]
620	L1_L2_intersect=L1_buf.Intersection(L2_buf)
621	if(L1_L2_intersect.GetGeometryCount()==1):
622	if(not check_polygon_is_small(L1_L2_intersect, buf, tol2)):
623	msg = 'line intersection is not allowed. ' + \
624	'Envelope %s '% str(L1_L2_intersect.GetEnvelope())
625	raise Exception(msg)
626	# Seems to need special treatment with only 1 intersection point
627	intersectionPts=[L1_L2_intersect.Centroid().GetPoint()]
628	else:
629	intersectionPts=[]
630	for feature in L1_L2_intersect:
631	if(not check_polygon_is_small(feature, buf, tol2)):
632	print feature.GetEnvelope()
633	raise Exception, 'line intersection is not allowed'
634	intersectionPts.append(feature.Centroid().GetPoint())
635
636	if(verbose):
637	print nameFlag
638	print ' Treating intersections in ', len(intersectionPts) , ' locations'
639	print intersectionPts
640
641	# Insert the points into the line segments
642	for i in range(len(intersectionPts)):
643	L1_pts = insert_intersection_point(intersectionPts[i], L1_pts,
644	point_movement_threshold, verbose=verbose)
645	L2_pts = insert_intersection_point(intersectionPts[i], L2_pts,
646	point_movement_threshold, verbose=verbose)
647
648	# Convert to the input format
649	L1_pts=ListPts2Wkb(L1_pts,geometry_type='line')
650	L2_pts=ListPts2Wkb(L2_pts,geometry_type='line')
651
652	return [L1_pts, L2_pts]
653	else:
654	return [L1, L2]
655
656	###########################################################
657	def getRasterExtent(rasterFile, asPolygon=False):
658	"""
659	Sometimes we need to know the extent of a raster
660	i.e. the minimum x, maximum x, minimum y, and maximum y values
661
662	INPUT:
663	rasterFile -- a gdal compatible rasterfile
664	asPolygon -- if False, return [xmin,xmax,ymin,ymax].
665	If True, return [ [xmin,ymin],[xmax,ymin],[xmax,ymax],[xmin,ymax]]
666	OUTPUT
667	The extent as defined above
668
669	"""
670	raster = gdal.Open(rasterFile)
671	transform=raster.GetGeoTransform()
672	xOrigin = transform[0]
673	yOrigin = transform[3]
674	xPixels=raster.RasterXSize
675	yPixels=raster.RasterYSize
676
677	# Compute the other extreme corner
678	x2=xOrigin + xPixelstransform[1]+yPixelstransform[2]
679	y2=yOrigin + xPixelstransform[4]+yPixelstransform[5]
680
681	xmin=min(xOrigin,x2)
682	xmax=max(xOrigin,x2)
683
684	ymin=min(yOrigin,y2)
685	ymax=max(yOrigin,y2)
686
687	if(asPolygon):
688	return [ [xmin,ymin], [xmax,ymin], [xmax,ymax], [xmin,ymax]]
689	else:
690	return [xmin,xmax,ymin,ymax]
691
692
693	###########################################################
694	def rasterValuesAtPoints(
695	xy,
696	rasterFile,
697	band=1,
698	nodata_rel_tol = 1.0e-08,
699	interpolation = 'pixel'):
700	"""
701	Get raster values at point locations.
702	Can be used to e.g. set quantity values
703
704	INPUT:
705	@param xy = numpy array with point locations
706
707	@param rasterFile = Filename of the gdal-compatible raster
708
709	@param band = band of the raster to get
710
711	@param nodata_rel_tol = Values are treated as nodata if
712	( abs(elev - nodataval) < nodata_rel_tol*abs(nodataval) )
713	This allows for truncation errors in nodata values which seem
714	to be introduced by some file-type conversions
715
716	@param interpolation 'pixel' or 'bilinear' determines how the
717	raster cell values are used to set the point value
718
719	OUTPUT:
720	1d numpy array with raster values at xy
721
722	"""
723	# Raster info
724	raster = gdal.Open(rasterFile)
725	rasterBand=raster.GetRasterBand(band)
726	rasterBandType=gdal.GetDataTypeName(rasterBand.DataType)
727	nodataval = rasterBand.GetNoDataValue()
728
729	# Projection info
730	transform=raster.GetGeoTransform()
731	xOrigin = transform[0]
732	yOrigin = transform[3]
733	pixelWidth = transform[1]
734	pixelHeight = transform[5] # Negative
735
736	# Get coordinates in pixel values
737	px = (xy[:,0] - xOrigin) / pixelWidth
738	py = (xy[:,1] - yOrigin) / pixelHeight
739
740	# Hold elevation
741	elev = px*0.
742
743	# Get the right character for struct.unpack
744	if (rasterBandType == 'Int16'):
745	CtypeName='h'
746	elif (rasterBandType == 'Float32'):
747	CtypeName='f'
748	elif (rasterBandType == 'Float64'):
749	CtypeName='d'
750	elif (rasterBandType == 'Byte'):
751	CtypeName='B'
752	elif (rasterBandType == 'Int32'):
753	CtypeName='i'
754	else:
755	print 'unrecognized DataType:', rasterBandType
756	print 'You might need to edit this code to read the data type'
757	raise Exception, 'Stopping'
758
759	# Upper bounds for pixel values, so we can fail gracefully
760	xMax = raster.RasterXSize
761	yMax = raster.RasterYSize
762	if(px.max()<xMax and px.min()>=0 and py.max()<yMax and py.min()>=0):
763	pass
764	else:
765	msg = 'Trying to extract point values that exceed the raster extent'
766	raise Exception, msg
767
768	# Get values -- seems we have to loop, but it is efficient enough
769	for i in range(len(px)):
770
771	if(interpolation=='pixel'):
772	# Pixel coordinates
773	xC=int(numpy.floor(px[i]))
774	yC=int(numpy.floor(py[i]))
775
776	structval = rasterBand.ReadRaster(xC,yC,1,1,
777	buf_type=rasterBand.DataType)
778	elev[i] = struct.unpack(CtypeName, structval)[0]
779
780	elif(interpolation=='bilinear'):
781	# Pixel coordinates
782	xl = int(numpy.floor(px[i]))
783	yl = int(numpy.floor(py[i]))
784
785	# Find neighbours required for bilinear interpolation
786	# l = lower, u = upper
787	if(px[i] - xl > 0.5):
788	xu = min(xl + 1, xMax - 1)
789	else:
790	# Swap xl for xu
791	xu = xl + 0
792	xl = max(xu - 1, 0)
793
794	if(py[i] - yl > 0.5):
795	yu = min(yl + 1, yMax - 1)
796	else:
797	yu = yl + 0
798	yl = max(yu - 1, 0)
799
800	# Map x,y to unit square
801	if(xu > xl):
802	x = px[i] - (xl + 0.5)
803	else:
804	x = 0.
805
806	if(yu > yl):
807	y = py[i] - (yl + 0.5)
808	else:
809	y = 0.
810
811	if not ( (x>=0.) & (x<=1.)):
812	print 'x-values error: ', x, xl, xu, px[i], xMax
813	raise Exception('x out of bounds')
814
815	if not ( (y>=0.) & (y<=1.)):
816	print 'y-values error: ', y, yl, yu, py[i]
817	raise Exception('y out of bounds')
818
819	# Lower-left
820	structval = rasterBand.ReadRaster(xl,yl,1,1,
821	buf_type=rasterBand.DataType)
822	r00 = struct.unpack(CtypeName, structval)[0]
823	# Upper left
824	structval = rasterBand.ReadRaster(xl,yu,1,1,
825	buf_type=rasterBand.DataType)
826	r01 = struct.unpack(CtypeName, structval)[0]
827	# Lower-right
828	structval = rasterBand.ReadRaster(xu,yl,1,1,
829	buf_type=rasterBand.DataType)
830	r10 = struct.unpack(CtypeName, structval)[0]
831	# Upper right
832	structval = rasterBand.ReadRaster(xu,yu,1,1,
833	buf_type=rasterBand.DataType)
834	r11 = struct.unpack(CtypeName, structval)[0]
835
836	# Bilinear interpolation
837	elev[i] = r00(1.-x)(1.-y) + r01(1.-x)y +\
838	r10x(1.-y) + r11xy
839
840	# Deal with nodata. This needs to be in the loop
841	# Just check if any of the pixels are nodata
842	if nodataval is not None:
843	if numpy.isfinite(nodataval):
844	rij = numpy.array([r00, r01, r10, r11])
845	rel_tol = ( abs(rij - nodataval) < nodata_rel_tol*abs(nodataval) )
846	missing = (rel_tol).nonzero()[0]
847	if len(missing) > 0:
848	elev[i] = numpy.nan
849	else:
850	raise Exception('Unknown value of "interpolation"')
851
852	# Deal with nodata for pixel based interpolation [efficient treatment
853	# outside of loop]
854	if (interpolation == 'pixel'):
855	if nodataval is not None:
856	if numpy.isfinite(nodataval):
857	rel_tol = ( abs(elev - nodataval) < nodata_rel_tol*abs(nodataval) )
858	missing = (rel_tol).nonzero()[0]
859	if len(missing) > 0:
860	elev[missing] = numpy.nan
861
862	return elev
863
864
865	def gridPointsInPolygon(polygon, approx_grid_spacing=[1.,1.], eps=1.0e-06):
866	"""
867	Get a 'grid' of points inside a polygon. Could be used with rasterValuesAtPoints to
868	get a range of raster values inside a polygon
869
870	Approach: A 'trial-grid' of points is created which is 'almost'
871	covering the range of the polygon (xmin-xmax,ymin-ymax).
872
873	(Actually it is just inside this region, to avoid polygon-boundary issues, see below)
874
875	Then we find those points which are actually inside the polygon.
876
877	The x/y point spacing on the trial-grid will be close to
878	approx_grid_spacing, but we ensure there are at least 4x4 points on the trial grid.
879	Also, we reduce the spacing so that the min_x+R and max_x-R
880	are both similarly close to the polygon extents [see definition of R below]
881
882	INPUTS:
883	polygon -- the polygon in ANUGA format (list of lists of ordered xy points)
884
885	approx_grid_spacing -- the approximate x,y grid spacing
886
887	eps -- 'trial-grid' of points has x range from min_polygon_x+R to
888	max_polygon_x - R, where R = (max_polygon_x-min_polygon_x)*eps
889	( and similarly for y).
890
891	This makes it more likely that points are inside the
892	polygon, not on the boundaries. Points on the boundaries can confuse the
893	point-in-polygon routine
894
895	OUTPUTS: A n x 2 numpy array of points in the polygon
896	"""
897
898	# Get polygon extent
899	polygonArr=numpy.array(polygon)
900	poly_xmin=polygonArr[:,0].min()
901	poly_xmax=polygonArr[:,0].max()
902	poly_ymin=polygonArr[:,1].min()
903	poly_ymax=polygonArr[:,1].max()
904
905	# Make a 'grid' of points which covers the polygon
906	xGridCount=max( numpy.ceil( (poly_xmax-poly_xmin)/approx_grid_spacing[0]+1. ).astype(int), 4)
907	R=(poly_xmax-poly_xmin)*eps
908	Xvals=numpy.linspace(poly_xmin+R,poly_xmax-R, xGridCount)
909	yGridCount=max( numpy.ceil( (poly_ymax-poly_ymin)/approx_grid_spacing[1]+1. ).astype(int), 4)
910	R=(poly_ymax-poly_ymin)*eps
911	Yvals=numpy.linspace(poly_ymin+R,poly_ymax-R, yGridCount)
912
913	xGrid,yGrid=numpy.meshgrid(Xvals,Yvals)
914	Grid=numpy.vstack([xGrid.flatten(),yGrid.flatten()]).transpose()
915
916	keepers = inside_polygon(Grid, polygon)
917	if(len(keepers)==0):
918	raise Exception('No points extracted from polygon')
919	xyInside=Grid[keepers,:]
920
921	return(xyInside)
922
923	#########################################################################
924	# Function to search for pattern matches in a string (turns out to be useful)
925	def matchInds(pattern, stringList):
926	"""
927	Find indices in stringList which match pattern
928	"""
929	#matches=[ (pattern in stringList[i]) for i in range(len(stringList))]
930	matches=[]
931	for i in range(len(stringList)):
932	if pattern in stringList[i]:
933	matches.append(i)
934	return matches
935
936
937	###########################################################################
938	#
939	# Less generic utilities below
940	#
941	# These are more 'anuga-specific' than above, aiming to make nice interfaces
942	# in ANUGA scripts
943	#
944	############################################################################
945
946
947	def add_intersections_to_domain_features(
948	bounding_polygonIn,
949	breakLinesIn={ },
950	riverWallsIn={ },
951	point_movement_threshold=0.,
952	verbose=True):
953	"""
954	If bounding polygon / breaklines /riverwalls intersect with each
955	other, then add intersection points.
956
957	INPUTS:
958	bounding_polygonIn -- the bounding polygon in ANUGA format
959	breakLinesIn -- the breaklines dictionary
960	riverWallsIn -- the riverWalls dictionary
961	point_movement_threshold -- if points on lines
962	are < this distance from intersection points, then they are
963	replaced with the intersection point. This can prevent
964	overly close points from breaking the mesh generation
965
966	OUTPUT:
967	List with bounding_polygon,breakLines,riverwalls
968	"""
969
970	bounding_polygon = copy.copy(bounding_polygonIn)
971	breakLines = copy.copy(breakLinesIn)
972	riverWalls = copy.copy(riverWallsIn)
973
974	# Quick exit
975	if (breakLines == {}) and (riverWalls == {}):
976	return [bounding_polygon, breakLines, riverWalls]
977
978	# Clean intersections of breakLines with itself
979	if(verbose):
980	print 'Cleaning breakline intersections'
981	if(len(breakLines)>0):
982	kbl = breakLines.keys()
983	for i in range(len(kbl)):
984	n1 = kbl[i]
985	for j in range(len(kbl)):
986	if(i >= j):
987	continue
988	n2 = kbl[j]
989	# Convert breaklines to wkb format
990	bl1 = ListPts2Wkb(breakLines[n1],geometry_type='line')
991	bl2 = ListPts2Wkb(breakLines[n2],geometry_type='line')
992	# Add intersection points
993	bl1, bl2 = addIntersectionPtsToLines(
994	bl1, bl2,
995	point_movement_threshold=point_movement_threshold,
996	verbose=verbose, nameFlag=n1+' intersects '+ n2)
997	breakLines[n1] = Wkb2ListPts(bl1)
998	breakLines[n2] = Wkb2ListPts(bl2)
999
1000
1001	# Clean intersections of riverWalls with itself
1002	if(verbose):
1003	print 'Cleaning riverWall intersections'
1004	if(len(riverWalls)>0):
1005	krw=riverWalls.keys()
1006	for i in range(len(krw)):
1007	n1=krw[i]
1008	for j in range(len(krw)):
1009	if(i>=j):
1010	continue
1011	n2=krw[j]
1012	# Convert breaklines to wkb format
1013	rw1=ListPts2Wkb(riverWalls[n1],geometry_type='line')
1014	rw2=ListPts2Wkb(riverWalls[n2],geometry_type='line')
1015	# Add intersection points
1016	rw1, rw2 =addIntersectionPtsToLines(rw1, rw2,\
1017	point_movement_threshold=point_movement_threshold,\
1018	verbose=verbose, nameFlag=n1+' intersects '+ n2)
1019	riverWalls[n1]=Wkb2ListPts(rw1)
1020	riverWalls[n2]=Wkb2ListPts(rw2)
1021
1022	# Clean intersections of breaklines with riverwalls
1023	if(verbose):
1024	print 'Cleaning breakLine-riverWall intersections'
1025	if( (len(riverWalls)>0) and (len(breakLines)>0)):
1026	krw=riverWalls.keys()
1027	kbl=breakLines.keys()
1028	for i in range(len(krw)):
1029	n1=krw[i]
1030	for j in range(len(kbl)):
1031	n2=kbl[j]
1032	# Convert breaklines to wkb format
1033	rw1=ListPts2Wkb(riverWalls[n1],geometry_type='line')
1034	bw2=ListPts2Wkb(breakLines[n2],geometry_type='line')
1035	# Add intersection points
1036	rw1, bw2 =addIntersectionPtsToLines(rw1, bw2,\
1037	point_movement_threshold=point_movement_threshold,\
1038	verbose=verbose, nameFlag=n1+' intersects '+ n2)
1039	riverWalls[n1]=Wkb2ListPts(rw1)
1040	breakLines[n2]=Wkb2ListPts(bw2)
1041
1042
1043	# Clean intersections of bounding polygon and riverwalls
1044	if(verbose):
1045	print 'Cleaning bounding_poly-riverWall intersections'
1046	if( (len(riverWalls)>0)):
1047	krw=riverWalls.keys()
1048	for i in range(len(krw)):
1049	n1=krw[i]
1050	# Convert breaklines to wkb format
1051	rw1=ListPts2Wkb(riverWalls[n1],geometry_type='line')
1052	bp2=ListPts2Wkb(bounding_polygon,geometry_type='line', appendFirstOnEnd=True)
1053	# Add intersection points
1054	rw1, bp2 =addIntersectionPtsToLines(rw1, bp2,\
1055	point_movement_threshold=point_movement_threshold,\
1056	verbose=verbose, nameFlag='Bounding Pol intersects '+ n1)
1057	riverWalls[n1]=Wkb2ListPts(rw1)
1058	# Since the bounding polygon is a loop, the first/last points are the same
1059	# If one of these was moved, the other should be moved too. Since we
1060	# will drop the last bounding_polygon point, we only need to worry about the first
1061	bounding_polygon=Wkb2ListPts(bp2,removeLast=False)
1062	if(bounding_polygon[-1] is not bounding_polygon[0]):
1063	bounding_polygon[0]=bounding_polygon[-1]
1064	# Drop the last point
1065	bounding_polygon=bounding_polygon[:-1]
1066
1067	# Clean intersections of bounding polygon and breaklines
1068	if(verbose):
1069	print 'Cleaning bounding_poly-breaklines intersections'
1070	if( (len(breakLines)>0)):
1071	kbl=breakLines.keys()
1072	for i in range(len(kbl)):
1073	n1=kbl[i]
1074	# Convert breaklines to wkb format
1075	bl1=ListPts2Wkb(breakLines[n1],geometry_type='line')
1076	bp2=ListPts2Wkb(bounding_polygon,geometry_type='line', appendFirstOnEnd=True)
1077	# Add intersection points
1078	bl1, bp2 =addIntersectionPtsToLines(bl1, bp2,\
1079	point_movement_threshold=point_movement_threshold,
1080	verbose=verbose, nameFlag='Bounding Pol intersects '+n1)
1081	breakLines[n1]=Wkb2ListPts(bl1)
1082	# Since the bounding polygon is a loop, the first/last points are the same
1083	# If one of these was moved, the other should be moved too. Since we
1084	# will drop the last bp2 point, we only need to worry about the first
1085	bounding_polygon=Wkb2ListPts(bp2,removeLast=False)
1086	if(bounding_polygon[-1] is not bounding_polygon[0]):
1087	bounding_polygon[0]=bounding_polygon[-1]
1088	# Drop the last point
1089	bounding_polygon=bounding_polygon[:-1]
1090
1091	# Remove the extra 0.0 from bounding polygon [this cannot have 3 coordinates]
1092	bounding_polygon = [ bounding_polygon[i][0:2] for i in range(len(bounding_polygon))]
1093	# Same for breaklines [although might not matter]
1094	for n1 in breakLines.keys():
1095	breakLines[n1] = [breakLines[n1][i][0:2] for i in range(len(breakLines[n1]))]
1096
1097	# Check that all mesh-boundary points are inside the bounding polygon
1098	from anuga.geometry.polygon import outside_polygon
1099	for blCat in [riverWalls, breakLines]:
1100	for n1 in blCat.keys():
1101	l=len(blCat[n1])
1102	# Test every point -- means we can strip 3rd coordinate if needed
1103	for j in range(l):
1104	isOut=outside_polygon(blCat[n1][j][0:2], bounding_polygon)
1105	if(len(isOut)>0):
1106	msg='Breakline/riverwall point '+str(blCat[n1][j][0:2]) +' on '+ n1+\
1107	' is outside the bounding polygon.\n'+\
1108	'Check that it exceeds the bounding polygon'+\
1109	' by a distance < point_movement_threshold \n'+\
1110	' so it can be moved back onto the polygon'
1111	print 'Polygon\n '
1112	print bounding_polygon
1113	print 'Line \n'
1114	print blCat[n1]
1115	raise Exception, msg
1116
1117	return [bounding_polygon, breakLines, riverWalls]
1118
1119	###################################################################
1120
1121	def readRegionPtAreas(shapefile, convert_length_to_area=False):
1122	"""
1123	Read a point shapefile to define the ANUGA mesh resoutions.
1124
1125	MUST HAVE A SINGLE ATTRIBUTE REPRESENTING THE LENGTHS OF TRIANGLES IN
1126	REGIONS
1127
1128	INPUT: shapefile -- name of the input shapefile
1129	convert_length_to_area -- if True, res values are assumed to
1130	represent triangle side lengths, and are converted to areas with 0.5res0res0
1131	Note that this might not ensure that the max triangle side length really is res0, but
1132	it will be of similar magnitude
1133	If False, attribute values are assumed to represent triangle areas
1134
1135	OUTPUT: list of the form [ [x0,y0,res0], [x1, y1, res1], ...]
1136	"""
1137
1138	ptData=readShpPtsAndAttributes(shapefile)
1139
1140	# Must have only 1 attribute
1141	assert len(ptData[2])==1
1142
1143	numPts=len(ptData[0])
1144	outData=[]
1145	for i in range(numPts):
1146	if(convert_length_to_area):
1147	newDat=[ptData[0][i][0], ptData[0][i][1], 0.5float(ptData[1][i])*2]
1148	else:
1149	newDat=[ptData[0][i][0], ptData[0][i][1], float(ptData[1][i])]
1150	outData.append(newDat)
1151
1152	return outData
1153
1154	#########################################
1155	def readListOfBreakLines(fileList):
1156	"""
1157	Take a list with the names of shapefiles or anuga_polygon csv files
1158
1159	They are assumed to be '2D breaklines', so we just read their
1160	coordinates into a dict with their names
1161
1162	Read them in
1163
1164	INPUT: fileList -- a list of shapefile and/or anuga_polygon csv filenames
1165	[e.g. from glob.glob('GIS/Breaklines/*.shp')]
1166
1167	OUTPUT: dictionary with breaklines [filenames are keys]
1168	"""
1169
1170	allBreakLines = {}
1171	for shapefile in fileList:
1172	allBreakLines[shapefile] = read_polygon(shapefile,
1173	close_polygon_shapefiles=True)
1174
1175	return allBreakLines
1176
1177	#########################################
1178	def readListOfRiverWalls(rwfileList):
1179	"""
1180	Take a list with the names of riverwall input files
1181	[should be comma-separated x,y,elevation files]
1182
1183	The input file can OPTIONALLY have a first line defining some
1184	hydraulic parameters. A valid example is
1185
1186	Qfactor: 1.5, s1: 0.94
1187	200., 300., 0.5
1188	300., 400., 0.7
1189	....and so on..
1190
1191	Read their coordinates into a dict with their names, read for use by ANUGA
1192
1193	INPUT: rwfileList -- a list of riverwall filenames
1194	[e.g. from glob.glob('GIS/RiverWalls/*.csv')]
1195
1196	OUTPUT:
1197	dictionary with riverwalls [filenames are keys] AND
1198	dictionary with hydraulic parameters [filenames are keys]
1199	"""
1200	import numpy
1201
1202	allRiverWalls = {}
1203	allRiverWallPar = {}
1204	for rwfile in rwfileList:
1205	f = open(rwfile)
1206	firstLine = f.readline()
1207	f.close()
1208	# If the top line has any letters, assume it is a hydraulic
1209	# variables line
1210	hasLetters = any(c.isalpha() for c in firstLine)
1211	if(not hasLetters):
1212	allRiverWalls[rwfile] = \
1213	numpy.genfromtxt(rwfile,delimiter=",").tolist()
1214	allRiverWallPar[rwfile] = {}
1215	else:
1216	# Get the wall geometry
1217	allRiverWalls[rwfile] = \
1218	numpy.genfromtxt(rwfile,delimiter=",",skip_header=1).tolist()
1219	# Get the hydraulic par
1220	firstLine = firstLine.replace(' ', '') # No whitespace
1221	wallPar = firstLine.split(',')
1222	allRiverWallPar[rwfile] = {}
1223	for wp in wallPar:
1224	keyNameValue = wp.split(':')
1225	allRiverWallPar[rwfile][keyNameValue[0]] = \
1226	float(keyNameValue[1])
1227
1228	return allRiverWalls, allRiverWallPar
1229
1230	############################################################################
1231
1232	def combine_breakLines_and_riverWalls_for_mesh(breakLines, riverWalls):
1233	"""
1234	Combine breaklines and riverwalls for input in mesh generation,
1235	ensuring both have 2 coordinates only
1236	"""
1237	mesh_breakLines=riverWalls.values()+breakLines.values()
1238	for i in range(len(mesh_breakLines)):
1239	mesh_breakLines[i] =\
1240	[mesh_breakLines[i][j][0:2] for j in range(len(mesh_breakLines[i]))]
1241	return mesh_breakLines
1242
1243	############################################################################
1244	def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
1245	""" We sometimes have breaklines defining 2 edges of a channel,
1246	and wish to make a polygon from them
1247
1248	Can do this with the current function
1249
1250	INPUTS: pattern == character string containing pattern which
1251	is inside exactly 2 keys in breaklines
1252
1253	breakLinesIn = the breakLines dictionary
1254
1255	reverse2nd = True/False or None. Reverse the order of the
1256	2nd set of edges before making the polygon.
1257	If None, then we compute the distance between the
1258	first point on breakline1 and the first/last
1259	points on breakline2, and reverse2nd if the
1260	'distance from the first point' <
1261	'distance from the last point'
1262
1263	OUTPUT: Polygon made with the 2 breaklines
1264	"""
1265
1266	breakLines=copy.copy(breakLinesIn)
1267	bk=breakLines.keys()
1268
1269	# They can be pathnames from glob, and sometimes / and \\ get mixed up
1270	# Fix that here
1271	pattern_norm = os.path.normpath(pattern)
1272	bk_norm = [ os.path.normpath(bk_i) for bk_i in bk ]
1273
1274	matchers=matchInds(pattern_norm, bk_norm)
1275
1276	if(len(matchers)==0):
1277	msg = 'Cannot match ' + pattern + ' in breaklines file names'
1278	raise Exception, msg
1279
1280	if(len(matchers)!=2):
1281	print 'Need exactly 2 matches, but pattern matched these', bk[matchers]
1282
1283	# There are 2 matches
1284
1285	if(reverse2nd is None):
1286	# Automatically compute whether we should reverse the 2nd breakline
1287	bl1_0=breakLines[bk[matchers[0]]][0]
1288	bl2_0=breakLines[bk[matchers[1]]][0]
1289	bl2_N=breakLines[bk[matchers[1]]][-1]
1290	d0 = ((bl1_0[0]-bl2_0[0])2 + (bl1_0[1]-bl2_0[1])2)
1291	dN = ((bl1_0[0]-bl2_N[0])2 + (bl1_0[1]-bl2_N[1])2)
1292	if(d0<dN):
1293	reverse2nd = True
1294	else:
1295	reverse2nd = False
1296
1297	if(reverse2nd):
1298	breakLines[bk[matchers[1]]].reverse()
1299	polyOut=breakLines[bk[matchers[0]]] + breakLines[bk[matchers[1]]]
1300
1301	# If the breakLines values have > 2 entries (i.e. like for riverwalls),
1302	# remove the third
1303	if(len(polyOut[0])>2):
1304	polyOut=[polyOut[i][0:2] for i in range(len(polyOut))]
1305
1306	return polyOut
1307	###################
1308
1309	else: # gdal_available == False
1310	msg='Failed to import gdal/ogr modules --'\
1311	+ 'perhaps gdal python interface is not installed.'
1312
1313
1314
1315	def readShp_1PolyGeo(shapefile, dropLast=True):
1316	raise ImportError, msg
1317
1318	def readShp_1LineGeo(shapefile):
1319	raise ImportError, msg
1320
1321	def read_csv_optional_header(filename):
1322	raise ImportError, msg
1323
1324	def read_polygon(filename):
1325	raise ImportError, msg
1326
1327	def readShpPtsAndAttributes(shapefile):
1328	raise ImportError, msg
1329
1330	def read_points(filename):
1331	raise ImportError, msg
1332
1333	def ListPts2Wkb( ptsIn, geometry_type='line', appendFirstOnEnd=None):
1334	raise ImportError, msg
1335
1336	def Wkb2ListPts(wkb_geo, removeLast=False, drop_third_dimension=False):
1337	raise ImportError, msg
1338
1339	def compute_squared_distance_to_segment(pt, line):
1340	raise ImportError, msg
1341
1342	def find_nearest_segment(pt, segments):
1343	raise ImportError, msg
1344
1345	def shift_point_on_line(pt, lineIn, nearest_segment_index):
1346	raise ImportError, msg
1347
1348	def insert_intersection_point(intersectionPt, line_pts,
1349	point_movement_threshold,verbose=False):
1350	raise ImportError, msg
1351
1352	def check_polygon_is_small(intersection, buf, tol2=100.):
1353	raise ImportError, msg
1354
1355	def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0,
1356	buf=1.0e-06, tol2 = 100,
1357	verbose=True, nameFlag=''):
1358	raise ImportError, msg
1359
1360	def getRasterExtent(rasterFile, asPolygon=False):
1361	raise ImportError, msg
1362
1363	def rasterValuesAtPoints(xy, rasterFile, band=1):
1364	raise ImportError, msg
1365
1366
1367	def gridPointsInPolygon(polygon, approx_grid_spacing=[1.,1.], eps=1.0e-06):
1368	raise ImportError, msg
1369
1370
1371	def matchInds(pattern, stringList):
1372	raise ImportError, msg
1373
1374
1375	def add_intersections_to_domain_features(bounding_polygonIn,
1376	breakLinesIn={ }, riverWallsIn={ }, point_movement_threshold=0.,
1377	verbose=True):
1378	raise ImportError, msg
1379
1380
1381	def readRegionPtAreas(shapefile, convert_length_to_area=False):
1382	raise ImportError, msg
1383
1384	def readListOfBreakLines(shapefileList):
1385	raise ImportError, msg
1386
1387	def combine_breakLines_and_riverwalls_for_mesh(breakLines, riverWalls):
1388	raise ImportError, msg
1389
1390	def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
1391	raise ImportError, msg
1392	###################
1393
1394

Note: See TracBrowser for help on using the repository browser.

Download in other formats: