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

Last change on this file since 5543 was 5543, checked in by ole, 14 years ago

Rerunning boxing day validation work

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