[1784] | 1 | """Create mesh for lwru2 validation of the Oshiri island tsunami |
---|
| 2 | """ |
---|
| 3 | |
---|
| 4 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 5 | # Assume that the root AnuGA dir is included in your pythonpath |
---|
| 6 | #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 7 | |
---|
[1832] | 8 | |
---|
[1784] | 9 | from pyvolution.coordinate_transforms.geo_reference import Geo_reference |
---|
[1796] | 10 | from pyvolution.coordinate_transforms.redfearn import redfearn |
---|
[1784] | 11 | |
---|
[1796] | 12 | |
---|
[1855] | 13 | def convert_points_from_latlon_to_utm(polygon, refzone): |
---|
| 14 | points = [] |
---|
| 15 | for point in polygon: |
---|
| 16 | zone, easting, northing = redfearn(point[0], point[1]) #FIXME: Use point.latitude etc |
---|
| 17 | assert zone == refzone |
---|
| 18 | points.append([easting, northing]) |
---|
[1796] | 19 | |
---|
[1855] | 20 | return points |
---|
| 21 | |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | def create_region(polygon, tags, refzone): |
---|
| 27 | """Create dictionary representation of region bounded by given polygon for use with pmesh |
---|
| 28 | """ |
---|
| 29 | |
---|
| 30 | |
---|
| 31 | # |
---|
| 32 | dict = {} |
---|
| 33 | dict['points'] = polygon |
---|
| 34 | |
---|
| 35 | #Create segments |
---|
| 36 | #E.g. [[0,1], [1,2], [2,3], [3,0]] |
---|
| 37 | segments = [] |
---|
| 38 | N = len(polygon) |
---|
| 39 | for i in range(N): |
---|
| 40 | lo = i |
---|
| 41 | hi = (lo + 1) % N |
---|
| 42 | segments.append( [lo, hi] ) |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | dict['segments'] = segments |
---|
| 46 | |
---|
| 47 | |
---|
| 48 | #Create tags |
---|
| 49 | #E.g. ['wall', 'wall', 'ocean', 'wall'] |
---|
| 50 | |
---|
| 51 | if tags is not None: |
---|
| 52 | segment_tags = [0]*N |
---|
| 53 | for key in tags: |
---|
| 54 | indices = tags[key] |
---|
| 55 | for i in indices: |
---|
| 56 | segment_tags[i] = key |
---|
| 57 | |
---|
| 58 | dict['segment_tags'] = segment_tags |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | return dict |
---|
| 63 | |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | def create_mesh(bounding_polygon, boundary_tags, resolution, |
---|
| 67 | filename = None, interior_regions = None): |
---|
[1832] | 68 | """Create mesh from bounding polygon, tags for all segments and resolution. |
---|
[1784] | 69 | |
---|
[1832] | 70 | Polygon is a list of points in latitude and longitudes - decimal degrees |
---|
[1784] | 71 | |
---|
[1832] | 72 | |
---|
| 73 | Tags is a list of symbolic tags |
---|
| 74 | |
---|
| 75 | Resolution is the maximal area |
---|
[1784] | 76 | |
---|
[1855] | 77 | Interior_regions is a list of tuples consisting of (polygon, resolution) |
---|
| 78 | |
---|
| 79 | This function does not allow segments to share points - use underlying |
---|
| 80 | pmesh functionality for that |
---|
| 81 | |
---|
[1832] | 82 | FIXME - FUTURE: |
---|
[1837] | 83 | Use class Point_set |
---|
[1832] | 84 | Take multiple polygons |
---|
[1784] | 85 | |
---|
[1832] | 86 | Accept deg, min, sec, e.g. |
---|
| 87 | [ [(-20,30,55), (116, 20, 00)], ...] |
---|
[1855] | 88 | |
---|
| 89 | |
---|
| 90 | FIXME: This whole thing needs to be completely rethought |
---|
[1832] | 91 | |
---|
| 92 | """ |
---|
[1784] | 93 | |
---|
[1832] | 94 | from pmesh.mesh import * |
---|
| 95 | from pyvolution.coordinate_transforms.redfearn import redfearn |
---|
[1855] | 96 | from pyvolution.util import populate_polygon |
---|
[1832] | 97 | |
---|
[1796] | 98 | |
---|
[1855] | 99 | #Make georef |
---|
| 100 | #FIXME: Pass in geo or create automatically somehow |
---|
[1832] | 101 | import project |
---|
| 102 | refzone = project.refzone #FIXME |
---|
| 103 | mesh_origin = project.mesh_origin |
---|
| 104 | geo = Geo_reference(xllcorner = mesh_origin[1], #From dem |
---|
| 105 | yllcorner = mesh_origin[2], |
---|
[1784] | 106 | zone = refzone) |
---|
[1818] | 107 | |
---|
[1832] | 108 | |
---|
[1855] | 109 | print "***********************" |
---|
| 110 | print "geo ref", geo |
---|
| 111 | print "***********************" |
---|
| 112 | |
---|
| 113 | m = Mesh(geo_reference=geo) |
---|
[1832] | 114 | |
---|
| 115 | |
---|
[1855] | 116 | |
---|
| 117 | #Convert to UTM |
---|
| 118 | bounding_polygon = convert_points_from_latlon_to_utm(bounding_polygon, refzone) |
---|
[1784] | 119 | |
---|
[1855] | 120 | |
---|
| 121 | #Do bounding polygon |
---|
| 122 | region_dict = create_region(bounding_polygon, boundary_tags, refzone) |
---|
| 123 | m.addVertsSegs(region_dict) |
---|
[1832] | 124 | |
---|
| 125 | |
---|
[1855] | 126 | #Find one point inside region automatically |
---|
| 127 | if interior_regions is not None: |
---|
| 128 | excluded_polygons = [] |
---|
| 129 | for P, res in interior_regions: |
---|
| 130 | polygon = convert_points_from_latlon_to_utm(P, refzone) |
---|
| 131 | excluded_polygons.append( polygon ) |
---|
| 132 | else: |
---|
| 133 | excluded_polygons = None |
---|
[1832] | 134 | |
---|
[1855] | 135 | from Numeric import array |
---|
[1784] | 136 | |
---|
[1855] | 137 | [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons) |
---|
| 138 | print 'I:', inner_point, resolution |
---|
| 139 | inner = m.addRegionEN(inner_point[0], inner_point[1]) |
---|
| 140 | inner.setMaxArea(resolution) |
---|
[1784] | 141 | |
---|
| 142 | |
---|
[1855] | 143 | |
---|
| 144 | #Do interior regions |
---|
| 145 | if interior_regions is not None: |
---|
| 146 | for P, res in interior_regions: |
---|
| 147 | polygon = convert_points_from_latlon_to_utm(P, refzone) |
---|
| 148 | region_dict = create_region(polygon, None, refzone) |
---|
| 149 | |
---|
| 150 | m.addVertsSegs(region_dict) |
---|
| 151 | |
---|
| 152 | [inner_point] = populate_polygon(polygon, 1) |
---|
| 153 | print 'I', inner_point, res |
---|
| 154 | X = m.addRegionEN(inner_point[0], inner_point[1]) |
---|
| 155 | X.setMaxArea(res) |
---|
| 156 | |
---|
| 157 | |
---|
[1784] | 158 | |
---|
[1832] | 159 | m.generateMesh('pzq28.0za1000000a') |
---|
| 160 | |
---|
| 161 | m.export_mesh_file(filename) |
---|
| 162 | |
---|
[1784] | 163 | |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | |
---|
| 167 | |
---|