source: production/karratha_2005/create_mesh.py @ 1920

Last change on this file since 1920 was 1868, checked in by ole, 18 years ago

Ongoing work with karratha

File size: 4.2 KB
Line 
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
8
9from pyvolution.coordinate_transforms.geo_reference import Geo_reference
10from pyvolution.coordinate_transforms.redfearn import redfearn
11
12
13def 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])
19
20    return points   
21
22
23
24
25
26def 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
66def create_mesh(bounding_polygon, boundary_tags, resolution,
67                filename = None, interior_regions = None):
68    """Create mesh from bounding polygon, tags for all segments and resolution.
69
70    Polygon is a list of points in latitude and longitudes - decimal degrees
71
72
73    Tags is a list of symbolic tags
74
75    Resolution is the maximal area
76
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
82    FIXME - FUTURE:
83      Use class Point_set
84      Take multiple polygons
85
86      Accept deg, min, sec, e.g.
87      [ [(-20,30,55), (116, 20, 00)], ...]
88
89
90    FIXME: This whole thing needs to be completely rethought
91   
92    """
93
94    from pmesh.mesh import Mesh
95    from pyvolution.coordinate_transforms.redfearn import redfearn
96    from pyvolution.util import populate_polygon
97   
98
99    #Make georef
100    #FIXME: Pass in geo or create automatically somehow
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],
106                        zone = refzone)
107
108
109
110   
111    print "***********************"
112    print "geo ref", geo
113    print "***********************"
114   
115    m = Mesh(geo_reference=geo)
116
117
118
119    #Convert to UTM
120    bounding_polygon = convert_points_from_latlon_to_utm(bounding_polygon, refzone)   
121   
122   
123    #Do bounding polygon
124    region_dict = create_region(bounding_polygon, boundary_tags, refzone)   
125    m.addVertsSegs(region_dict)
126
127
128    #Find one point inside region automatically
129    if interior_regions is not None:
130        excluded_polygons = []       
131        for P, res in interior_regions:
132            polygon = convert_points_from_latlon_to_utm(P, refzone)           
133            excluded_polygons.append( polygon )
134    else:
135        excluded_polygons = None
136
137    from Numeric import array
138   
139    [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons)
140    print 'I:', inner_point, resolution
141    inner = m.addRegionEN(inner_point[0], inner_point[1])
142    inner.setMaxArea(resolution)
143
144
145
146    #Do interior regions
147    if interior_regions is not None:   
148        for P, res in interior_regions:
149            polygon = convert_points_from_latlon_to_utm(P, refzone)                       
150            region_dict = create_region(polygon, None, refzone)
151
152            m.addVertsSegs(region_dict)
153
154            [inner_point] = populate_polygon(polygon, 1)
155            print 'I', inner_point, res
156            X = m.addRegionEN(inner_point[0], inner_point[1])
157            X.setMaxArea(res)
158           
159
160   
161    m.generateMesh('pzq28.0za1000000a')
162   
163    m.export_mesh_file(filename)
164   
165
166
167
168
169   
Note: See TracBrowser for help on using the repository browser.