Changeset 2261 for inundation/pmesh

Timestamp:

Jan 20, 2006, 10:59:35 AM (19 years ago)

Author:

duncan

Message:

added geo-ref to poly to region method

Location:

inundation/pmesh

Files:

• mesh.py (modified) (25 diffs)
• test_mesh.py (modified) (4 diffs)

inundation/pmesh/mesh.py

@@ -358 +358 @@
         return self.vertices
     
+    def get_vertices(self):
+        """
+        Return a list of the vertices.  The x and y values will be absolute
+        Easting and Northings for the zone of the current geo_ref.
+        """
+        return self.vertices
+    
     def calcArea(self):
         ax = self.vertices[0].x
@@ -571 +578 @@
         
         return (dic.__cmp__(dic_other))
-        
-    def generateMesh(self, mode = None, maxArea = None, isRegionalMaxAreas = True):
+    
+    def addUserPoint(self, pointType, x,y):
+        if pointType == Vertex:
+            point = self.addUserVertex(x,y)
+        if pointType == Hole:
+            point = self._addHole(x,y)
+        if pointType == Region:
+            point = self._addRegion(x,y)
+        return point
+    
+    def addUserVertex(self, x,y):
+        v=Vertex(x, y)
+        self.userVertices.append(v)
+        return v
+
+    def _addHole(self, x,y):
+        h=Hole(x, y)
+        self.holes.append(h)
+        return h
+   
+    def _addRegion(self, x,y):
+        h=Region(x, y)
+        self.regions.append(h)
+        return h
+   
+    def add_region(self, x,y, geo_reference=None):
+        [x,y] = self.geo_reference.change_points_geo_ref([x,y],
+                                                 points_geo_ref=geo_reference)
+        return self._addRegion(x, y)
+
+    # Depreciated
+    def addRegionEN(self, x,y):
+        print "depreciated, use add_region"
+        return self.add_region(x,y)
+    
+    def add_region_from_polygon(self, polygon, tags=None,
+                             max_area=None, geo_reference=None):
+        """
+        Add a polygon with tags to the current mesh, as a region.
+        The maxArea of the region can be specified.
+
+        If a geo_ref is given, this is used.
+        If not;
+        The x,y info is assumed to be Easting and Northing, absolute,
+        for the meshes zone.
+
+        polygon a list of points, in meters that describe the polygon
+             (e.g. [[x1,y1],[x2,y2],...]
+        tags (e.g.{'wall':[0,1,3],'ocean':[2]})
+
+        This returns the region instance, so if the user whats to modify
+        it they can.
+        
+        """
+        #FIXME: take into account georef on mesh side
+        polygon = self.geo_reference.change_points_geo_ref(polygon,
+                                                 points_geo_ref=geo_reference)
+        #print "polygon - should be relative to mesh geo_ref",polygon 
+        #create points, segs and tags
+        region_dict = {}
+        region_dict['points'] = polygon
+        
+        #Create segments
+        #E.g. [[0,1], [1,2], [2,3], [3,0]]
+        #from polygon
+        #[0,1,2,3]
+        segments = []
+        N = len(polygon)
+        for i in range(N):
+            lo = i
+            hi = (lo + 1) % N
+            segments.append( [lo, hi] ) 
+        region_dict['segments'] = segments
+
+
+        #Create tags
+        #E.g. ['wall', 'wall', 'ocean', 'wall']
+        # from a dic
+        #{'wall':[0,1,3],'ocean':[2]}
+        segment_tags = ['']*N
+        if tags is not None:
+            for key in tags:
+                indices = tags[key]
+                for i in indices:
+                    segment_tags[i] = key
+        region_dict['segment_tags'] = segment_tags
+    
+        self.addVertsSegs(region_dict) #this is assuming absolute geos
+    
+        #get inner point
+        inner_point = point_in_polygon(polygon)
+        inner = self.add_region(inner_point[0], inner_point[1],
+                               geo_reference=self.geo_reference)
+
+        if max_area is not None:
+            inner.setMaxArea(max_area)
+        return inner
+        
+    def getUserVertices(self):
+        return self.userVertices
+    
+    def getUserSegments(self):
+        allSegments = self.userSegments + self.alphaUserSegments
+        #print "self.userSegments",self.userSegments
+        #print "self.alphaUserSegments",self.alphaUserSegments
+        #print "allSegments",allSegments
+        return allSegments
+    
+    def deleteUserSegments(self,seg):
+        if self.userSegments.count(seg) == 0:
+            self.alphaUserSegments.remove(seg)
+            pass
+        else:
+            self.userSegments.remove(seg)
+            
+    def clearUserSegments(self):
+        self.userSegments = []
+        self.alphaUserSegments = []
+                
+    def getTriangulation(self):
+        return self.meshTriangles
+    
+    def getMeshVertices(self):
+        return self.meshVertices
+ 
+    def getMeshSegments(self):
+        return self.meshSegments
+    
+    def getHoles(self):
+        return self.holes
+    
+    def getRegions(self):
+        return self.regions
+    
+    def isTriangulation(self):
+        if self.meshVertices == []:
+            return False 
+        else:
+            return True
+    
+    def addUserSegment(self, v1,v2):
+        """
+        PRECON: A segment between the two vertices is not already present.
+        Check by calling isUserSegmentNew before calling this function.
+        
+        """
+        s=Segment( v1,v2)
+        self.userSegments.append(s)
+        return s
+        
+    def generateMesh(self, mode = None, maxArea = None,
+                     isRegionalMaxAreas = True):
         """
         Based on the current user vaules, holes and regions
@@ -592 +749 @@
             # it's more comlex than this.  eg holes
         if not re.match('z',self.mode):
-            self.mode += 'z' # z - Number all items starting from zero (rather than one)
+            self.mode += 'z' # z - Number all items starting from zero
+                             # (rather than one)
         if not re.match('n',self.mode):
             self.mode += 'n' # n - output a list of neighboring triangles 
@@ -631 +789 @@
         #print "generated",generatedMesh
         generatedMesh['generatedsegmentmarkerlist'] = \
-                     segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
+             segment_ints2strings(generatedMesh['generatedsegmentmarkerlist'],
                                   segconverter)
         #print "processed gen",generatedMesh['generatedsegmentmarkerlist']
@@ -641 +799 @@
         if len(generatedMesh['generatedpointattributelist'][0])==0:
             self.attributeTitles = []
-        generatedMesh['generatedpointattributetitlelist']=self.attributeTitles     
+        generatedMesh['generatedpointattributetitlelist']= \
+                                            self.attributeTitles
 
         self.setTriangulation(generatedMesh)
-    
-    def addUserPoint(self, pointType, x,y):
-        if pointType == Vertex:
-            point = self.addUserVertex(x,y)
-        if pointType == Hole:
-            point = self.addHole(x,y)
-        if pointType == Region:
-            point = self.addRegion(x,y)
-        return point
-    
-    def addUserVertex(self, x,y):
-        v=Vertex(x, y)
-        self.userVertices.append(v)
-        return v
-
-    def addHole(self, x,y):
-        h=Hole(x, y)
-        self.holes.append(h)
-        return h
-   
-    def addRegion(self, x,y):
-        h=Region(x, y)
-        self.regions.append(h)
-        return h
-
-    #FIXME(DSG-DSG) remove EN, have a relative flag. 
-    def addRegionEN(self, x,y):
-        h=Region(x-self.geo_reference.xllcorner,
-                 y-self.geo_reference.yllcorner)
-        self.regions.append(h)
-        return h
-    def addRegionFromPolygon(self, polygon, tags=None,
-                             maxArea=None, Geo_ref=None):
-        #take into account georef
-        #create points, segs and tags
-        region_dict = {}
-        region_dict['points'] = polygon
-
-        #Create segments
-        #E.g. [[0,1], [1,2], [2,3], [3,0]]
-        segments = []
-        N = len(polygon)
-        for i in range(N):
-            lo = i
-            hi = (lo + 1) % N
-            segments.append( [lo, hi] ) 
-        region_dict['segments'] = segments
-
-
-        #Create tags
-        #E.g. ['wall', 'wall', 'ocean', 'wall']
-
-        segment_tags = ['']*N
-        if tags is not None:
-            for key in tags:
-                indices = tags[key]
-                for i in indices:
-                    segment_tags[i] = key
-        region_dict['segment_tags'] = segment_tags
-    
-        self.addVertsSegs(region_dict) #this is assuming absolute geos
-    
-        #get inner point
-        inner_point = point_in_polygon(polygon)
-        inner = self.addRegionEN(inner_point[0], inner_point[1])
-
-        if maxArea is not None:
-            inner.setMaxArea(maxArea)
-        
-    def getUserVertices(self):
-        return self.userVertices
-    
-    def getUserSegments(self):
-        allSegments = self.userSegments + self.alphaUserSegments
-        #print "self.userSegments",self.userSegments
-        #print "self.alphaUserSegments",self.alphaUserSegments
-        #print "allSegments",allSegments
-        return allSegments
-    
-    def deleteUserSegments(self,seg):
-        if self.userSegments.count(seg) == 0:
-            self.alphaUserSegments.remove(seg)
-            pass
-        else:
-            self.userSegments.remove(seg)
-            
-    def clearUserSegments(self):
-        self.userSegments = []
-        self.alphaUserSegments = []
-                
-    def getTriangulation(self):
-        return self.meshTriangles
-    
-    def getMeshVertices(self):
-        return self.meshVertices
- 
-    def getMeshSegments(self):
-        return self.meshSegments
-    
-    def getHoles(self):
-        return self.holes
-    
-    def getRegions(self):
-        return self.regions
-    
-    def isTriangulation(self):
-        if self.meshVertices == []:
-            return False 
-        else:
-            return True
-    
-    def addUserSegment(self, v1,v2):
-        """
-        PRECON: A segment between the two vertices is not already present.
-        Check by calling isUserSegmentNew before calling this function.
-        
-        """
-        s=Segment( v1,v2)
-        self.userSegments.append(s)
-        return s
     
     def clearTriangulation(self):
@@ -777 +816 @@
         """
         assert self.getUserSegments() == []
