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

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

Name changes and change to good bathymetry

File size: 7.8 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,basename_out=base_name,
118        ordering_filename=order_filename,
119        mean_stage=tide,
120        verbose=False)
121#else:
122#    print 'sts file exists'
123
124# Read in boundary from ordered sts file
125urs_boundary=create_sts_boundary(base_name)
126
127# Checks that the boundary read from the sts file is the same as the one
128# defined by the function create_sts_boundary_order_file
129assert allclose(boundary_utm,urs_boundary)
130
131#############################################################
132# Add urs boundary segment to user defined boundary polygon #
133#############################################################
134
135# Read in specific user defined boundary
136dir='mesh_polygons'
137extent = 'extent.csv'
138user_boundary_polygon = read_polygon(dir+sep+extent)
139
140# The following is specific to this scenario
141indices=inside_polygon(urs_boundary, user_boundary_polygon, closed=True, verbose=False)
142bounding_polygon=[]
143for i in indices[:-2]:
144    #disregard last two entries because they are not on boundary
145    bounding_polygon.append(urs_boundary[i])
146
147# Append clipped sts boundary to user defined domain boundary
148bounding_polygon.extend(user_boundary_polygon[2:len(user_boundary_polygon)])
149
150###############################
151# Interior region definitions #
152###############################
153
154contour20m = read_polygon(dir+sep+'20m.csv')
155contour20m[0][1]+=3000
156contour20m[-2][0]+=3000
157contour20m[-2][1]-=13000
158contour20m[-1][1]+=9000
159patong = read_polygon(dir+sep+'patong.csv')
160bay = read_polygon(dir+sep+'bay.csv')
161
162##########################################
163# Plot urs boundary and internal regions #
164##########################################
165
166plot=False
167if plot:
168    from pylab import plot,show,hold
169    bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only
170    #boundary=ensure_numeric(boundary)
171    #plot(boundary[:,0],boundary[:,1],'o')
172    plot(bounding_polygon[:,0],bounding_polygon[:,1],'o')
173    plot(bounding_polygon[:,0],bounding_polygon[:,1])
174    contour20m = ensure_numeric(contour20m)
175    plot(contour20m[:,0],contour20m[:,1])
176    patong=ensure_numeric(patong)
177    plot(patong[:,0],patong[:,1])
178    show()
179
180   
181#######################
182# For MOST SIMULATION #
183#######################
184east  = 97.7
185west  = 96.7
186north = 9.0
187south = 5.9
188
189
190bay_west = 417348.0
191bay_east = 425656.0
192bay_south = 870565.0
193bay_north = 877008.0
194
195bay_grid_poly=[[bay_west,bay_south],[bay_west,bay_north],[bay_east,bay_north],[bay_east,bay_south]]
196
197check_ferret_extent=False
198if check_ferret_extent:
199    east  = east
200    west  = west
201    north = north
202    south = south
203   
204    from anuga.coordinate_transforms.redfearn import redfearn
205    zone0, easting0, northing0 = redfearn(north,east)
206    print zone0,easting0,northing0
207    zone1, easting1, northing1 = redfearn(south,east)
208    print zone1,easting1,northing1
209    zone2, easting2, northing2 = redfearn(south,west)
210    print zone2,easting2,northing2
211    zone3, easting3, northing3 = redfearn(north,west)
212    print zone3,easting3,northing3
213   
214    assert zone0==zone1==zone2==zone3
215   
216    ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]]
217    points_outside=[]
218    from anuga.utilities.polygon import is_inside_polygon
219    for i,point in enumerate(bounding_polygon):
220        if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False):
221            print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1])
222            points_outside.append(point)
223
224    from pylab import plot,show
225    ferret_polygon=ensure_numeric(ferret_polygon)
226    points_outside=ensure_numeric(points_outside)
227    plot(ferret_polygon[:,0],ferret_polygon[:,1])
228    plot(points_outside[:,0],points_outside[:,1],'o')
229    show()
230    print ferret_polygon[0]
231    print ferret_polygon[1]
232    print ferret_polygon[2]
233    print ferret_polygon[3]
234    print points_outside
235   
Note: See TracBrowser for help on using the repository browser.