"""Create mesh for lwru2 validation of the Oshiri island tsunami """ #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # Assume that the root AnuGA dir is included in your pythonpath #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! from pyvolution.coordinate_transforms.geo_reference import Geo_reference from pyvolution.coordinate_transforms.redfearn import redfearn def convert_points_from_latlon_to_utm(polygon, refzone): points = [] for point in polygon: zone, easting, northing = redfearn(point[0], point[1]) #FIXME: Use point.latitude etc assert zone == refzone points.append([easting, northing]) return points def create_region(polygon, tags, refzone): """Create dictionary representation of region bounded by given polygon for use with pmesh """ # dict = {} 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] ) dict['segments'] = segments #Create tags #E.g. ['wall', 'wall', 'ocean', 'wall'] if tags is not None: segment_tags = [0]*N for key in tags: indices = tags[key] for i in indices: segment_tags[i] = key dict['segment_tags'] = segment_tags return dict def create_mesh(bounding_polygon, boundary_tags, resolution, filename = None, interior_regions = None): """Create mesh from bounding polygon, tags for all segments and resolution. Polygon is a list of points in latitude and longitudes - decimal degrees Tags is a list of symbolic tags Resolution is the maximal area Interior_regions is a list of tuples consisting of (polygon, resolution) This function does not allow segments to share points - use underlying pmesh functionality for that FIXME - FUTURE: Use class Point_set Take multiple polygons Accept deg, min, sec, e.g. [ [(-20,30,55), (116, 20, 00)], ...] FIXME: This whole thing needs to be completely rethought """ from pmesh.mesh import Mesh from pyvolution.coordinate_transforms.redfearn import redfearn from pyvolution.util import populate_polygon #Make georef #FIXME: Pass in geo or create automatically somehow import project refzone = project.refzone #FIXME mesh_origin = project.mesh_origin geo = Geo_reference(xllcorner = mesh_origin[1], #From dem yllcorner = mesh_origin[2], zone = refzone) print "***********************" print "geo ref", geo print "***********************" m = Mesh(geo_reference=geo) #Convert to UTM bounding_polygon = convert_points_from_latlon_to_utm(bounding_polygon, refzone) #Do bounding polygon region_dict = create_region(bounding_polygon, boundary_tags, refzone) m.addVertsSegs(region_dict) #Find one point inside region automatically if interior_regions is not None: excluded_polygons = [] for P, res in interior_regions: polygon = convert_points_from_latlon_to_utm(P, refzone) excluded_polygons.append( polygon ) else: excluded_polygons = None from Numeric import array [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons) print 'I:', inner_point, resolution inner = m.addRegionEN(inner_point[0], inner_point[1]) inner.setMaxArea(resolution) #Do interior regions if interior_regions is not None: for P, res in interior_regions: polygon = convert_points_from_latlon_to_utm(P, refzone) region_dict = create_region(polygon, None, refzone) m.addVertsSegs(region_dict) [inner_point] = populate_polygon(polygon, 1) print 'I', inner_point, res X = m.addRegionEN(inner_point[0], inner_point[1]) X.setMaxArea(res) m.generateMesh('pzq28.0za1000000a') m.export_mesh_file(filename)