-        self.userVertices, counter =  self.removeDuplicatedVertices(self.userVertices)
+        self.userVertices, counter =  self.removeDuplicatedVertices(
+            self.userVertices)
         return counter
     
@@ -852 +892 @@
             
     def isUserSegmentNew(self, v1,v2):
-        identicalSegs= [x for x in self.getUserSegments() if (x.vertices[0] == v1 and x.vertices[1] == v2) or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
+        identicalSegs= [x for x in self.getUserSegments() \
+                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
+        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
         
         return len(identicalSegs) == 0
@@ -897 +939 @@
                         regions=None):
         """
-        Convert the Mesh to a dictionary of the lists needed for the triang modul;
+        Convert the Mesh to a dictionary of the lists needed for the
+        triang module
         points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
         pointattributelist: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
         segment list: [(point1,point2),(p3,p4),...] (Tuples of integers) 
-        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
+        hole list: [(x1,y1),...](Tuples of doubles, one inside each hole)
         regionlist: [ (x1,y1,tag, max area),...] (Tuple of 3-4 doubles)
         
@@ -960 +1003 @@
     def Mesh2MeshList(self):
         """
-        Convert the Mesh to a dictionary of lists describing the triangulation variables;
+        Convert the Mesh to a dictionary of lists describing the
+        triangulation variables;
         generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)
