[2276] | 1 | |
---|
| 2 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 3 | # Assume that the root AnuGA dir (inundation) is included in your pythonpath |
---|
| 4 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 5 | |
---|
[2402] | 6 | from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE |
---|
[3182] | 7 | from utilities.polygon import point_in_polygon ,populate_polygon |
---|
[2402] | 8 | from utilities.numerical_tools import ensure_numeric |
---|
| 9 | from Numeric import Float |
---|
[3051] | 10 | from utilities.polygon import inside_polygon |
---|
[2397] | 11 | |
---|
| 12 | # This is due to pmesh being a package and a module and |
---|
[2282] | 13 | # the current dir being unknown |
---|
| 14 | try: |
---|
| 15 | from pmesh.mesh import Mesh |
---|
| 16 | except ImportError: |
---|
| 17 | from mesh import Mesh |
---|
[2402] | 18 | |
---|
[3056] | 19 | import exceptions |
---|
| 20 | class PolygonError(exceptions.Exception): pass |
---|
[2402] | 21 | |
---|
[2276] | 22 | def create_mesh_from_regions(bounding_polygon, |
---|
| 23 | boundary_tags, |
---|
[3193] | 24 | maximum_triangle_area=None, |
---|
[2276] | 25 | filename=None, |
---|
| 26 | interior_regions=None, |
---|
| 27 | poly_geo_reference=None, |
---|
| 28 | mesh_geo_reference=None, |
---|
[3193] | 29 | minimum_triangle_angle=28.0, |
---|
| 30 | verbose=True): |
---|
[2276] | 31 | """Create mesh from bounding polygons, and resolutions. |
---|
| 32 | |
---|
[2402] | 33 | bounding_polygon is a list of points in Eastings and Northings, |
---|
| 34 | relative to the poly_geo_reference. |
---|
[2276] | 35 | |
---|
| 36 | Boundary tags is a dictionary of symbolic tags. For every tag there |
---|
[2391] | 37 | is a list of indices referring to segments associated with that tag |
---|
[2276] | 38 | |
---|
[3193] | 39 | maximum_triangle_area is the maximal area per triangle |
---|
| 40 | for the bounding polygon, excluding the interior regions. |
---|
[2276] | 41 | |
---|
| 42 | Interior_regions is a list of tuples consisting of (polygon, resolution) |
---|
| 43 | for each region to be separately refined. |
---|
| 44 | |
---|
| 45 | This function does not allow segments to share points - use underlying |
---|
| 46 | pmesh functionality for that |
---|
| 47 | |
---|
[3193] | 48 | poly_geo_reference is the geo_reference of the bounding polygon and |
---|
| 49 | the interior polygons. |
---|
[2402] | 50 | If none, assume absolute. Please pass one though, since absolute |
---|
| 51 | references have a zone. |
---|
[2276] | 52 | |
---|
| 53 | mesh_geo_reference is the geo_reference of the mesh to be created. |
---|
[2402] | 54 | If none is given one will be automatically generated. It was use |
---|
| 55 | the lower left hand corner of bounding_polygon (absolute) |
---|
| 56 | as the x and y values for the geo_ref. |
---|
[2276] | 57 | |
---|
| 58 | Returns the mesh instance if no finename is given |
---|
[2430] | 59 | |
---|
| 60 | Note, interior regions should be fully nested, as overlaps may cause |
---|
| 61 | unintended resolutions. |
---|
[2276] | 62 | |
---|
| 63 | """ |
---|
[2402] | 64 | #FIXME (OLE-DSG) |
---|
[2276] | 65 | # check the segment indexes - throw an error if they are out of bounds |
---|
[3056] | 66 | #(DSG) Yes! |
---|
| 67 | |
---|
[2391] | 68 | #In addition I reckon the polygons could be of class Geospatial_data |
---|
[3056] | 69 | #(DSG) Yes! |
---|
[2391] | 70 | |
---|
[3038] | 71 | # First check that interior polygons are fully contained in bounding polygon |
---|
[3056] | 72 | #Note, Both poly's have the same geo_ref, therefore don't take into account |
---|
| 73 | # geo_ref |
---|
| 74 | if interior_regions is not None: |
---|
[3183] | 75 | # Test that all the interior polygons are inside the bounding_poly |
---|
[3056] | 76 | for interior_polygon, res in interior_regions: |
---|
| 77 | indices = inside_polygon(interior_polygon, bounding_polygon, |
---|
| 78 | closed = True, verbose = False) |
---|
[3051] | 79 | |
---|
[3056] | 80 | if len(indices) <> len(interior_polygon): |
---|
| 81 | msg = 'Interior polygon %s is outside bounding polygon: %s'\ |
---|
| 82 | %(str(interior_polygon), str(bounding_polygon)) |
---|
| 83 | raise PolygonError, msg |
---|
[3038] | 84 | |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | # Resolve geo referencing |
---|
[2402] | 88 | if mesh_geo_reference is None: |
---|
| 89 | bounding_polygon = ensure_numeric(bounding_polygon, Float) |
---|
| 90 | xllcorner = min(bounding_polygon[:,0]) |
---|
| 91 | yllcorner = min(bounding_polygon[:,1]) |
---|
| 92 | # |
---|
| 93 | if poly_geo_reference is None: |
---|
| 94 | zone = DEFAULT_ZONE |
---|
| 95 | else: |
---|
| 96 | zone = poly_geo_reference.get_zone() |
---|
| 97 | [(xllcorner,yllcorner)] = poly_geo_reference.get_absolute( \ |
---|
| 98 | [(xllcorner,yllcorner)]) |
---|
| 99 | # create a geo_ref, based on the llc of the bounding_polygon |
---|
| 100 | mesh_geo_reference = Geo_reference(xllcorner = xllcorner, |
---|
[2408] | 101 | yllcorner = yllcorner, |
---|
| 102 | zone = zone) |
---|
[2402] | 103 | |
---|
[2276] | 104 | m = Mesh(geo_reference=mesh_geo_reference) |
---|
| 105 | |
---|
[3038] | 106 | # Do bounding polygon |
---|
[2276] | 107 | m.add_region_from_polygon(bounding_polygon, |
---|
[2391] | 108 | tags=boundary_tags, |
---|
| 109 | geo_reference=poly_geo_reference) |
---|
[2276] | 110 | |
---|
[3038] | 111 | # Find one point inside region automatically |
---|
[2276] | 112 | if interior_regions is not None: |
---|
| 113 | excluded_polygons = [] |
---|
| 114 | for polygon, res in interior_regions: |
---|
| 115 | excluded_polygons.append( polygon ) |
---|
| 116 | else: |
---|
| 117 | excluded_polygons = None |
---|
| 118 | |
---|
[3038] | 119 | # Convert bounding poly to absolute values |
---|
[2276] | 120 | # this sort of thing can be fixed with the geo_points class |
---|
| 121 | if poly_geo_reference is not None: |
---|
| 122 | bounding_polygon_absolute = \ |
---|
| 123 | poly_geo_reference.get_absolute(bounding_polygon) |
---|
| 124 | else: |
---|
| 125 | bounding_polygon_absolute = bounding_polygon |
---|
[3182] | 126 | |
---|
| 127 | inner_point = point_in_polygon(bounding_polygon_absolute) |
---|
[2276] | 128 | inner = m.add_region(inner_point[0], inner_point[1]) |
---|
| 129 | inner.setMaxArea(maximum_triangle_area) |
---|
| 130 | |
---|
[3038] | 131 | # Do interior regions |
---|
[2276] | 132 | if interior_regions is not None: |
---|
| 133 | for polygon, res in interior_regions: |
---|
| 134 | m.add_region_from_polygon(polygon, |
---|
| 135 | geo_reference=poly_geo_reference) |
---|
| 136 | # convert bounding poly to absolute values |
---|
| 137 | if poly_geo_reference is not None: |
---|
| 138 | polygon_absolute = \ |
---|
| 139 | poly_geo_reference.get_absolute(polygon) |
---|
| 140 | else: |
---|
| 141 | polygon_absolute = polygon |
---|
[3183] | 142 | inner_point = point_in_polygon(polygon_absolute) |
---|
[2276] | 143 | region = m.add_region(inner_point[0], inner_point[1]) |
---|
| 144 | region.setMaxArea(res) |
---|
[2391] | 145 | |
---|
[2276] | 146 | |
---|
| 147 | if filename is None: |
---|
| 148 | return m |
---|
| 149 | else: |
---|
[3042] | 150 | m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle, |
---|
[3193] | 151 | verbose=verbose) |
---|
[2276] | 152 | m.export_mesh_file(filename) |
---|
[2402] | 153 | |
---|
| 154 | |
---|