source: anuga_work/development/boxingday08/project.py @ 5545

Last change on this file since 5545 was 5545, checked in by ole, 16 years ago

Updated get_gauges for boxing day simulation

File size: 7.9 KB
Line 
1from anuga.utilities.numerical_tools import ensure_numeric
2from Scientific.IO.NetCDF import NetCDFFile
3from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
4    compress,zeros,fabs,allclose,ones
5from anuga.utilities.polygon import inside_polygon,read_polygon
6from os import sep
7from time import localtime, strftime, gmtime
8from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
9import os
10
11def create_sts_boundary_order_file(urs_filename,order_filename,sts_filename,verbose=False):
12    """
13    Note This is not a robust algorithm. Be Careful when applying that
14    polygon does not intersect itself and that resulting boundary actually
15    follows the path it is supposed too.
16    """
17    from anuga.shallow_water.data_manager import read_mux2_py
18    times, y, x, elevation, quantity, starttime = read_mux2_py([urs_filename],weights=ones(1,Float))
19   
20    #find indices of the southernmost gauges
21    #if more than one find gauges furthest left
22    #will fail if domain crosses 180 degrees
23    #i.e. if left most most point is not west i.e x=1 and other points 170<x<180
24    #or even in utm just a diffenet zone
25    south=min(y)
26    south_indices=compress(y-south<1e-6,arange(len(y)))
27    west=x[south_indices[0]]
28    west_index=south_indices[0]
29    if south_indices>1:
30         #find most western gauge
31         for i,south_index in enumerate(south_indices[1:]):
32              if x[south_index]<west:
33                   west_index=south_index
34   
35    n=len(x)
36    boundary=zeros((n,2),Float)
37    indices=arange(n).tolist()
38    boundary[0][0]=x[west_index]
39    boundary[0][1]=y[west_index]
40    indices.remove(west_index)
41    order=[]
42    for i in range(n-1):
43         min_dist=1e12
44         for j,k in enumerate(indices):
45              dist=sqrt((boundary[i][0]-x[k])**2+(boundary[i][1]-y[k])**2)
46              if dist<min_dist:
47                   min_dist=dist
48                   index=k
49              elif fabs(dist-min_dist)<1e-6:
50                   msg='Two or more points are equal distance apart. Please resolve'
51                   raise msg
52         indices.remove(index)
53         boundary[i+1][0]=x[index]
54         boundary[i+1][1]=y[index]
55         order.append(index)
56    #ensure boundary is a polygon i.e no intesections
57    #e.g.
58    #msg='Boundary is not a regular polygon. Intersections found'
59    #assert boundary.is_polygon, msg
60             
61    #So far i cannot tell if the wrong point is closer.
62   
63    boundary=ensure_numeric(boundary)
64    from anuga.coordinate_transforms import redfearn,convert_from_latlon_to_utm
65    boundary_utm,zone=convert_from_latlon_to_utm(longitudes=list(boundary[:,0]),latitudes=list(boundary[:,1]))
66
67    d=","
68    fid=open(order_filename,'w')
69    header='index, longitude, latitude\n'
70    fid.write(header)
71    line=str(west_index)+d+str(x[west_index])+d+\
72        str(y[west_index])+"\n"
73    fid.write(line)
74    for i in order:
75            line=str(i)+d+str(x[i])+d+\
76                str(y[i])+"\n"
77            fid.write(line)
78    fid.close()
79    return boundary,boundary_utm
80
81######################
82# Domain definitions #
83######################
84
85# Bathymetry and topography filenames
86dir='data'
87
88time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
89poor_combined_dir_name = dir+sep+'poor_phuket_bathymetry'
90good_combined_dir_name = dir+sep+'good_phuket_bathymetry'
91meshname = 'phuket_mesh'
92mesh_elevname = 'phuket_mesh_elev'
93boundary_most_in = dir+sep+'out'
94#boundary_most_out = dir+sep+'most_boundary_condition'
95boundary_most_out = dir+sep+'tide_most_boundary_condition'
96
97dir='mesh_polygons'
98extent = 'extent.csv'
99patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
100   
101tide=0.35
102
103######################################
104# Create urs boundary from mux2files #
105######################################
106urs_filename='polyline'
107base_name='tide_polyline'
108order_filename='boxingday_boundary_order.txt'
109
110# Create ordering file
111# Usually this would be provided before hand
112boundary,boundary_utm=create_sts_boundary_order_file(urs_filename+'-z-mux2',order_filename,base_name,verbose=False)
113
114# Create ordered sts file
115#if not os.path.exists(base_name+'.sts'):
116print 'creating sts file'
117urs2sts(urs_filename,
118        basename_out=base_name,
119        ordering_filename=order_filename,
120        mean_stage=tide,
121        verbose=True)
122#else:
123#    print 'sts file exists'
124print 'done'
125
126# Read in boundary from ordered sts file
127urs_boundary=create_sts_boundary(base_name)
128
129# Checks that the boundary read from the sts file is the same as the one
130# defined by the function create_sts_boundary_order_file
131assert allclose(boundary_utm,urs_boundary)
132
133#############################################################
134# Add urs boundary segment to user defined boundary polygon #
135#############################################################
136
137# Read in specific user defined boundary
138dir='mesh_polygons'
139extent = 'extent.csv'
140user_boundary_polygon = read_polygon(dir+sep+extent)
141
142# The following is specific to this scenario
143indices=inside_polygon(urs_boundary, user_boundary_polygon, closed=True, verbose=False)
144bounding_polygon=[]
145for i in indices[:-2]:
146    #disregard last two entries because they are not on boundary
147    bounding_polygon.append(urs_boundary[i])
148
149# Append clipped sts boundary to user defined domain boundary
150bounding_polygon.extend(user_boundary_polygon[2:len(user_boundary_polygon)])
151
152###############################
153# Interior region definitions #
154###############################
155
156contour20m = read_polygon(dir+sep+'20m.csv')
157contour20m[0][1]+=3000
158contour20m[-2][0]+=3000
159contour20m[-2][1]-=13000
160contour20m[-1][1]+=9000
161patong = read_polygon(dir+sep+'patong.csv')
162bay = read_polygon(dir+sep+'bay.csv')
163
164##########################################
165# Plot urs boundary and internal regions #
166##########################################
167
168plot=False
169if plot:
170    from pylab import plot,show,hold
171    bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only
172    #boundary=ensure_numeric(boundary)
173    #plot(boundary[:,0],boundary[:,1],'o')
174    plot(bounding_polygon[:,0],bounding_polygon[:,1],'o')
175    plot(bounding_polygon[:,0],bounding_polygon[:,1])
176    contour20m = ensure_numeric(contour20m)
177    plot(contour20m[:,0],contour20m[:,1])
178    patong=ensure_numeric(patong)
179    plot(patong[:,0],patong[:,1])
180    show()
181
182   
183#######################
184# For MOST SIMULATION #
185#######################
186east  = 97.7
187west  = 96.7
188north = 9.0
189south = 5.9
190
191
192bay_west = 417348.0
193bay_east = 425656.0
194bay_south = 870565.0
195bay_north = 877008.0
196
197bay_grid_poly=[[bay_west,bay_south],[bay_west,bay_north],[bay_east,bay_north],[bay_east,bay_south]]
198
199check_ferret_extent=False
200if check_ferret_extent:
201    east  = east
202    west  = west
203    north = north
204    south = south
205   
206    from anuga.coordinate_transforms.redfearn import redfearn
207    zone0, easting0, northing0 = redfearn(north,east)
208    print zone0,easting0,northing0
209    zone1, easting1, northing1 = redfearn(south,east)
210    print zone1,easting1,northing1
211    zone2, easting2, northing2 = redfearn(south,west)
212    print zone2,easting2,northing2
213    zone3, easting3, northing3 = redfearn(north,west)
214    print zone3,easting3,northing3
215   
216    assert zone0==zone1==zone2==zone3
217   
218    ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]]
219    points_outside=[]
220    from anuga.utilities.polygon import is_inside_polygon
221    for i,point in enumerate(bounding_polygon):
222        if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False):
223            print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1])
224            points_outside.append(point)
225
226    from pylab import plot,show
227    ferret_polygon=ensure_numeric(ferret_polygon)
228    points_outside=ensure_numeric(points_outside)
229    plot(ferret_polygon[:,0],ferret_polygon[:,1])
230    plot(points_outside[:,0],points_outside[:,1],'o')
231    show()
232    print ferret_polygon[0]
233    print ferret_polygon[1]
234    print ferret_polygon[2]
235    print ferret_polygon[3]
236    print points_outside
237   
Note: See TracBrowser for help on using the repository browser.