-        generated point attribute list: [(a11,a12,...),(a21,a22),...] (Tuples of doubles)
-        generated point attribute title list:[A1Title, A2Title ...] (list of strings)
-        generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers) 
+        generated point attribute list: [(a11,a12,...),(a21,a22),...]
+            (Tuples of doubles)
+        generated point attribute title list:[A1Title, A2Title ...]
+            (list of strings)
+        generated segment list: [(point1,point2),(p3,p4),...]
+            (Tuples of integers) 
         generated segment tag list: [tag,tag,...] list of strings
 
@@ -971 +1018 @@
         generated triangle attribute list: [s1,s2,...] list of strings
 
-        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....] tuple of triangles
+        generated triangle neighbor list: [(t1,t2,t3), (t4,t5,t6),....]
+            tuple of triangles
         
         Used to produce .tsh file
@@ -1010 +1058 @@
         triangleneighborlist = []
         for tri in self.meshTriangles:  
-            trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index))  
+            trianglelist.append((tri.vertices[0].index,tri.vertices[1].index,
+                                 tri.vertices[2].index))  
             triangleattributelist.append([tri.attribute])
             neighborlist = [-1,-1,-1]
@@ -1046 +1095 @@
         returned from the triang module       
         generated point list: [(x1,y1),(x2,y2),...] (Tuples of doubles)  
