Changeset 2200

Timestamp:

Jan 11, 2006, 5:31:49 PM (19 years ago)

Author:

duncan

Message:

in pmesh.mesh, can now add a region specified as a polygon. (This isn't fully implemented yet)

Location:

inundation

Files:

• pmesh/mesh.py (modified) (8 diffs)
• pmesh/test_mesh.py (modified) (3 diffs)
• pyvolution/coordinate_transforms/geo_reference.py (modified) (1 diff)
• utilities/polygon.py (modified) (1 diff)
• utilities/test_polygon.py (modified) (1 diff)

inundation/pmesh/mesh.py

@@ -28 +28 @@
 import load_mesh
 from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
+from utilities.polygon import point_in_polygon 
 
 SET_COLOUR='red'
@@ -185 +186 @@
     Defined by a point in the hole enclosed by segments.
     """
+
     HOLECORNERLENGTH = 6
+    
     def draw(self, canvas, tags, colour = 'purple',scale = 1, xoffset = 0, yoffset =0, ):
         x =  scale*(self.x + xoffset)
@@ -396 +399 @@
         Null represents no neighbor
         """
-        self.attribute = attribute
+        self.attribute = attribute #this is a string
         
     def __repr__(self):
@@ -501 +504 @@
         self.meshVertices=[]
 
-        #Peters stuff
-        # FIXME (DSG) Sets of what?  
         self.setID={}
         #a dictionary of names.
@@ -667 +668 @@
         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,
@@ -673 +675 @@
         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
@@ -1312 +1351 @@
         segs2delete = self.alphaUserSegments
         #FIXME(DSG-DSG) this algorithm needs comments
-        #FIXME(DSG-DSG) can it be sped up?  It's slow
         new_segs = {}
         #alpha_segs = []
@@ -1603 +1641 @@
         yoffset = HEIGHT - b
         return [SCALE, xoffset, yoffset]
-            
-    def plotMeshTriangle(self,tag = 0,WIDTH = 400,HEIGHT = 400):
-        """
-        Plots all node connections.
-        tag = 0 (no node numbers), tag = 1 (node numbers)
-        """
-
-        try:
-            from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
-
-            [SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
-            
-            root = Tk()
-            frame = Frame(root)
-            frame.pack()
-            button = Button(frame, text="OK", fg="red", command=frame.quit)
-            button.pack(side=BOTTOM)
-            canvas = Canvas(frame,bg="white", width=WIDTH, height=HEIGHT)
-            canvas.pack()
-            text = Label(frame, width=20, height=10, text='triangle mesh')
-            text.pack()
-
-            #print self.meshTriangles
-            for triangle in self.meshTriangles:
-                triangle.draw(canvas,1,
-                              scale = SCALE,
-                              xoffset = xoffset,
-                              yoffset = yoffset )
-                
-            root.mainloop()
-
-        except:
-            print "Unexpected error:", sys.exc_info()[0]
-            raise
-
-            #print """
-            #node::plot Failed.
-            #Most probably, the Tkinter module is not available.
-            #"""
-
-    def plotUserSegments(self,tag = 0,WIDTH = 400,HEIGHT = 400):
-        """
-        Plots all node connections.
-        tag = 0 (no node numbers), tag = 1 (node numbers)
-        """
-
-        try:
-            from Tkinter import Tk, Frame, Button, Canvas, BOTTOM, Label
-            
-            [SCALE, xoffset, yoffset] = self.scaleoffset( WIDTH, HEIGHT)
-
-            root = Tk()
-            frame = Frame(root)
-            frame.pack()
-            button = Button(frame, text="OK", fg="red", command=frame.quit)
-            button.pack(side=BOTTOM)
-            canvas = Canvas(frame, bg="white", width=WIDTH, height=HEIGHT)
-            canvas.pack()
-            text = Label(frame, width=20, height=10, text='user segments')
-            text.pack()
-           
-            for segment in self.getUserSegments():
-                segment.draw(canvas,SCALE, xoffset, yoffset )
-                
-            root.mainloop()
-
-        except:
-            print "Unexpected error:", sys.exc_info()[0]
-            raise
-
-            #print """
-            #node::plot Failed.
-            #Most probably, the Tkinter module is not available.
-            #"""
-            
-
+
+         
     def exportASCIIobj(self,ofile):
         """

inundation/pmesh/test_mesh.py

@@ -1198 +1198 @@
                         'Wrong segment tag.')
         
+    def test_addVertsSegs2(self):
+        geo = Geo_reference(56,5,10)
+        m = Mesh(geo_reference=geo)
+        dict = {}
+        dict['points'] = [[2.0, 1.0], [3.0, 1.0], [2.0, 2.0]]
+        dict['segments'] = [[0, 1], [1, 2], [2,0]]
+        dict['segment_tags'] = ['','do-op','']
+        m.addVertsSegs(dict)
+
     def test_exportASCIImeshfile(self):
     
@@ -1490 +1499 @@
 
 
-#___________begining of Peters tests
+#___________beginning of Peters tests
 
     def test_set_stuff(self):
@@ -1644 +1653 @@
         assert not range(17) in a
 
+#___________end of Peters tests
+
+    def test_addRegionFromPolygon(self):
+        m=Mesh()
+        m.addRegionFromPolygon([[0,0],[1,0],[0,1]])
+        self.failUnless(len(m.regions)==1,
+                        'FAILED!')
+        self.failUnless(len(m.getUserSegments())==3,
+                        'FAILED!')
+        self.failUnless(len(m.userVertices)==3,
+                        'FAILED!')
+       
+    def test_addRegionFromPolygon2(self):
+        m=Mesh()
+        m.addRegionFromPolygon([[0,0],[1,0],[1,1],[0,1]],
+                               {'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!') 
+        
+        
 
 def list_comp(A,B):

inundation/pyvolution/coordinate_transforms/geo_reference.py

@@ -115 +115 @@
         """
         Change the reference of a list or Numeric array of points.
-        If the points do not have ageo ref, assume 'absolute' values
+        If the points do not have a geo ref, assume 'absolute' values
         """
         is_list = False

inundation/utilities/polygon.py

@@ -461 +461 @@
     return points
 
-
+def point_in_polygon(polygon, delta=1e-8):
+    """Return a point inside a given polygon which will be close to the
+    polygon edge.
+
+    Input:
+       polygon - list of vertices of polygon
+       delta - the square root of 2 * delta is the maximum distance from the
+       polygon points and the returned point.
+    Output:
+       points - a point inside polygon
+
+       searches in all diagonals and up and down (not left and right)
+    """
+    import exceptions
+    class Found(exceptions.Exception): pass
+
+    point_in = False
+    while not point_in:
+        try:
+            for poly_point in polygon: #[1:]:
+                for x_mult in range (-1,2):
+                    for y_mult in range (-1,2):
+                        x = poly_point[0]
+                        y = poly_point[1]
+                        if x == 0:
+                            x_delta = x_mult*delta
+                        else:
+                            x_delta = x+x_mult*x*delta
+                            
+                        if y == 0:
+                            y_delta = y_mult*delta
+                        else:
+                            y_delta = y+y_mult*y*delta
+                            
+                        point = [x_delta, y_delta]
+                        #print "point",point
+                        if inside_polygon( point, polygon, closed = False ):
+                            raise Found
+        except Found:  
+            point_in = True
+        else:
+            delta = delta*0.1  
+    return point
 
 

inundation/utilities/test_polygon.py

@@ -379 +379 @@
             assert not inside_polygon(point, ex_poly), '%s' %str(point)                        
 
-
+    def test_point_in_polygon(self):
+        polygon = [[0,0], [1,0], [1,1], [0,1]]
+        point = point_in_polygon(polygon)
+        assert inside_polygon(point, polygon)
+
+        #this get's into a vicious loop
+        polygon = [[1e32,1e54], [1,0], [1,1], [0,1]]
+        point = point_in_polygon(polygon)
+        assert inside_polygon(point, polygon)
+
+
+        polygon = [[0,0], [1,0], [1,1], [1e8,1e8]]
+        point = point_in_polygon(polygon)
+        assert inside_polygon(point, polygon)
+        
+        polygon = [[1e32,1e54], [-1e32,1e54],[1e32,-1e54]]
+        point = point_in_polygon(polygon)
+        assert inside_polygon(point, polygon)
+        
+        polygon = [[1e32,1e54], [1,0], [0,1]]
+        point = point_in_polygon(polygon)
+        assert inside_polygon(point, polygon)
 
 #-------------------------------------------------------------