| 1 | |
|---|
| 2 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 3 | # Assume that the root AnuGA dir (inundation) is included in your pythonpath |
|---|
| 4 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 5 | |
|---|
| 6 | from coordinate_transforms.geo_reference import Geo_reference |
|---|
| 7 | from utilities.polygon import populate_polygon |
|---|
| 8 | |
|---|
| 9 | # This is due to pmesh being a package and a module and |
|---|
| 10 | # the current dir being unknown |
|---|
| 11 | try: |
|---|
| 12 | from pmesh.mesh import Mesh |
|---|
| 13 | except ImportError: |
|---|
| 14 | from mesh import Mesh |
|---|
| 15 | |
|---|
| 16 | def create_mesh_from_regions(bounding_polygon, |
|---|
| 17 | boundary_tags, |
|---|
| 18 | maximum_triangle_area, |
|---|
| 19 | filename=None, |
|---|
| 20 | interior_regions=None, |
|---|
| 21 | poly_geo_reference=None, |
|---|
| 22 | mesh_geo_reference=None, |
|---|
| 23 | minimum_triangle_angle=28.0): |
|---|
| 24 | """Create mesh from bounding polygons, and resolutions. |
|---|
| 25 | |
|---|
| 26 | Polygon is a list of points in Eastings and Northings, absolute |
|---|
| 27 | |
|---|
| 28 | Boundary tags is a dictionary of symbolic tags. For every tag there |
|---|
| 29 | is a list of indices referring to segments associated with that tag |
|---|
| 30 | |
|---|
| 31 | Resolution is the maximal area per triangle for the bounding polygon |
|---|
| 32 | (excluding interior regions, see later) |
|---|
| 33 | |
|---|
| 34 | Interior_regions is a list of tuples consisting of (polygon, resolution) |
|---|
| 35 | for each region to be separately refined. |
|---|
| 36 | |
|---|
| 37 | This function does not allow segments to share points - use underlying |
|---|
| 38 | pmesh functionality for that |
|---|
| 39 | |
|---|
| 40 | poly_geo_reference is the geo_reference of the polygons. |
|---|
| 41 | If none, assume absolute. |
|---|
| 42 | |
|---|
| 43 | mesh_geo_reference is the geo_reference of the mesh to be created. |
|---|
| 44 | If none, assume absolute. |
|---|
| 45 | |
|---|
| 46 | Returns the mesh instance if no finename is given |
|---|
| 47 | |
|---|
| 48 | """ |
|---|
| 49 | # To do make maximum_triangle_area optional? |
|---|
| 50 | # check the segment indexes - throw an error if they are out of bounds |
|---|
| 51 | |
|---|
| 52 | #IDEA (Ole): Let us work out what the mesh_geo_reference is automatically: |
|---|
| 53 | #If bounding_polygon is absolute it would look as follows: |
|---|
| 54 | # (assuming bounding_polygon is a Numeric array) |
|---|
| 55 | #xllcorner = min(bounding_polygon[:,0]) |
|---|
| 56 | #yllcorner = min(bounding_polygon[:,1]) |
|---|
| 57 | # |
|---|
| 58 | #mesh_geo_reference = Geo_reference(xllcorner = xllcorner, |
|---|
| 59 | # yllcorner = yllcorner, |
|---|
| 60 | # zone = "zone from bounding polygon") |
|---|
| 61 | # |
|---|
| 62 | #In addition I reckon the polygons could be of class Geospatial_data |
|---|
| 63 | |
|---|
| 64 | m = Mesh(geo_reference=mesh_geo_reference) |
|---|
| 65 | |
|---|
| 66 | #Do bounding polygon |
|---|
| 67 | m.add_region_from_polygon(bounding_polygon, |
|---|
| 68 | tags=boundary_tags, |
|---|
| 69 | geo_reference=poly_geo_reference) |
|---|
| 70 | |
|---|
| 71 | #Find one point inside region automatically |
|---|
| 72 | if interior_regions is not None: |
|---|
| 73 | excluded_polygons = [] |
|---|
| 74 | for polygon, res in interior_regions: |
|---|
| 75 | excluded_polygons.append( polygon ) |
|---|
| 76 | else: |
|---|
| 77 | excluded_polygons = None |
|---|
| 78 | |
|---|
| 79 | from Numeric import array |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | # convert bounding poly to absolute values |
|---|
| 83 | # this sort of thing can be fixed with the geo_points class |
|---|
| 84 | if poly_geo_reference is not None: |
|---|
| 85 | bounding_polygon_absolute = \ |
|---|
| 86 | poly_geo_reference.get_absolute(bounding_polygon) |
|---|
| 87 | else: |
|---|
| 88 | bounding_polygon_absolute = bounding_polygon |
|---|
| 89 | |
|---|
| 90 | [inner_point] = populate_polygon(bounding_polygon_absolute, |
|---|
| 91 | 1, exclude = excluded_polygons) |
|---|
| 92 | inner = m.add_region(inner_point[0], inner_point[1]) |
|---|
| 93 | inner.setMaxArea(maximum_triangle_area) |
|---|
| 94 | |
|---|
| 95 | #Do interior regions |
|---|
| 96 | if interior_regions is not None: |
|---|
| 97 | for polygon, res in interior_regions: |
|---|
| 98 | #polygon = convert_points_from_latlon_to_utm(P, refzone) |
|---|
| 99 | m.add_region_from_polygon(polygon, |
|---|
| 100 | geo_reference=poly_geo_reference) |
|---|
| 101 | # convert bounding poly to absolute values |
|---|
| 102 | if poly_geo_reference is not None: |
|---|
| 103 | polygon_absolute = \ |
|---|
| 104 | poly_geo_reference.get_absolute(polygon) |
|---|
| 105 | else: |
|---|
| 106 | polygon_absolute = polygon |
|---|
| 107 | [inner_point] = populate_polygon(polygon_absolute, 1) |
|---|
| 108 | region = m.add_region(inner_point[0], inner_point[1]) |
|---|
| 109 | region.setMaxArea(res) |
|---|
| 110 | |
|---|
| 111 | m.generate_mesh(minimum_triangle_angle=minimum_triangle_angle, |
|---|
| 112 | maximum_triangle_area=maximum_triangle_area) |
|---|
| 113 | |
|---|
| 114 | if filename is None: |
|---|
| 115 | return m |
|---|
| 116 | else: |
|---|
| 117 | m.export_mesh_file(filename) |
|---|