-        generated point attribute list:[(P1att1,P1attt2, ...),(P2att1,P2attt2,...),...]
-        generated point attribute title list:[A1Title, A2Title ...] (list of strings)
-        generated segment list: [(point1,point2),(p3,p4),...] (Tuples of integers)
+        generated point attribute list:[(P1att1,P1attt2, ...),
+            (P2att1,P2attt2,...),...]
+        generated point attribute title list:[A1Title, A2Title ...]
+            (list of strings)
+        generated segment list: [(point1,point2),(p3,p4),...]
+            (Tuples of integers)
         generated segment marker list: [S1Tag, S2Tag, ...] (list of ints)
-        triangle list:  [(point1,point2, point3),(p5,p4, p1),...] (Tuples of integers)
-        triangle neighbor list: [(triangle1,triangle2, triangle3),(t5,t4, t1),...] (Tuples of integers) -1 means there's no triangle neighbor
-        triangle attribute list: [(T1att), (T2att), ...] (list of a list of strings)
+        triangle list:  [(point1,point2, point3),(p5,p4, p1),...]
+            (Tuples of integers)
+        triangle neighbor list: [(triangle1,triangle2, triangle3),
+            (t5,t4, t1),...] (Tuples of integers)
+            -1 means there's no triangle neighbor
+        triangle attribute list: [(T1att), (T2att), ...]
+            (list of a list of strings)
         """
         #Clear the current generated mesh values
@@ -1074 +1130 @@
 
         index = 0
-        for seg,marker in map(None,genDict['generatedsegmentlist'],genDict['generatedsegmentmarkerlist']):
+        for seg,marker in map(None,genDict['generatedsegmentlist'],
+                              genDict['generatedsegmentmarkerlist']):
             segObject = Segment( self.meshVertices[seg[0]],
                            self.meshVertices[seg[1]], tag = marker )
@@ -1115 +1172 @@
                 else:
                     ObjectNeighbor.append(None)
-            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
+            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
+                                                   ObjectNeighbor[1],
+                                                   ObjectNeighbor[2])
             index += 1
 
@@ -1150 +1209 @@
 
         #index = 0
-        for seg,tag in map(None,genDict['segmentlist'],genDict['segmenttaglist']):
+        for seg,tag in map(None,genDict['segmentlist'],
+                           genDict['segmenttaglist']):
             segObject = Segment( self.userVertices[seg[0]],
                            self.userVertices[seg[1]], tag = tag )
@@ -1425 +1485 @@
     
     def representedAlphaUserSegment(self, v1,v2):
-        identicalSegs= [x for x in self.alphaUserSegments if (x.vertices[0] == v1 and x.vertices[1] == v2) or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
+        identicalSegs= [x for x in self.alphaUserSegments \
+                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
+        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
 
         if identicalSegs == []:
@@ -1434 +1496 @@
     
     def representedUserSegment(self, v1,v2):
-        identicalSegs= [x for x in self.userSegments if (x.vertices[0] == v1 and x.vertices[1] == v2) or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
+        identicalSegs= [x for x in self.userSegments \
+                        if (x.vertices[0] == v1 and x.vertices[1] == v2)
+        or (x.vertices[0] == v2 and x.vertices[1] == v1) ]
 
         if identicalSegs == []:
@@ -1515 +1579 @@
     def boxsizeVerts(self):
         """
