source: anuga_work/development/boxingday08/create_boundary.py @ 5393

Last change on this file since 5393 was 5393, checked in by jakeman, 16 years ago

updating scripts for boxingday validation

File size: 4.3 KB
Line 
1def create_sts_boundary(filename,verbose=False):
2    from Scientific.IO.NetCDF import NetCDFFile
3    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
4        compress,zeros,fabs
5    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
6    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
7    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
8    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
9   
10    #find indices of the southernmost gauges
11    #if more than one find gauges furthest left
12    #will fail if domain crosses 180 degrees
13    #i.e. if left most most point is not west i.e x=1 and other points 170<x<180
14    #or even in utm just a diffenet zone
15    south=min(y)
16    south_indices=compress(y-south<1e-6,arange(len(y)))
17    west=x[south_indices[0]]
18    west_index=south_indices[0]
19    if south_indices>1:
20         #find most western gauge
21         for i,south_index in enumerate(south_indices[1:]):
22              if x[south_index]<west:
23                   west_index=south_index
24   
25    n=len(x)
26    boundary=zeros((n,2),Float)
27    indices=arange(n).tolist()
28    boundary[0][0]=x[west_index]
29    boundary[0][1]=y[west_index]
30    indices.remove(west_index)
31    for i in range(n-1):
32         min_dist=1e12
33         for j,k in enumerate(indices):
34              dist=sqrt((boundary[i][0]-x[k])**2+(boundary[i][1]-y[k])**2)
35              if dist<min_dist:
36                   min_dist=dist
37                   index=k
38              elif fabs(dist-min_dist)<1e-6:
39                   msg='Two or more points are equal distance apart. Please resolve'
40                   raise msg
41         indices.remove(index)
42         boundary[i+1][0]=x[index]
43         boundary[i+1][1]=y[index]
44    #ensure boundary is a polygon i.e no intesections
45    #e.g.
46    #msg='Boundary is not a regular polygon. Intersections found'
47    #assert boundary.is_polygon, msg
48             
49    #So far i cannot tell if the wrong point is closer.
50
51    #if verbose:
52         #from pylab import plot,show
53         #plot(x,y,'o')
54         #plot(boundary[:,0],boundary[:,1])
55         #show()
56    return boundary
57
58from anuga.shallow_water.data_manager import urs2sts
59base_name='polyline'
60import os
61if os.path.exists(base_name+'.sts'):
62        print 'sts boundary file already created.'
63else:
64    urs2sts(base_name,mean_stage=0.0,verbose=False)
65boundary=create_sts_boundary(base_name,verbose=True)
66bounding_polygon=boundary.tolist()
67del bounding_polygon[0:2] #delete first two elements
68n=len(bounding_polygon)
69del bounding_polygon[n-8:n]
70
71import project
72bp=project.bounding_polygon
73bounding_polygon.insert(0,bp[0])
74bounding_polygon.extend(bp[len(bp)-5:len(bp)])
75from anuga.utilities.numerical_tools import ensure_numeric
76bounding_polygon=ensure_numeric(bounding_polygon)
77#from pylab import plot,show
78#plot(bounding_polygon[:,0],bounding_polygon[:,1])
79#bp=ensure_numeric(bp)
80#plot(bp[:,0],bp[:,1])
81#show()
82
83check_ferret_extent=False
84if check_ferret_extent:
85    east  = project.east
86    west  = project.west
87    north = project.north
88    south = project.south
89   
90    from anuga.coordinate_transforms.redfearn import redfearn
91    zone0, easting0, northing0 = redfearn(north,east)
92    print zone0,easting0,northing0
93    zone1, easting1, northing1 = redfearn(south,east)
94    print zone1,easting1,northing1
95    zone2, easting2, northing2 = redfearn(south,west)
96    print zone2,easting2,northing2
97    zone3, easting3, northing3 = redfearn(north,west)
98    print zone3,easting3,northing3
99   
100    assert zone0==zone1==zone2==zone3
101   
102    ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]]
103    points_outside=[]
104    from anuga.utilities.polygon import is_inside_polygon
105    for i,point in enumerate(bounding_polygon):
106        if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False):
107            print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1])
108            points_outside.append(point)
109
110    from pylab import plot,show
111    ferret_polygon=ensure_numeric(ferret_polygon)
112    points_outside=ensure_numeric(points_outside)
113    plot(ferret_polygon[:,0],ferret_polygon[:,1])
114    plot(points_outside[:,0],points_outside[:,1],'o')
115    show()
116    print ferret_polygon[0]
117    print ferret_polygon[1]
118    print ferret_polygon[2]
119    print ferret_polygon[3]
120    print points_outside
121   
Note: See TracBrowser for help on using the repository browser.