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

Last change on this file since 9242 was 9242, checked in by davies, 10 years ago

Improvements + tests for polygon intersections in spatialInputUtils
  • Property svn:executable set to *
File size: 43.4 KB
Line 
1"""
2
3Routines to ease the import of spatial data to ANUGA
4
5Key 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    readShpPointsAndAttributes -- read a multi-point shapefile with its attributes into ANUGA
9
10    ListPts2Wkb -- (Probably for internal use) Convert a list of points to a
11                    Wkb geometry, allowing us to use GDALs geometry tools
12    Wbk2ListPts -- reverse of ListPts2Wkb
13   
14    addIntersectionPtsToLines -- (Probably for internal use, see add_intersections_to_domain_features)
15                                 Given 2 intersecting lines, add their intersection points to them, or
16                                 move existing points if they are within a given distance of the intersection.
17
18    add_intersections_to_domain_features -- Add intersections to bounding_polygon, breaklines, riverwalls
19                                            This sets them up for being passed to the mesh generator.
20                                            Points within a given distance of an intersection can be replaced
21                                            with the intersection point if desired
22
23    rasterValuesAtPoints -- (Probably for internal use, see quantityRasterFun in quantity_setting_functions)
24                            Quite efficiently get raster cell values at points in any gdal-compatible raster
25                            [gridPointsInPolygon could in future be linked with this to get raster values in a region,
26                             if we develop a version of ANUGA with sub-grid topography]
27
28    readRegionPtAreas -- read a shapefile containin regionPtAreas -- xy coordinates + 1 attribute, which is
29                         the mesh triangle side length (or area) limit. Can be passed as regionPtAreas in
30                         the mesh generation stage.
31
32    readListOfBreakLines -- takes a list of shapefile names, each containing a single line geometry, and reads
33                            it in as a dictionary of breaklines.
34
35    polygon_from_matching_breaklines -- given a pattern (string) matching exactly 2 breaklines in a directory,
36                                        convert them to a single polygon. This is useful to e.g. make
37                                        a polygon defining the extent of a channel that is defined from 2 breaklines
38
39   
40"""
41import sys
42import os
43import copy
44import struct
45import numpy
46from matplotlib import path
47
48try:
49    import gdal
50    import ogr
51    gdal_available = True
52    #import osr # Not needed here but important in general
53except ImportError, err:
54    gdal_available = False
55
56
57##################################### 
58if gdal_available:
59
60
61    def readShp_1PolyGeo(shapefile, dropLast=True):
62        """
63            Read a "single polygon" shapefile into an ANUGA_polygon object
64            (a list of lists, each containing a polygon coordinate)
65   
66            The attribute table is ignored, and there can only be a single geometry in the shapefile
67           
68            INPUTS: shapefile -- full path to shapefile name
69                    dropLast -- Logical. If so, cut the last point (which closes
70                                the polygon). ANUGA uses this by default for e.g.
71                                bounding polygons
72        """
73   
74        # Get the data
75        driver=ogr.GetDriverByName("ESRI Shapefile")
76        dataSrc=driver.Open(shapefile, 0)
77        #dataSrc=ogr.Open(shapefile)
78        layer=dataSrc.GetLayer()
79       
80        # Need a single polygon
81        try:
82            assert(len(layer)==1)
83        except:
84            print shapefile
85       
86        boundary_poly=[]
87        for feature in layer:
88            #geom=feature.GetGeometryReg()
89            boundary=feature.GetGeometryRef().Boundary().GetPoints()
90            boundary=[list(pts) for pts in boundary]
91            boundary_poly.extend(boundary)
92   
93        if(dropLast):
94            # Return a list of points, without last point [= first point]
95            return boundary_poly[:-1]
96        else:
97            return boundary_poly
98   
99    ####################
100   
101    def readShp_1LineGeo(shapefile):
102        """
103            Read a "single-line" shapefile into a list of lists (each containing a point),
104                resembling an ANUGA polygon object
105   
106            The attribute table is ignored, and there can only be a single geometry in the shapefile
107           
108            INPUTS: shapefile -- full path to shapefile name
109        """
110   
111        driver=ogr.GetDriverByName("ESRI Shapefile")
112        dataSrc=driver.Open(shapefile, 0)
113        #dataSrc=ogr.Open(shapefile)
114        layer=dataSrc.GetLayer()
115       
116        # Need a single line
117        try:
118            assert(len(layer)==1)
119        except:
120            print shapefile
121       
122        line_all=[]
123        for feature in layer:
124            line=feature.GetGeometryRef().GetPoints()
125            line=[list(pts) for pts in line]
126            line_all.extend(line)
127   
128        return line_all
129           
130    ####################
131   
132    def readShpPtsAndAttributes(shapefile):
133        """
134            Read a point shapefile with an attribute table into a list
135   
136            INPUT: shapefile -- name of shapefile to read
137           
138            OUTPUT: List with [ list_of_points, list_of_attributes, names_of_attributes]
139        """
140       
141        driver=ogr.GetDriverByName("ESRI Shapefile")
142        dataSrc=driver.Open(shapefile, 0)
143        #dataSrc=ogr.Open(shapefile)
144        layer=dataSrc.GetLayer()
145   
146        pts=[]
147        attribute=[]
148        for i, feature in enumerate(layer):
149            if(i==0):
150                attributeNames=feature.keys()
151            pt=feature.GetGeometryRef().GetPoints()
152            pt=[list(p) for p in pt]
153            pts.extend(pt)
154            att=[feature[i] for i in attributeNames]
155            attribute.extend(att)
156   
157        return [pts, attribute, attributeNames]
158   
159    ########################################
160    def ListPts2Wkb( ptsIn, geometry_type='line', appendFirstOnEnd=None):
161        """
162            Convert list of points to a GDAl Wkb format
163            Can be either points, lines, or polygon
164   
165            Purpose is that once data in in Wkb format, we can use GDAL's geometric operations
166                (e.g. to test for intersections, etc)
167           
168            INPUT: ptsIn -- list of points in the format [[x0,y0], [x1, y1], ..., [xn, yn]]
169                          Actually it is also ok if [x0,y0], .. is a tuple instead
170                   geometry_type -- 'point' or 'line' or 'polygon'
171   
172                   appendFirstOnEnd -- logical. If true, add the first point to the
173                        end. Probably wanted for polygons when they are unclosed in ANUGA
174   
175            OUTPUT:
176                    The points as a Wkb Geometry.
177                    geometry_type='point' produces a MULTIPOINT Geometry
178                                 ='line' produces a LINESTRING Geometry
179                                 ='polygon' produces a POLYGON Geometry
180   
181            FIXME: This might not sensibly-use the gdal geometry types (although it
182                    passes our tests) -- consider revising
183        """
184        # Avoid modifying ptsIn
185        pts=copy.copy(ptsIn)
186   
187        if appendFirstOnEnd is None:
188            if(geometry_type=='polygon'):
189                appendFirstOnEnd=True
190            else:
191                appendFirstOnEnd=False
192   
193        if appendFirstOnEnd:
194            pts.append(pts[0])   
195       
196        if(geometry_type=='point'):
197            data=ogr.Geometry(ogr.wkbMultiPoint)
198        elif(geometry_type=='line'):
199            data=ogr.Geometry(ogr.wkbLineString)
200        elif(geometry_type=='polygon'):
201            data=ogr.Geometry(ogr.wkbLinearRing)
202        else:
203            raise Exception, "Type must be either 'point' or 'line' or 'polygon'"
204         
205        for i in range(len(pts)):
206            if(len(pts[i])==2):
207                if(geometry_type=='point'):
208                    newPt=ogr.Geometry(ogr.wkbPoint)
209                    newPt.AddPoint(pts[i][0], pts[i][1]) 
210                    data.AddGeometryDirectly(newPt)
211                else:
212                    data.AddPoint(pts[i][0], pts[i][1])
213            elif(len(pts[i])==3):
214                if(geometry_type=='point'):
215                    newPt=ogr.Geometry(ogr.wkbPoint)
216                    newPt.AddPoint(pts[i][0], pts[i][1], pts[i][2]) 
217                    data.AddGeometryDirectly(newPt)
218                else:
219                    data.AddPoint(pts[i][0], pts[i][1], pts[i][2])
220            else:
221                raise Exception, 'Points must be either 2 or 3 dimensional'
222   
223        if(geometry_type=='polygon'):   
224            poly = ogr.Geometry(ogr.wkbPolygon)
225            poly.AddGeometry(data)
226            data=poly
227       
228        return(data)
229   
230    ############################################################################
231    def Wkb2ListPts(wkb_geo, removeLast=False, drop_third_dimension=False):
232        """
233            Reverse of ListPts2Wkb
234        """
235        if(wkb_geo.GetGeometryName()=='POLYGON'):
236            X=wkb_geo.GetBoundary()
237            new=[ list(X.GetPoints()[i]) for i in range(len(X.GetPoints())) ]
238        elif(wkb_geo.GetGeometryName()=='MULTIPOINT'):
239            new=[ [feature.GetX(), feature.GetY(),feature.GetZ()] for feature in wkb_geo]
240        elif(wkb_geo.GetGeometryName()=='LINESTRING'):
241            new=[ list(wkb_geo.GetPoints()[i]) for i in range(len(wkb_geo.GetPoints())) ]
242        else:
243            raise Exception, 'Geometry type not supported'
244       
245        if(removeLast):
246            new=new[:-1]
247        if(drop_third_dimension):
248            new = [new[i][0:2] for i in range(len(new))]
249        return new
250   
251    ############################################################################
252    def compute_squared_distance_to_segment(pt, line):
253        """
254            Compute the squared distance between a point and a [finite] line segment
255   
256            INPUT: pt -- [x,y]
257                   line -- [[x0,y0],[x1,y1]] -- 2 points defining a line segment
258   
259            OUTPUT: The distance^2 of pt to the line segment
260   
261        """
262        p0 = line[0]
263        p1 = line[1]
264        #
265        # Get unit vector along segment
266        seg_unitVec_x=float(p1[0]-p0[0])
267        seg_unitVec_y=float(p1[1]-p0[1])
268        segLen=(seg_unitVec_x**2+seg_unitVec_y**2)**0.5
269        if(segLen==0.):
270            #print line
271            #print 'Pt'
272            #print pt
273            raise Exception, 'Line has repeated points: Line %s Pt %s' % (str(line),str(pt))
274        seg_unitVec_x=seg_unitVec_x/segLen
275        seg_unitVec_y=seg_unitVec_y/segLen
276        #
277        # Get vector from pt to p0
278        pt_p0_vec_x=float(pt[0]-p0[0])
279        pt_p0_vec_y=float(pt[1]-p0[1])
280        pt_p0_vec_len_squared=(pt_p0_vec_x**2 + pt_p0_vec_y**2)
281        # Get dot product of above vector with unit vector
282        pt_dot_segUnitVec=(pt_p0_vec_x)*seg_unitVec_x+(pt_p0_vec_y)*seg_unitVec_y
283        #
284        if( (pt_dot_segUnitVec<segLen) and (pt_dot_segUnitVec > 0.)):
285            # The nearest point on the line is actually between p0 and p1, so we have a 'real' candidate
286            # Get distance^2
287            output = pt_p0_vec_len_squared - pt_dot_segUnitVec**2.
288        else:
289            # Distance is the min distance from p0 and p1.
290            output = min( pt_p0_vec_len_squared,  (float(pt[0]-p1[0])**2+float(pt[1]-p1[1])**2))
291        if(output < -1.0e-06):
292            raise Exception, 'round-off in compute_squared_distance_to_segment'
293        if(output < 0.):
294            output=0.
295        return output
296   
297    ############################################################################
298   
299    def find_nearest_segment(pt, segments):
300        """
301            Given a point and a line, find the line segment nearest to the line
302   
303            NOTE: The answer can be ambiguous if one of the segment endpoints is
304                the nearest point. In that case, the behaviour is determined by the behaviour
305                of numpy.argmin. Won't be a problem for this application
306   
307            INPUT: pt -- [x,y]
308                   segments -- [[x0,y0], [x1,y1], ...]
309                             A list of points, consecutive points are interpreted
310                             as joined and so defining line segments
311           
312            OUTPUT: The squared distance, and the index i of the segment [x_i,y_i],[x_i+1,y_i+1] in segments which is closest to pt
313        """
314        ll=len(segments)
315        if(ll<=1):
316            raise Exception, 'Segments must have length > 1 in find_nearest_segment'
317       
318        ptDist_sq=numpy.zeros(ll-1) # Hold the squared distance from the point to the line segment
319        for i in range(len(segments)-1):
320            # Compute distance from segment
321            ptDist_sq[i]=compute_squared_distance_to_segment(pt, [segments[i],segments[i+1]])
322       
323        return [ptDist_sq.min(), ptDist_sq.argmin()]
324   
325    ######################################################
326   
327    def shift_point_on_line(pt, lineIn, nearest_segment_index):
328        """
329            Support pt is a point, which is near to the 'nearest_segment_index' segment of the line 'lineIn'
330   
331            This routine moves the nearest end point of that segment on line to pt.
332   
333            INPUTS: pt -- [x,y] point
334                    lineIn -- [ [x0, y0], [x1, y1], ..., [xN,yN]]
335                    nearest_segment_index = index where the distance of pt to the line from [x_i,y_i] to [x_i+1,y_i+1] is minimum
336   
337            OUTPUT: The new line
338        """
339        # Avoid Changing line
340        line=copy.copy(lineIn)
341   
342        # Replace the nearest point on L1 with the intersection point
343        p0 = line[nearest_segment_index]
344        p1 = line[nearest_segment_index+1]
345        d_p0 = ( (pt[0]-p0[0])**2 + (pt[1]-p0[1])**2)
346        d_p1 = ( (pt[0]-p1[0])**2 + (pt[1]-p1[1])**2)
347        changeP1=(d_p1<d_p0)
348        line[nearest_segment_index+changeP1][0] = pt[0]
349        line[nearest_segment_index+changeP1][1] = pt[1]
350   
351        return line
352   
353    #################################################################################
354    def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold, verbose=False):
355        """
356            Add intersectionPt to line_pts, either by inserting it, or if a point on line_pts is
357            closer than point_movement_threshold, then by moving that point to the intersection point
358   
359            INPUTS:
360                    intersectionPt -- the intersection point [known to lie on line_pts]
361                    line_pts -- ordered list of [x,y] points making a line.
362                    point_movement_threshold -- used to decide to move or add intersectionPt
363   
364            OUTPUT:
365                new version of lint_pts with the point added
366        """
367        # Avoid pointer/copy issues
368        L1_pts = copy.copy(line_pts)
369        iP=copy.copy(intersectionPt)
370   
371        # Adjust L1
372        tmp = find_nearest_segment(intersectionPt, L1_pts)
373        # Compute the distance from the end points of the segment to the
374        # intersection point. Based on this we decide to add or move the point
375        p0 = L1_pts[tmp[1]]
376        p1 = L1_pts[tmp[1]+1]
377        endPt_Dist_Sq=min( ( (p0[0]-iP[0])**2 + (p0[1]-iP[1])**2), ( (p1[0]-iP[0])**2 + (p1[1]-iP[1])**2))
378        #
379        if(endPt_Dist_Sq>point_movement_threshold**2):
380            # Insert the intersection point. We do this in a tricky way to
381            # account for the possibility of L1_pts having > 2 coordinates
382            # (riverWalls)
383            if verbose:
384                print '      Inserting new point'
385            dummyPt=copy.copy(L1_pts[tmp[1]])
386            L1_pts.insert(tmp[1]+1,dummyPt) 
387            L1_pts[tmp[1]+1][0]=iP[0]
388            L1_pts[tmp[1]+1][1]=iP[1]
389            if(len(L1_pts[tmp[1]+1])==3):
390                # Treat 3rd coordinate
391                # Find distance of inserted point from neighbours, and
392                # Set 3rd coordinate as distance-weighted average of the others
393                d0=((L1_pts[tmp[1]][0]-L1_pts[tmp[1]+1][0])**2.+\
394                   (L1_pts[tmp[1]][1]-L1_pts[tmp[1]+1][1])**2.)**0.5
395                d1=((L1_pts[tmp[1]+2][0]-L1_pts[tmp[1]+1][0])**2.+\
396                   (L1_pts[tmp[1]+2][1]-L1_pts[tmp[1]+1][1])**2.)**0.5
397                L1_pts[tmp[1]+1][2] = (d0*L1_pts[tmp[1]+2][2] + d1*L1_pts[tmp[1]][2])/(d0+d1)
398   
399        else:
400            if verbose:
401                print '      Shifting existing point'
402            # Move a point already on L1
403            L1_pts=shift_point_on_line(iP, L1_pts, tmp[1])
404   
405        return L1_pts
406   
407    #######################################################################################################
408    def check_polygon_is_small(intersection, buf, tol2=100.):
409        """
410           
411            Elsewhere in the code, we check whether lines intersect by buffering them
412            to polygons with a small buffer = buf, then getting the intersection.
413             [since intersection with polygons is supported by gdal, but apparently
414            not directly with lines]. 
415           
416            The logic of our code only works with point intersections, 
417            and it will fails if 2 lines overlap in a line.
418   
419            We crudely check for this situation here, by ensuring that the intersection polygon is 'small'
420   
421            WARNING: The gdal geometry routines may be a bit rough (?)
422                     Intersections not very precise, etc (?)
423   
424            INPUT: intersection -- intersection of the 2 lines [gdal geometry]
425                   buf -- a length scale giving the size of the intersection extent that we expect for a point
426                   tol2 -- Throw an error if the x or y extent is greater than buf*tol2. Seems this needs to be
427                           large sometimes -- this might reflect the stated weaknesses of GDALs geometry routines?
428            OUTPUT: True/False
429                    False should suggest that the intersection is not describing a point
430        """
431   
432        extent=intersection.GetEnvelope()
433        assert(len(extent)==4) # Make sure this assumption is valid
434        if( (abs(extent[0]-extent[1])>buf*tol2) or (abs(extent[2]-extent[3]) > buf*tol2)):
435            return False
436        else:
437            return True
438   
439    #######################################################################################################
440   
441    def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0, buf=1.0e-05, tol2 = 100,
442                                  verbose=True, nameFlag=''):
443        """
444            Add intersection points to lines L1 and L2 if they intersect each other
445   
446            This is useful e.g. so that intersections can be exact (important for
447                 mesh generation in ANUGA)
448   
449            It currently only supports point intersection of 2 lines.
450                Line intersections should fail gracefully
451           
452            INPUTS:  L1, L2 = Wkb LineString geometries
453   
454                     point_movement_threshold -- if the distance from the nearest
455                        point on L1 or L2 to the intersection is < point_movement_threshold, then the
456                        nearest point has its coordinates replaced with the intersection point. This is
457                        to prevent points being too close to each other
458       
459                     buf = tolerence that is used to buffer lines to find
460                            intersections. Probably doesn't need modification
461   
462                     tol2 = [see check_polygon_is_small] Probably doesn't need to change
463   
464                     nameFlag = This will be printed if intersection occurs. Convenient way to display intersecting filenames
465           
466            OUTPUTS: L1,L2 with intersection points added in the right places
467        """
468   
469        if(L1.Intersects(L2)):
470            # Get points on the lines
471            L1_pts=Wkb2ListPts(L1) #L1.GetPoints()
472            L2_pts=Wkb2ListPts(L2) #L2.GetPoints()
473           
474            # Buffer lines by a small amount
475            L1_buf=L1.Buffer(buf)
476            L2_buf=L2.Buffer(buf)
477               
478            # Get intersection point[s]   
479            L1_L2_intersect=L1_buf.Intersection(L2_buf)
480            if(L1_L2_intersect.GetGeometryCount()==1):
481                if(not check_polygon_is_small(L1_L2_intersect, buf, tol2)):
482                    #print L1_L2_intersect.GetEnvelope()
483                    raise Exception, 'line intersection is not allowed. Envelope %s '% str(L1_L2_intersect.GetEnvelope())
484                # Seems to need special treatment with only 1 intersection point
485                intersectionPts=[L1_L2_intersect.Centroid().GetPoint()]
486            else:
487                intersectionPts=[]
488                for feature in L1_L2_intersect:
489                    if(not check_polygon_is_small(feature, buf, tol2)):
490                        print feature.GetEnvelope()
491                        raise Exception, 'line intersection is not allowed'
492                    intersectionPts.append(feature.Centroid().GetPoint())
493   
494            if(verbose):
495                print nameFlag
496                print '    Treating intersections in ', len(intersectionPts) , ' locations'
497                print intersectionPts
498   
499            # Insert the points into the line segments
500            for i in range(len(intersectionPts)):
501                L1_pts = insert_intersection_point(intersectionPts[i], L1_pts, point_movement_threshold, verbose=verbose)
502                L2_pts = insert_intersection_point(intersectionPts[i], L2_pts, point_movement_threshold, verbose=verbose)
503               
504            # Convert to the input format
505            L1_pts=ListPts2Wkb(L1_pts,geometry_type='line')
506            L2_pts=ListPts2Wkb(L2_pts,geometry_type='line')
507   
508            return [L1_pts, L2_pts]
509        else:
510            return [L1, L2]
511   
512    ###########################################################
513    def getRasterExtent(rasterFile, asPolygon=False):
514        """
515            Sometimes we need to know the extent of a raster
516            i.e. the minimum x, maximum x, minimum y, and maximum y values
517           
518            INPUT:
519                rasterFile -- a gdal compatible rasterfile
520                asPolygon -- if False, return [xmin,xmax,ymin,ymax].
521                             If True, return [ [xmin,ymin],[xmax,ymin],[xmax,ymax],[xmin,ymax]]
522            OUTPUT
523                The extent as defined above
524
525        """
526        raster = gdal.Open(rasterFile)
527        transform=raster.GetGeoTransform()
528        xOrigin = transform[0]
529        yOrigin = transform[3]
530        xPixels=raster.RasterXSize
531        yPixels=raster.RasterYSize
532
533        # Compute the other extreme corner
534        x2=xOrigin + xPixels*transform[1]+yPixels*transform[2]
535        y2=yOrigin + xPixels*transform[4]+yPixels*transform[5]
536       
537        xmin=min(xOrigin,x2) 
538        xmax=max(xOrigin,x2)
539
540        ymin=min(yOrigin,y2)
541        ymax=max(yOrigin,y2)
542
543        if(asPolygon):
544            return [ [xmin,ymin], [xmax,ymin], [xmax,ymax], [xmin,ymax]]
545        else:
546            return [xmin,xmax,ymin,ymax]
547
548 
549    ###########################################################
550    def rasterValuesAtPoints(xy, rasterFile, band=1):
551        """
552            Get raster values at point locations.
553            Can be used to e.g. set quantity values
554       
555            INPUT:
556            xy = numpy array with point locations
557            rasterFile = Filename of the gdal-compatible raster
558            band = band of the raster to get
559   
560            OUTPUT:
561            1d numpy array with raster values at xy
562   
563        """
564        # Raster info
565        raster = gdal.Open(rasterFile)
566        rasterBand=raster.GetRasterBand(band)
567        rasterBandType=gdal.GetDataTypeName(rasterBand.DataType)
568   
569        # Projection info
570        transform=raster.GetGeoTransform()
571        xOrigin = transform[0]
572        yOrigin = transform[3]
573        pixelWidth = transform[1]
574        pixelHeight = transform[5] # Negative
575       
576        # Get coordinates in pixel values
577        px = ((xy[:,0] - xOrigin) / pixelWidth).astype(int) #x
578        py = ((xy[:,1] - yOrigin) / pixelHeight).astype(int) #y
579     
580        # Hold elevation
581        elev=px*0. 
582   
583        # Get the right character for struct.unpack
584        if (rasterBandType == 'Int16'):
585            CtypeName='h'
586        elif (rasterBandType == 'Float32'):
587            CtypeName='f'
588        elif (rasterBandType == 'Byte'):
589            CtypeName='B'
590        else:
591            print 'unrecognized DataType:', gdal.GetDataTypeName(band.DataType)
592            print 'You might need to edit this code to read the data type'
593            raise Exception, 'Stopping'
594 
595        # Upper bounds for pixel values, so we can fail gracefully
596        xMax=raster.RasterXSize
597        yMax=raster.RasterYSize
598        if(px.max()<xMax and px.min()>=0 and py.max()<yMax and py.min()>=0):
599            pass
600        else:
601            raise Exception, 'Trying to extract point values that exceed the raster extent'
602
603        # Get values -- seems we have to loop, but it is efficient enough
604        for i in range(len(px)):
605            xC=int(px[i])
606            yC=int(py[i])
607            structval=rasterBand.ReadRaster(xC,yC,1,1,buf_type=rasterBand.DataType)
608            elev[i] = struct.unpack(CtypeName, structval)[0]
609   
610        return elev
611   
612   
613    def gridPointsInPolygon(polygon, approx_grid_spacing=[1.,1.], eps=1.0e-06):
614        """
615            Get a 'grid' of points inside a polygon. Could be used with rasterValuesAtPoints to
616               get a range of raster values inside a polygon
617   
618            Approach: A 'trial-grid' of points is created which is 'almost' covering the range of the polygon (xmin-xmax,ymin-ymax).
619   
620                      (Actually it is just inside this region, to avoid polygon-boundary issues, see below)
621   
622                      Then we find those points which are actually inside the polygon.
623   
624                      The x/y point spacing on the trial-grid will be close to
625                      approx_grid_spacing, but we ensure there are at least 3x3 points on the trial grid.
626                      Also, we reduce the spacing so that the min_x+R and max_x-R
627                        are both similarly close to the polygon extents [see definition of R below]
628   
629            INPUTS:
630                polygon -- the polygon in ANUGA format (list of lists of ordered xy points)
631   
632                approx_grid_spacing -- the approximate x,y grid spacing
633   
634                eps -- 'trial-grid' of points has x range from min_polygon_x+R to
635                        max_polygon_x - R, where R = (max_polygon_x-min_polygon_x)*eps
636                        ( and similarly for y).
637   
638                       This makes it more likely that points are inside the
639                        polygon, not on the boundaries. Points on the boundaries can confuse the
640                        point-in-polygon routine
641   
642            OUTPUTS: A n x 2 numpy array of points in the polygon
643        """
644   
645        # Get polygon extent
646        polygonArr=numpy.array(polygon)
647        poly_xmin=polygonArr[:,0].min()
648        poly_xmax=polygonArr[:,0].max()
649        poly_ymin=polygonArr[:,1].min()
650        poly_ymax=polygonArr[:,1].max()
651   
652        # Make a 'grid' of points which covers the polygon
653        xGridCount=max( numpy.ceil( (poly_xmax-poly_xmin)/approx_grid_spacing[0]+1. ).astype(int), 3)
654        R=(poly_xmax-poly_xmin)*eps
655        Xvals=numpy.linspace(poly_xmin+R,poly_xmax-R, xGridCount)
656        yGridCount=max( numpy.ceil( (poly_ymax-poly_ymin)/approx_grid_spacing[1]+1. ).astype(int), 3)
657        R=(poly_ymax-poly_ymin)*eps
658        Yvals=numpy.linspace(poly_ymin+R,poly_ymax-R, yGridCount)
659   
660        xGrid,yGrid=numpy.meshgrid(Xvals,Yvals)
661        Grid=numpy.vstack([xGrid.flatten(),yGrid.flatten()]).transpose()
662   
663        # Now find the xy values which are actually inside the polygon
664        polpath=path.Path(polygonArr)
665        keepers=[]
666        for i in range(len(Grid[:,0])):
667           if(polpath.contains_point(Grid[i,:])):
668                keepers=keepers+[i]
669        #import pdb
670        #pdb.set_trace()
671        xyInside=Grid[keepers,:]
672   
673        return(xyInside)
674   
675    #########################################################################
676    # Function to search for pattern matches in a string (turns out to be useful)
677    def matchInds(pattern, stringList):
678        """
679            Find indices in stringList which match pattern
680        """
681        #matches=[ (pattern in stringList[i]) for i in range(len(stringList))]
682        matches=[]
683        for i in range(len(stringList)):
684            if pattern in stringList[i]:
685                matches.append(i)
686        return matches
687   
688   
689    ###########################################################################
690    #
691    # Less generic utilities below
692    #
693    # These are more 'anuga-specific' than above, aiming to make nice interfaces
694    # in ANUGA scripts
695    #
696    ############################################################################
697   
698   
699    def add_intersections_to_domain_features(bounding_polygonIn,
700                breakLinesIn={ }, riverWallsIn={ }, point_movement_threshold=0.,
701                verbose=True):
702        """
703            If bounding polygon / breaklines /riverwalls intersect with each other, then
704            add intersection points.
705   
706            INPUTS:
707                bounding_polygonIn -- the bounding polygon in ANUGA format
708                breakLinesIn -- the breaklines dictionary
709                riverWallsIn -- the riverWalls dictionary
710   
711                point_movement_threshold -- if points on lines are < this distance
712                        from intersection points, then they are replaced with the intersection point.
713                        This can prevent overly close points from breaking the mesh generation
714   
715            OUTPUT:
716                List with bounding_polygon,breakLines,riverwalls
717        """
718   
719        bounding_polygon=copy.copy(bounding_polygonIn)
720        breakLines=copy.copy(breakLinesIn)
721        riverWalls=copy.copy(riverWallsIn)
722   
723        # Clean intersections of breakLines with itself
724        if(verbose): 
725            print 'Cleaning breakline intersections'
726        if(len(breakLines)>0):
727            kbl=breakLines.keys()
728            for i in range(len(kbl)):
729                n1=kbl[i]
730                for j in range(len(kbl)):
731                    if(i>=j):
732                        continue
733                    n2=kbl[j]
734                    # Convert breaklines to wkb format
735                    bl1=ListPts2Wkb(breakLines[n1],geometry_type='line')
736                    bl2=ListPts2Wkb(breakLines[n2],geometry_type='line')
737                    # Add intersection points
738                    bl1, bl2 =addIntersectionPtsToLines(bl1, bl2,\
739                                    point_movement_threshold=point_movement_threshold,
740                                    verbose=verbose, nameFlag=n1+' intersects '+ n2)
741                    breakLines[n1]=Wkb2ListPts(bl1)
742                    breakLines[n2]=Wkb2ListPts(bl2)
743                   
744   
745        # Clean intersections of riverWalls with itself
746        if(verbose): 
747            print 'Cleaning riverWall intersections'
748        if(len(riverWalls)>0):
749            krw=riverWalls.keys()
750            for i in range(len(krw)):
751                n1=krw[i]
752                for j in range(len(krw)):
753                    if(i>=j):
754                        continue 
755                    n2=krw[j]
756                    # Convert breaklines to wkb format
757                    rw1=ListPts2Wkb(riverWalls[n1],geometry_type='line')
758                    rw2=ListPts2Wkb(riverWalls[n2],geometry_type='line')
759                    # Add intersection points
760                    rw1, rw2 =addIntersectionPtsToLines(rw1, rw2,\
761                                    point_movement_threshold=point_movement_threshold,\
762                                    verbose=verbose, nameFlag=n1+' intersects '+ n2)
763                    riverWalls[n1]=Wkb2ListPts(rw1)
764                    riverWalls[n2]=Wkb2ListPts(rw2)
765   
766        # Clean intersections of breaklines with riverwalls
767        if(verbose): 
768            print 'Cleaning breakLine-riverWall intersections'
769        if( (len(riverWalls)>0) and (len(breakLines)>0)):
770            krw=riverWalls.keys()
771            kbl=breakLines.keys()
772            for i in range(len(krw)):
773                n1=krw[i]
774                for j in range(len(kbl)):
775                    n2=kbl[j]
776                    # Convert breaklines to wkb format
777                    rw1=ListPts2Wkb(riverWalls[n1],geometry_type='line')
778                    bw2=ListPts2Wkb(breakLines[n2],geometry_type='line')
779                    # Add intersection points
780                    rw1, bw2 =addIntersectionPtsToLines(rw1, bw2,\
781                                    point_movement_threshold=point_movement_threshold,\
782                                    verbose=verbose, nameFlag=n1+' intersects '+ n2)
783                    riverWalls[n1]=Wkb2ListPts(rw1)
784                    breakLines[n2]=Wkb2ListPts(bw2)
785                   
786   
787        # Clean intersections of bounding polygon and riverwalls
788        if(verbose): 
789            print 'Cleaning bounding_poly-riverWall intersections'
790        if( (len(riverWalls)>0)):
791            krw=riverWalls.keys()
792            for i in range(len(krw)):
793                n1=krw[i]
794                # Convert breaklines to wkb format
795                rw1=ListPts2Wkb(riverWalls[n1],geometry_type='line')
796                bp2=ListPts2Wkb(bounding_polygon,geometry_type='line', appendFirstOnEnd=True)
797                # Add intersection points
798                rw1, bp2 =addIntersectionPtsToLines(rw1, bp2,\
799                                point_movement_threshold=point_movement_threshold,\
800                                verbose=verbose, nameFlag='Bounding Pol intersects '+ n1)
801                riverWalls[n1]=Wkb2ListPts(rw1)
802                # Since the bounding polygon is a loop, the first/last points are the same
803                # If one of these was moved, the other should be moved too. Since we
804                # will drop the last bounding_polygon point, we only need to worry about the first
805                bounding_polygon=Wkb2ListPts(bp2,removeLast=False)
806                if(bounding_polygon[-1] is not bounding_polygon[0]):
807                    bounding_polygon[0]=bounding_polygon[-1]
808                # Drop the last point
809                bounding_polygon=bounding_polygon[:-1]
810   
811        # Clean intersections of bounding polygon and breaklines
812        if(verbose):
813            print 'Cleaning bounding_poly-breaklines intersections'
814        if( (len(breakLines)>0)):
815            kbl=breakLines.keys()
816            for i in range(len(kbl)):
817                n1=kbl[i]
818                # Convert breaklines to wkb format
819                bl1=ListPts2Wkb(breakLines[n1],geometry_type='line')
820                bp2=ListPts2Wkb(bounding_polygon,geometry_type='line', appendFirstOnEnd=True)
821                # Add intersection points
822                bl1, bp2 =addIntersectionPtsToLines(bl1, bp2,\
823                                point_movement_threshold=point_movement_threshold,
824                                verbose=verbose, nameFlag='Bounding Pol intersects '+n1)
825                breakLines[n1]=Wkb2ListPts(bl1)
826                # Since the bounding polygon is a loop, the first/last points are the same
827                # If one of these was moved, the other should be moved too. Since we
828                # will drop the last bp2 point, we only need to worry about the first
829                bounding_polygon=Wkb2ListPts(bp2,removeLast=False)
830                if(bounding_polygon[-1] is not bounding_polygon[0]):
831                    bounding_polygon[0]=bounding_polygon[-1]
832                # Drop the last point
833                bounding_polygon=bounding_polygon[:-1]
834
835        # Remove the extra 0.0 from bounding polygon [this cannot have 3 coordinates]
836        bounding_polygon = [ bounding_polygon[i][0:2] for i in range(len(bounding_polygon))]
837        # Same for breaklines [although might not matter]
838        for n1 in breakLines.keys():
839            breakLines[n1] = [breakLines[n1][i][0:2] for i in range(len(breakLines[n1]))]
840
841        # Check that all mesh-boundary points are inside the bounding polygon
842        from anuga.geometry.polygon import outside_polygon
843        for blCat in [riverWalls, breakLines]:
844            for n1 in blCat.keys():
845                l=len(blCat[n1])
846                # Test every point -- means we can strip 3rd coordinate if needed
847                for j in range(l):
848                    isOut=outside_polygon(blCat[n1][j][0:2], bounding_polygon)
849                    if(len(isOut)>0):
850                        msg='Breakline/riverwall point '+str(blCat[n1][j][0:2]) +' on '+ n1+\
851                            ' is outside the bounding polygon.\n'+\
852                            'Check that it exceeds the bounding polygon by a distance < point_movement_threshold \n'+\
853                            ' so it can be moved back onto the polygon'
854                        print 'Polygon\n '
855                        print bounding_polygon
856                        print 'Line \n'
857                        print blCat[n1]
858                        raise Exception, msg
859     
860        return [bounding_polygon, breakLines, riverWalls]
861   
862    ###################################################################
863   
864    def readRegionPtAreas(shapefile, convert_length_to_area=False):
865        """
866            Read a point shapefile to define the ANUGA mesh resoutions.
867   
868            MUST HAVE A SINGLE ATTRIBUTE REPRESENTING THE LENGTHS OF TRIANGLES IN
869             REGIONS
870   
871            INPUT: shapefile -- name of the input shapefile
872                   convert_length_to_area -- if True, res values are assumed to
873                          represent triangle side lengths, and are converted to areas with 0.5*res0*res0
874                          Note that this might not ensure that the max triangle side length really is res0, but
875                          it will be of similar magnitude
876                          If False, attribute values are assumed to represent triangle areas
877   
878            OUTPUT: list of the form  [ [x0,y0,res0], [x1, y1, res1], ...]
879        """
880   
881        ptData=readShpPtsAndAttributes(shapefile)
882   
883        # Must have only 1 attribute
884        assert len(ptData[2])==1
885   
886        numPts=len(ptData[0])
887        outData=[]
888        for i in range(numPts):
889            if(convert_length_to_area):
890                newDat=[ptData[0][i][0], ptData[0][i][1], 0.5*float(ptData[1][i])**2]
891            else:
892                newDat=[ptData[0][i][0], ptData[0][i][1], float(ptData[1][i])]
893            outData.append(newDat)
894   
895        return outData
896   
897    #########################################
898    def readListOfBreakLines(shapefileList):
899        """
900            Take a list with the names of shapefiles
901       
902            They are assumed to be '2D breaklines', so we just read their
903                coordinates into a dict with their names
904   
905            Read them in
906           
907            INPUT: shapefileList -- a list of shapefile names [e.g. from glob.glob('GIS/Breaklines/*.shp')]
908   
909            OUTPUT: dictionary with breaklines [filenames are keys]
910        """
911   
912        allBreakLines={}
913        for shapefile in shapefileList:
914            allBreakLines[shapefile]=readShp_1LineGeo(shapefile)
915       
916        return allBreakLines
917   
918    ############################################################################
919    def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
920        """
921            We sometimes have breaklines defining 2 edges of a channel,
922            and wish to make a polygon from them
923   
924            Can do this with the current function
925   
926            INPUTS: pattern == character string containing pattern which is inside exactly 2 keys in breaklines
927       
928                    breakLines = the breakLinesIn dictionary
929                   
930                    reverse2nd = True/False or None. Reverse the order of the 2nd set of edges before making the polygon.
931                                 If None, then we compute the distance between the
932                                 first point on breakline1 and the first/last points on breakline2, and reverse2nd if
933                                 the 'distance from the first point' < 'distance from the last point'
934       
935            OUTPUT: Polygon made with the 2 breaklines
936        """
937   
938        breakLines=copy.copy(breakLinesIn)
939        bk=breakLines.keys()
940        matchers=matchInds(pattern, bk)
941   
942        if(len(matchers)==0):
943            msg = 'Cannot match ' + pattern + 'in breaklines file names'
944            raise Exception, msg
945   
946        if(len(matchers)!=2):
947            print 'Need exactly 2 matches, but pattern matched these', bk[matchers]
948   
949        # There are 2 matches
950       
951        if(reverse2nd is None):
952            # Automatically compute whether we should reverse the 2nd breakline
953            bl1_0=breakLines[bk[matchers[0]]][0]
954            bl2_0=breakLines[bk[matchers[1]]][0]
955            bl2_N=breakLines[bk[matchers[1]]][-1]
956            d0 = ((bl1_0[0]-bl2_0[0])**2 + (bl1_0[1]-bl2_0[1])**2)
957            dN = ((bl1_0[0]-bl2_N[0])**2 + (bl1_0[1]-bl2_N[1])**2)
958            if(d0<dN):
959                reverse2nd = True
960            else:
961                reverse2nd = False
962   
963        if(reverse2nd):
964            breakLines[bk[matchers[1]]].reverse()
965        polyOut=breakLines[bk[matchers[0]]] +  breakLines[bk[matchers[1]]]
966        #polyOut=polyOut+[polyOut[0]]
967        return polyOut
968    ###################
969   
970else: # gdal_available == False
971    msg='Failed to import gdal/ogr modules --'\
972        + 'perhaps gdal python interface is not installed.'
973
974
975   
976    def readShp_1PolyGeo(shapefile, dropLast=True):
977        raise ImportError, msg
978   
979    def readShp_1LineGeo(shapefile):
980        raise ImportError, msg
981   
982    def readShpPtsAndAttributes(shapefile):
983        raise ImportError, msg
984   
985    def ListPts2Wkb( ptsIn, geometry_type='line', appendFirstOnEnd=None):
986        raise ImportError, msg
987   
988    def Wkb2ListPts(wkb_geo, removeLast=False, drop_third_dimension=False):
989        raise ImportError, msg
990   
991    def compute_squared_distance_to_segment(pt, line):
992        raise ImportError, msg
993   
994    def find_nearest_segment(pt, segments):
995        raise ImportError, msg
996   
997    def shift_point_on_line(pt, lineIn, nearest_segment_index):
998        raise ImportError, msg
999   
1000    def insert_intersection_point(intersectionPt, line_pts, point_movement_threshold,verbose=False):
1001        raise ImportError, msg
1002
1003    def check_polygon_is_small(intersection, buf, tol2=100.):
1004        raise ImportError, msg
1005   
1006    def addIntersectionPtsToLines(L1,L2, point_movement_threshold=0.0, buf=1.0e-06, tol2 = 100,
1007                                  verbose=True, nameFlag=''):
1008        raise ImportError, msg
1009   
1010
1011    def rasterValuesAtPoints(xy, rasterFile, band=1):
1012        raise ImportError, msg
1013   
1014   
1015    def gridPointsInPolygon(polygon, approx_grid_spacing=[1.,1.], eps=1.0e-06):
1016        raise ImportError, msg
1017   
1018
1019    def matchInds(pattern, stringList):
1020        raise ImportError, msg
1021   
1022   
1023    def add_intersections_to_domain_features(bounding_polygonIn,
1024                breakLinesIn={ }, riverWallsIn={ }, point_movement_threshold=0.,
1025                verbose=True):
1026        raise ImportError, msg
1027   
1028   
1029    def readRegionPtAreas(shapefile, convert_length_to_area=False):
1030        raise ImportError, msg
1031   
1032    def readListOfBreakLines(shapefileList):
1033        raise ImportError, msg
1034   
1035    def polygon_from_matching_breaklines(pattern,breakLinesIn, reverse2nd=None):
1036        raise ImportError, msg
1037    ###################   
1038
1039
Note: See TracBrowser for help on using the repository browser.