-        Returns a list of verts denoting a box or triangle that contains verts on the xmin, ymin, xmax and ymax axis.
+        Returns a list of verts denoting a box or triangle that contains
+        verts on the xmin, ymin, xmax and ymax axis.
         Structure: list of verts 
         """
@@ -1549 +1614 @@
                 ymax = vertex.y
                 ymaxVert = vertex
-        verts, count = self.removeDuplicatedVertices([xminVert,xmaxVert,yminVert,ymaxVert])
+        verts, count = self.removeDuplicatedVertices([xminVert,
+                                                      xmaxVert,
+                                                      yminVert,
+                                                      ymaxVert])
          
         return verts
@@ -1555 +1623 @@
     def boxsize(self):
         """
-        Returns a list denoting a box that contains the entire structure of vertices
+        Returns a list denoting a box that contains the entire
+        structure of vertices
         Structure: [xmin, ymin, xmax, ymax] 
         """
@@ -1589 +1658 @@
     def maxMinVertAtt(self, iatt):
         """
-        Returns a list denoting a box that contains the entire structure of vertices
+        Returns a list denoting a box that contains the entire structure
+        of vertices
         Structure: [xmin, ymin, xmax, ymax] 
         """
@@ -1798 +1868 @@
         triangle_neighbors = []
         for tri in self.meshTriangles:  
-            triangles.append([tri.vertices[0].index,tri.vertices[1].index,tri.vertices[2].index])  
+            triangles.append([tri.vertices[0].index,
+                              tri.vertices[1].index,
+                              tri.vertices[2].index])  
             triangle_tags.append(tri.attribute)
             neighborlist = [-1,-1,-1]
@@ -1952 +2024 @@
                 else:
                     ObjectNeighbor.append(None)
-            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],ObjectNeighbor[1],ObjectNeighbor[2])
+            self.meshTriangles[index].setNeighbors(ObjectNeighbor[0],
+                                                   ObjectNeighbor[1],
+                                                   ObjectNeighbor[2])
             index += 1
 
@@ -1980 +2054 @@
 
         #index = 0
-        for seg,tag in map(None,genDict['outline_segments'],genDict['outline_segment_tags']):
+        for seg,tag in map(None,genDict['outline_segments'],
+                           genDict['outline_segment_tags']):
             segObject = Segment( self.userVertices[seg[0]],
                            self.userVertices[seg[1]], tag = tag )
@@ -3504 +3579 @@
 # Create a clear interface. eg
 # have the interface methods more at the top of this file and add comments
+# for the interface functions/methods, use function_name (not functionName),
+
+#Currently
+#function_name methods assume absolute values.  Geo-refs can be passed in.
+#
+
+# instead of functionName
 if __name__ == "__main__":
     #from mesh import *

inundation/pmesh/test_mesh.py

@@ -1655 +1655 @@
 #___________end of Peters tests
 
-    def test_addRegionFromPolygon(self):
+    def test_add_region_from_polygon(self):
         m=Mesh()
-        m.addRegionFromPolygon([[0,0],[1,0],[0,1]])
+        m.add_region_from_polygon([[0,0],[1,0],[0,1]])
         self.failUnless(len(m.regions)==1,
                         'FAILED!')
@@ -1665 +1665 @@
                         'FAILED!')
        
-    def test_addRegionFromPolygon2(self):
+    def test_add_region_from_polygon2(self):
         m=Mesh()
-        m.addRegionFromPolygon([[0,0],[1,0],[1,1],[0,1]],
+        m.add_region_from_polygon([[0,0],[1,0],[1,1],[0,1]],
                                {'tagin':[0,1],'bom':[2]})
         self.failUnless(len(m.regions)==1,
@@ -1685 +1685 @@
         self.failUnless(segs[3].tag=='',
                         'FAILED!') 
-        
-        
-
+       
+    def test_add_region_from_polygon3(self):
+        x=0
+        y=0
+        m=Mesh(geo_reference=Geo_reference(56,x,y))
+        polygon = [[0,0],[1,0],[1,1],[0,1]]
+        x_p = 1000
+        y_p = 4000
+        geo_ref_poly = Geo_reference(56, x_p, y_p)
+        polygon = geo_ref_poly.change_points_geo_ref(polygon)
+        #print "polygon", polygon
+        m.add_region_from_polygon(polygon,
+                               {'tagin':[0,1],'bom':[2]},
+                                  geo_reference=geo_ref_poly)
+        self.failUnless(len(m.regions)==1,
+                        'FAILED!')
+        segs = m.getUserSegments()
+        self.failUnless(len(segs)==4,
+                        'FAILED!')
+        self.failUnless(len(m.userVertices)==4,
+                        'FAILED!') 
+        self.failUnless(segs[0].tag=='tagin',
+                        'FAILED!')  
+        self.failUnless(segs[1].tag=='tagin',
+                        'FAILED!') 
+         
+        self.failUnless(segs[2].tag=='bom',
+                        'FAILED!') 
+        self.failUnless(segs[3].tag=='',
+                        'FAILED!')
+        verts = m.getUserVertices()
+        #print "User verts",verts
+        #print 'polygon',polygon
+        #vert values are relative
+        for point,new_point in map(None,polygon,verts):
+            point_x = point[0] + geo_ref_poly.get_xllcorner()
+            new_point_x = new_point.x + m.geo_reference.get_xllcorner()
+            point_y = point[1] + geo_ref_poly.get_yllcorner()
+            new_point_y = new_point.y + m.geo_reference.get_yllcorner()
+            #print "point_y",point_y
+            #print "new_point_y",new_point_y 
+            
+            self.failUnless(point_x == new_point_x, ' failed')
+            self.failUnless(point_y == new_point_y, ' failed')
+            
+         
+        #REJIG
+    def test_add_region_from_polygon4(self):
+        x=0
+        y=0
+        m=Mesh(geo_reference=Geo_reference(56,x,y))
+        polygon = [[0,0],[1,0],[1,1],[0,1]]
+        
+        m.add_region_from_polygon(polygon,
+                               {'tagin':[0,1],'bom':[2]})
+        self.failUnless(len(m.regions)==1,
+                        'FAILED!')
+        segs = m.getUserSegments()
+        self.failUnless(len(segs)==4,
+                        'FAILED!')
+        self.failUnless(len(m.userVertices)==4,
+                        'FAILED!') 
+        self.failUnless(segs[0].tag=='tagin',
+                        'FAILED!')  
+        self.failUnless(segs[1].tag=='tagin',
+                        'FAILED!') 
+         
+        self.failUnless(segs[2].tag=='bom',
+                        'FAILED!') 
+        self.failUnless(segs[3].tag=='',
+                        'FAILED!')
+        verts = m.getUserVertices()
+        #print "User verts",verts
+        #print 'polygon',polygon
+        #vert values are relative
+        for point,new_point in map(None,polygon,verts):
+            point_x = point[0] 
+            new_point_x = new_point.x + m.geo_reference.get_xllcorner()
+            #print "point_x",point_x
+            #print "new_point_x",new_point_x 
+            point_y = point[1] 
+            new_point_y = new_point.y + m.geo_reference.get_yllcorner()
+            
+            self.failUnless(point_x == new_point_x, ' failed')
+            self.failUnless(point_y == new_point_y, ' failed')
+
+
+       
 def list_comp(A,B):
     yes = len(A)==len(B)
@@ -1695 +1780 @@
     return yes
 
-#___________end of Peters tests
-            
 #-------------------------------------------------------------
 if __name__ == "__main__":