# source:anuga_work/development/boxingday08/project.py@5458

Last change on this file since 5458 was 5458, checked in by jakeman, 15 years ago

updates of scripts used for boxingday validation to include boundary ordering

File size: 6.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
6from os import sep
7from time import localtime, strftime, gmtime
8
9def create_sts_boundary_order_file(order_filename,stsfilename,verbose=False):
10    """
11    Note This is not a robust algorithm. Be Careful when applying that
12    polygon does not intersect itself and that resulting boundary actually
13    follows the path it is supposed too.
14    """
15    fid = NetCDFFile(stsfilename+'.sts', 'r')    #Open existing file for read
16    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
17    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
18    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
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    d=","
63    fid=open(order_filename,'w')
66    line=str(west_index)+d+str(x[west_index])+d+\
67        str(y[west_index])+"\n"
68    fid.write(line)
69    for i in order:
70            line=str(i)+d+str(x[i])+d+\
71                str(y[i])+"\n"
72            fid.write(line)
73    fid.close()
74    return boundary
75
76
78from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
79urs_name='polyline'
80#base_name=urs_name
81tide=0.35
82base_name='tide_polyline'
83#base_name='most_polyline'
84import os
85if os.path.exists(base_name+'.sts'):
86    print 'sts boundary file already created.'
87else:
88    print 'creatin sts file'
89    urs2sts(urs_name,base_name,mean_stage=tide,verbose=False)
90
91order_filename='sts_order.txt'
92boundary=create_sts_boundary_order_file(order_filename,base_name,verbose=False)
93urs_boundary=create_sts_boundary(order_filename,base_name,lat_long=False)
94assert allclose(boundary.tolist(),urs_boundary)
95
97dir='mesh_polygons'
98extent = 'extent.csv'
100
101#Now order of boundary is correct clip sts boundary to small region for
102#this particular scenario
103indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
104bounding_polygon=[]
105for i in indices[:-2]:
106    #disregard last two entries because they are not on boundary
107    bounding_polygon.append(urs_boundary[i])
108
109#Append clipped sts boundary to large domain boundary
110bounding_polygon.extend(bp[2:len(bp)])
111
112###############################
113# Domain definitions
114###############################
115
116# Bathymetry and topography filenames
117dir='data'
118
119time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
120poor_combined_dir_name = dir+sep+'poor_phuket_bathymetry'
121good_combined_dir_name = dir+sep+'good_phuket_bathymetry'
122meshname = 'phuket_mesh'
123mesh_elevname = 'phuket_mesh_elev'
124boundary_most_in = dir+sep+'out'
125#boundary_most_out = dir+sep+'most_boundary_condition'
126boundary_most_out = dir+sep+'tide_most_boundary_condition'
127
128dir='mesh_polygons'
129extent = 'extent.csv'
131
132###############################
133# Interior region definitions
134###############################
135
137contour20m[0][1]+=3000
138contour20m[-2][0]+=3000
139contour20m[-2][1]-=13000
140contour20m[-1][1]+=9000
143
144plot=False
145if plot:
146    from pylab import plot,show,hold
147    bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only
148    #boundary=ensure_numeric(boundary)
149    #plot(boundary[:,0],boundary[:,1],'o')
150    plot(bounding_polygon[:,0],bounding_polygon[:,1],'o')
151    plot(bounding_polygon[:,0],bounding_polygon[:,1])
152    contour20m = ensure_numeric(contour20m)
153    plot(contour20m[:,0],contour20m[:,1])
154    patong=ensure_numeric(patong)
155    plot(patong[:,0],patong[:,1])
156    show()
157
158east  = 97.7
159west  = 96.7
160north = 9.0
161south = 5.9
162
163
164bay_west = 417348.0
165bay_east = 425656.0
166bay_south = 870565.0
167bay_north = 877008.0
168
169bay_grid_poly=[[bay_west,bay_south],[bay_west,bay_north],[bay_east,bay_north],[bay_east,bay_south]]
170
171check_ferret_extent=False
172if check_ferret_extent:
173    east  = east
174    west  = west
175    north = north
176    south = south
177
178    from anuga.coordinate_transforms.redfearn import redfearn
179    zone0, easting0, northing0 = redfearn(north,east)
180    print zone0,easting0,northing0
181    zone1, easting1, northing1 = redfearn(south,east)
182    print zone1,easting1,northing1
183    zone2, easting2, northing2 = redfearn(south,west)
184    print zone2,easting2,northing2
185    zone3, easting3, northing3 = redfearn(north,west)
186    print zone3,easting3,northing3
187
188    assert zone0==zone1==zone2==zone3
189
190    ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]]
191    points_outside=[]
192    from anuga.utilities.polygon import is_inside_polygon
193    for i,point in enumerate(bounding_polygon):
194        if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False):
195            print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1])
196            points_outside.append(point)
197
198    from pylab import plot,show
199    ferret_polygon=ensure_numeric(ferret_polygon)
200    points_outside=ensure_numeric(points_outside)
201    plot(ferret_polygon[:,0],ferret_polygon[:,1])
202    plot(points_outside[:,0],points_outside[:,1],'o')
203    show()
204    print ferret_polygon[0]
205    print ferret_polygon[1]
206    print ferret_polygon[2]
207    print ferret_polygon[3]
208    print points_outside
209
Note: See TracBrowser for help on using the repository browser.