Changeset 2276 for inundation/pmesh
- Timestamp:
- Jan 25, 2006, 12:11:06 PM (19 years ago)
- Location:
- inundation/pmesh
- Files:
-
- 2 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pmesh/create_gong.py
r1086 r2276 33 33 m = Mesh(geo_reference=geo) 34 34 35 dict = {} 36 dict['points'] = [[west, north_beach_northing - border], #sw 35 factor = 1000 #low res 10000, high res 1000 36 37 poly = [[west, north_beach_northing - border], #sw 37 38 [west, bellambi_northing + border], #nw 38 39 [sea_boundary, … … 43 44 north_beach_northing - border] #s_sea_boundary 44 45 ] 45 dict['segments'] = [[0,1],[1,2],[2,3], 46 [3,4],[4,5],[5,0], # the outer boarder 47 [2,5] # the sea_boundary line 48 ] 46 tags = {'wall':[0,1,2,3,4,5]} 47 48 m.add_region_from_polygon(poly, tags, max_triangle_area=factor*1000) 49 poly = [ 50 [inner_west_south,north_beach_northing], #sw 51 [inner_west_north,bellambi_northing], #nw 52 [inner_east,bellambi_northing], #ne 53 [inner_east,north_beach_northing] #se 54 ] 55 m.add_region_from_polygon(poly, max_triangle_area=factor*100) 49 56 50 m.addVertsSegs(dict)51 52 dict = {}53 dict['points'] = [54 [inner_west_south,north_beach_northing], #sw55 [inner_west_north,bellambi_northing], #nw56 [inner_east,bellambi_northing], #ne57 [inner_east,north_beach_northing] #se58 ]59 dict['segments'] = [[0,1],[1,2],[2,3],[3,0]] # the inner boarder60 m.addVertsSegs(dict)61 57 62 58 #Even more inner points … … 68 64 north_beach_northing = 6183000.0 69 65 70 dict = {} 71 dict['points'] = [ 66 poly = [ 72 67 [flagstaff_x,flagstaff_y], #n 73 68 [west_x,west_y], #w 74 69 [north_beach_easting,north_beach_northing], #s 75 70 ] 76 dict['segments'] = [[0,1],[1,2],[2,0]] # the inner boarder 77 m.addVertsSegs(dict) 71 m.add_region_from_polygon(poly, max_triangle_area=factor*10) 78 72 79 factor = 1000 #low res 10000, high res 100080 low = m.addRegionEN(east-0.5, bellambi_northing + border-0.5)81 low.setMaxArea(100*factor)82 83 medium = m.addRegionEN(sea_boundary-0.5, bellambi_northing + border-0.5)84 medium.setMaxArea(10*factor)85 86 high = m.addRegionEN(inner_east-0.5,bellambi_northing-0.5)87 high.setMaxArea(5*factor) #doing 0.5, fac = 1000 gives degenerate error88 89 center = m.addRegionEN(flagstaff_x,flagstaff_y-0.5)90 center.setMaxArea(0.05*factor) #doing 0.5, fac = 1000 gives degenerate error91 73 #m.generateMesh() 92 74 m.generateMesh("pzq28.0za1000000a") -
inundation/pmesh/mesh.py
r2261 r2276 604 604 605 605 def add_region(self, x,y, geo_reference=None): 606 [x,y] = self.geo_reference.change_points_geo_ref([x,y], 606 """ 607 adds a point, which represents a region. 608 609 The point data can have it's own geo_refernece. 610 If geo_reference is None the data is asumed to be absolute 611 """ 612 [[x,y]] = self.geo_reference.change_points_geo_ref([x,y], 607 613 points_geo_ref=geo_reference) 608 614 return self._addRegion(x, y) … … 614 620 615 621 def add_region_from_polygon(self, polygon, tags=None, 616 max_area=None, geo_reference=None):622 max_triangle_area=None, geo_reference=None): 617 623 """ 618 624 Add a polygon with tags to the current mesh, as a region. 619 625 The maxArea of the region can be specified. 620 626 621 If a geo_ref is given, this is used.627 If a geo_reference of the polygon points is given, this is used. 622 628 If not; 623 629 The x,y info is assumed to be Easting and Northing, absolute, … … 632 638 633 639 """ 634 #FIXME: take into account georef on mesh side 635 polygon = self.geo_reference.change_points_geo_ref(polygon, 636 points_geo_ref=geo_reference) 637 #print "polygon - should be relative to mesh geo_ref",polygon 640 #take into account georef on mesh side 641 #polygon = self.geo_reference.change_points_geo_ref(polygon, 642 # points_geo_ref=geo_reference) 643 #print "polygon - should be relative to poly geo_ref",polygon 644 645 #get absolute values 646 if geo_reference is not None: 647 polygon = geo_reference.get_absolute(polygon) 648 # polygon is now absolute 649 #print "polygon should be absolute",polygon 638 650 #create points, segs and tags 639 651 region_dict = {} … … 666 678 667 679 self.addVertsSegs(region_dict) #this is assuming absolute geos 668 669 #get inner point 670 inner_point = point_in_polygon(polygon) 671 inner = self.add_region(inner_point[0], inner_point[1], 672 geo_reference=self.geo_reference) 673 674 if max_area is not None: 675 inner.setMaxArea(max_area) 680 681 inner = None 682 683 if max_triangle_area is not None: 684 #get inner point - absolute values 685 inner_point = point_in_polygon(polygon) 686 inner = self.add_region(inner_point[0], inner_point[1], 687 geo_reference=None) 688 inner.setMaxArea(max_triangle_area) 689 676 690 return inner 677 691 … … 728 742 return s 729 743 730 def generateMesh(self, mode = None, maxArea = None, 744 def generate_mesh(self, 745 maximum_triangle_area=None, 746 minimum_triangle_angle=None, 747 verbose=False): 748 if verbose is True: 749 silent = '' 750 else: 751 silent = 'Q' 752 self.generateMesh(mode = silent +"pzq"+str(minimum_triangle_angle) 753 +"a"+str(maximum_triangle_area) 754 +"a") 755 #The last a is so areas for regions will be used 756 757 def generateMesh(self, mode = None, maxArea = None, minAngle= None, 731 758 isRegionalMaxAreas = True): 732 759 """ … … 743 770 self.mode = mode 744 771 745 if not re.match('p',self.mode):772 if not self.mode.find('p'): 746 773 self.mode += 'p' #p - read a planar straight line graph. 747 774 #there must be segments to use this switch 748 775 # TODO throw an aception if there aren't seg's 749 776 # it's more comlex than this. eg holes 750 if not re.match('z',self.mode):777 if not self.mode.find('z'): 751 778 self.mode += 'z' # z - Number all items starting from zero 752 779 # (rather than one) 753 if not re.match('n',self.mode):780 if self.mode.find('n'): 754 781 self.mode += 'n' # n - output a list of neighboring triangles 755 if not re.match('A',self.mode):782 if not self.mode.find('A'): 756 783 self.mode += 'A' # A - output region attribute list for triangles 757 if not re.match('V',self.mode) and not re.match('Q',self.mode): 784 785 if not self.mode.find('V') and not self.mode.find('Q'): 758 786 self.mode += 'V' # V - output info about what Triangle is doing 759 787 788 if not self.mode.find('q') and minAngle is not None: 789 # print "**********8minAngle******** ",minAngle 790 min_angle = 'q' + str(minAngle) 791 self.mode += min_angle # z - Number all items starting from zero 792 # (rather than one) 760 793 if maxArea != None: 761 794 self.mode += 'a' + str(maxArea) 762 795 #FIXME (DSG-DSG) This isn't explained. I don't think it works 796 # well with maxArea = None as well. 763 797 if isRegionalMaxAreas: 764 798 self.mode += 'a' … … 1255 1289 1256 1290 Assume the values are in Eastings and Northings, with no reference 1257 point 1291 point. eg absolute 1258 1292 """ 1259 1293 if not outlineDict.has_key('segment_tags'): -
inundation/pmesh/test_mesh.py
r2261 r2276 6 6 from load_mesh.loadASCII import * 7 7 from coordinate_transforms.geo_reference import Geo_reference 8 from utilities.polygon import inside_polygon 8 9 9 10 class meshTestCase(unittest.TestCase): 10 11 def setUp(self): 11 12 pass 12 13 14 13 15 14 def tearDown(self): 16 15 pass … … 79 78 80 79 81 80 # FIXME add test for minAngle 82 81 def testgenerateMesh(self): 83 82 a = Vertex (0.0, 0.0) … … 1657 1656 def test_add_region_from_polygon(self): 1658 1657 m=Mesh() 1659 m.add_region_from_polygon([[0,0],[1,0],[0,1]]) 1658 region = m.add_region_from_polygon([[0,0],[1,0],[0,1]], 1659 max_triangle_area = 88) 1660 1660 self.failUnless(len(m.regions)==1, 1661 'FAILED!') 1662 self.failUnless(region.getMaxArea()==88, 1661 1663 'FAILED!') 1662 1664 self.failUnless(len(m.getUserSegments())==3, … … 1668 1670 m=Mesh() 1669 1671 m.add_region_from_polygon([[0,0],[1,0],[1,1],[0,1]], 1670 {'tagin':[0,1],'bom':[2]}) 1672 {'tagin':[0,1],'bom':[2]}, 1673 max_triangle_area=10) 1671 1674 self.failUnless(len(m.regions)==1, 1672 1675 'FAILED!') … … 1687 1690 1688 1691 def test_add_region_from_polygon3(self): 1689 x= 01690 y= 01692 x=-500 1693 y=-1000 1691 1694 m=Mesh(geo_reference=Geo_reference(56,x,y)) 1692 polygon = [[0,0],[1,0],[1,1],[0,1]] 1693 x_p = 1000 1694 y_p = 4000 1695 1696 # These are the absolute values 1697 polygon_absolute = [[0,0],[1,0],[1,1],[0,1]] 1698 1699 x_p = -10 1700 y_p = -40 1695 1701 geo_ref_poly = Geo_reference(56, x_p, y_p) 1696 polygon = geo_ref_poly.change_points_geo_ref(polygon) 1697 #print "polygon", polygon 1698 m.add_region_from_polygon(polygon, 1699 {'tagin':[0,1],'bom':[2]}, 1700 geo_reference=geo_ref_poly) 1702 polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute) 1703 1704 poly_point = m.add_region_from_polygon(polygon, 1705 {'tagin':[0,1],'bom':[2]}, 1706 geo_reference=geo_ref_poly, 1707 max_triangle_area=10) 1708 # poly_point values are relative to the mesh geo-ref 1709 # make them absolute 1710 self.failUnless(inside_polygon([poly_point.x+x,poly_point.y+y], 1711 polygon_absolute, closed = False), 1712 'FAILED!') 1713 1701 1714 self.failUnless(len(m.regions)==1, 1702 1715 'FAILED!') … … 1723 1736 new_point_x = new_point.x + m.geo_reference.get_xllcorner() 1724 1737 point_y = point[1] + geo_ref_poly.get_yllcorner() 1738 #print "new_point.y",new_point.y 1739 #print "m.geo_ref.get_yllcorner()",m.geo_reference.get_yllcorner() 1725 1740 new_point_y = new_point.y + m.geo_reference.get_yllcorner() 1726 1741 #print "point_y",point_y … … 1731 1746 1732 1747 1733 #REJIG1734 1748 def test_add_region_from_polygon4(self): 1735 x= 01736 y= 01749 x=50000 1750 y=1000 1737 1751 m=Mesh(geo_reference=Geo_reference(56,x,y)) 1738 1752 polygon = [[0,0],[1,0],[1,1],[0,1]] 1739 1753 1740 1754 m.add_region_from_polygon(polygon, 1741 {'tagin':[0,1],'bom':[2]}) 1755 {'tagin':[0,1],'bom':[2]}, 1756 max_triangle_area=10) 1742 1757 self.failUnless(len(m.regions)==1, 1743 1758 'FAILED!')
Note: See TracChangeset
for help on using the changeset viewer.