Context Navigation

spatialInputUtil.py @ 9582

Last change on this file since 9582 was 9553, checked in by davies, 10 years ago
Bugfix related to bilinear raster lookup
Property svn:executable set to ``*
File size: 55.8 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):
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)
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
1174	return allBreakLines
1175
1176	#########################################
1177	def readListOfRiverWalls(rwfileList):
1178	"""
1179	Take a list with the names of riverwall input files
1180	[should be comma-separated x,y,elevation files]
1181
1182	The input file can OPTIONALLY have a first line defining some
1183	hydraulic parameters. A valid example is
1184
1185	Qfactor: 1.5, s1: 0.94
1186	200., 300., 0.5
1187	300., 400., 0.7
1188	....and so on..
1189
1190	Read their coordinates into a dict with their names, read for use by ANUGA
1191
1192	INPUT: rwfileList -- a list of riverwall filenames
1193	[e.g. from glob.glob('GIS/RiverWalls/*.csv')]
1194
1195	OUTPUT:
1196	dictionary with riverwalls [filenames are keys] AND
1197	dictionary with hydraulic parameters [filenames are keys]
1198	"""
1199	import numpy
1200
1201	allRiverWalls = {}
1202	allRiverWallPar = {}
1203	for rwfile in rwfileList:
1204	f = open(rwfile)
1205	firstLine = f.readline()
1206	f.close()
1207	# If the top line has any letters, assume it is a hydraulic
1208	# variables line
1209	hasLetters = any(c.isalpha() for c in firstLine)
1210	if(not hasLetters):
1211	allRiverWalls[rwfile] = \
1212	numpy.genfromtxt(rwfile,delimiter=",").tolist()
1213	allRiverWallPar[rwfile] = {}
1214	else:
1215	# Get the wall geometry
1216	allRiverWalls[rwfile] = \
1217	numpy.genfromtxt(rwfile,delimiter=",",skip_header=1).tolist()
1218	# Get the hydraulic par
1219	firstLine = firstLine.replace(' ', '') # No whitespace
1220	wallPar = firstLine.split(',')
1221	allRiverWallPar[rwfile] = {}
1222	for wp in wallPar:
1223	keyNameValue = wp.split(':')
1224	allRiverWallPar[rwfile][keyNameValue[0]] = \
1225	float(keyNameValue[1])
1226
1227	return allRiverWalls, allRiverWallPar
1228
1229	############################################################################
1230
1231	def combine_breakLines_and_riverWalls_for_mesh(breakLines, riverWalls):
1232	"""
1233	Combine breaklines and riverwalls for input in mesh generation,
1234	ensuring both have 2 coordinates only
1235	"""
1236	mesh_breakLines=riverWalls.values()+breakLines.values()
1237	for i in range(len(mesh_breakLines)):
1238	mesh_breakLines[i] =\
1239	[mesh_breakLines[i][j][0:2] for j in range(len(mesh_breakLines[i]))]
1240	return mesh_breakLines
1241
1242	############################################################################
1243	def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
1244	""" We sometimes have breaklines defining 2 edges of a channel,
1245	and wish to make a polygon from them
1246
1247	Can do this with the current function
1248
1249	INPUTS: pattern == character string containing pattern which
1250	is inside exactly 2 keys in breaklines
1251
1252	breakLinesIn = the breakLines dictionary
1253
1254	reverse2nd = True/False or None. Reverse the order of the
1255	2nd set of edges before making the polygon.
1256	If None, then we compute the distance between the
1257	first point on breakline1 and the first/last
1258	points on breakline2, and reverse2nd if the
1259	'distance from the first point' <
1260	'distance from the last point'
1261
1262	OUTPUT: Polygon made with the 2 breaklines
1263	"""
1264
1265	breakLines=copy.copy(breakLinesIn)
1266	bk=breakLines.keys()
1267
1268	# They can be pathnames from glob, and sometimes / and \\ get mixed up
1269	# Fix that here
1270	pattern_norm = os.path.normpath(pattern)
1271	bk_norm = [ os.path.normpath(bk_i) for bk_i in bk ]
1272
1273	matchers=matchInds(pattern_norm, bk_norm)
1274
1275	if(len(matchers)==0):
1276	msg = 'Cannot match ' + pattern + ' in breaklines file names'
1277	raise Exception, msg
1278
1279	if(len(matchers)!=2):
1280	print 'Need exactly 2 matches, but pattern matched these', bk[matchers]
1281
1282	# There are 2 matches
1283
1284	if(reverse2nd is None):
1285	# Automatically compute whether we should reverse the 2nd breakline
1286	bl1_0=breakLines[bk[matchers[0]]][0]
1287	bl2_0=breakLines[bk[matchers[1]]][0]
1288	bl2_N=breakLines[bk[matchers[1]]][-1]
1289	d0 = ((bl1_0[0]-bl2_0[0])2 + (bl1_0[1]-bl2_0[1])2)
1290	dN = ((bl1_0[0]-bl2_N[0])2 + (bl1_0[1]-bl2_N[1])2)
1291	if(d0<dN):
1292	reverse2nd = True
1293	else:
1294	reverse2nd = False
1295
1296	if(reverse2nd):
1297	breakLines[bk[matchers[1]]].reverse()
1298	polyOut=breakLines[bk[matchers[0]]] + breakLines[bk[matchers[1]]]
1299
1300	# If the breakLines values have > 2 entries (i.e. like for riverwalls),
1301	# remove the third
1302	if(len(polyOut[0])>2):
1303	polyOut=[polyOut[i][0:2] for i in range(len(polyOut))]
1304
1305	return polyOut
1306	###################
1307
1308	else: # gdal_available == False
1309	msg='Failed to import gdal/ogr modules --'\
1310	+ 'perhaps gdal python interface is not installed.'
1311
1312
1313
1314	def readShp_1PolyGeo(shapefile, dropLast=True):
1315	raise ImportError, msg
1316
1317	def readShp_1LineGeo(shapefile):
1318	raise ImportError, msg
1319
1320	def read_csv_optional_header(filename):
1321	raise ImportError, msg
1322
1323	def read_polygon(filename):
1324	raise ImportError, msg
1325
1326	def readShpPtsAndAttributes(shapefile):
1327	raise ImportError, msg
1328
1329	def read_points(filename):
1330	raise ImportError, msg
1331
1332	def ListPts2Wkb( ptsIn, geometry_type='line', appendFirstOnEnd=None):
1333	raise ImportError, msg
1334
1335	def Wkb2ListPts(wkb_geo, removeLast=False, drop_third_dimension=False):
1336	raise ImportError, msg
1337
1338	def compute_squared_distance_to_segment(pt, line):
1339	raise ImportError, msg
1340
1341	def find_nearest_segment(pt, segments):
1342	raise ImportError, msg
1343
1344	def shift_point_on_line(pt, lineIn, nearest_segment_index):
1345	raise ImportError, msg
1346
1347	def insert_intersection_point(intersectionPt, line_pts,
1348	point_movement_threshold,verbose=False):
1349	raise ImportError, msg
1350
1351	def check_polygon_is_small(intersection, buf, tol2=100.):
1352	raise ImportError, msg
1353
1354	def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0,
1355	buf=1.0e-06, tol2 = 100,
1356	verbose=True, nameFlag=''):
1357	raise ImportError, msg
1358
1359	def getRasterExtent(rasterFile, asPolygon=False):
1360	raise ImportError, msg
1361
1362	def rasterValuesAtPoints(xy, rasterFile, band=1):
1363	raise ImportError, msg
1364
1365
1366	def gridPointsInPolygon(polygon, approx_grid_spacing=[1.,1.], eps=1.0e-06):
1367	raise ImportError, msg
1368
1369
1370	def matchInds(pattern, stringList):
1371	raise ImportError, msg
1372
1373
1374	def add_intersections_to_domain_features(bounding_polygonIn,
1375	breakLinesIn={ }, riverWallsIn={ }, point_movement_threshold=0.,
1376	verbose=True):
1377	raise ImportError, msg
1378
1379
1380	def readRegionPtAreas(shapefile, convert_length_to_area=False):
1381	raise ImportError, msg
1382
1383	def readListOfBreakLines(shapefileList):
1384	raise ImportError, msg
1385
1386	def combine_breakLines_and_riverwalls_for_mesh(breakLines, riverWalls):
1387	raise ImportError, msg
1388
1389	def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
1390	raise ImportError, msg
1391	###################
1392
1393

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

Context Navigation

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

Download in other formats: