1 | def 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 | |
---|
58 | from anuga.shallow_water.data_manager import urs2sts |
---|
59 | base_name='polyline' |
---|
60 | import os |
---|
61 | if os.path.exists(base_name+'.sts'): |
---|
62 | print 'sts boundary file already created.' |
---|
63 | else: |
---|
64 | urs2sts(base_name,mean_stage=0.0,verbose=False) |
---|
65 | boundary=create_sts_boundary(base_name,verbose=True) |
---|
66 | bounding_polygon=boundary.tolist() |
---|
67 | del bounding_polygon[0:2] #delete first two elements |
---|
68 | n=len(bounding_polygon) |
---|
69 | del bounding_polygon[n-8:n] |
---|
70 | |
---|
71 | import project |
---|
72 | bp=project.bounding_polygon |
---|
73 | bounding_polygon.insert(0,bp[0]) |
---|
74 | bounding_polygon.extend(bp[len(bp)-5:len(bp)]) |
---|
75 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
76 | bounding_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 | |
---|
83 | check_ferret_extent=False |
---|
84 | if 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 | |
---|