Changeset 5402
- Timestamp:
- Jun 13, 2008, 10:37:14 AM (15 years ago)
- Location:
- anuga_work/development/boxingday08
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/boxingday08/get_gauges.py
r5393 r5402 1 1 import os 2 2 from os import sep 3 import create_boundary 4 from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold 3 #import create_boundary 4 import project 5 try: 6 from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold 7 except: 8 print 'Cannot import pylab plotting will not work. Csv files are still created' 5 9 from anuga.utilities.polygon import inside_polygon 6 10 from anuga.utilities.numerical_tools import ensure_numeric … … 16 20 elevation=fid.variables['elevation'][:] 17 21 indices=inside_polygon(points, polygon, closed=True, verbose=False) 18 inside_points=[] 19 for i in indices: 20 inside_points.append(points[i]) 21 22 inside_points=take(points,indices,0) 23 elevation=take(elevation,indices,0) 22 24 if verbose: 23 25 inside_points_array=ensure_numeric(inside_points)#for plotting only 24 #plot(points[:,0],points[:,1],'o')26 plot(points[:,0],points[:,1],'o') 25 27 plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o') 28 polygon=ensure_numeric(polygon) 26 29 plot(polygon[:,0],polygon[:,1]) 27 30 show() 28 31 29 32 quantity_names=['stage','xmomentum','ymomentum'] 30 33 quantities = {} … … 38 41 fid.close() 39 42 40 return points,quantities,elevation,time43 return inside_points,quantities,elevation,time 41 44 42 45 from anuga.shallow_water.data_manager import urs2sts … … 48 51 urs2sts(base_name,mean_stage=0.0,verbose=False) 49 52 50 polygon=create_boundary.bounding_polygon 51 points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=True) 52 53 #polygon=create_boundary.bounding_polygon 54 polygon=project.bounding_polygon 55 points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=False) 56 assert len(points)==len(elevation)==len(quantities['stage'][0,:]) 53 57 54 58 def plot_gauge(time,quantity,quantity_name,elevation,gauge_id,pos): 55 figname=figname = '%s_%s%s.png' %(' urs_gauge',59 figname=figname = '%s_%s%s.png' %('tgs/urs_gauge', 56 60 quantity_name, 57 61 str(gauge_id)) … … 60 64 title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1])) 61 65 xlabel('Time (minutes)') 62 ylabel( 'Stage (m)')66 ylabel(quantity_name) 63 67 legend() 64 68 print 'saving figure here %s' %figname … … 67 71 points=points[::5] 68 72 elevation=elevation[::5] 69 print points.shape70 73 gauge_filename='urs_gauges' 71 74 d=',' … … 91 94 92 95 print 'generating graphs' 93 csv2timeseries_graphs(directories_dic={dir_name:[ 'dir_name',0, 0]},96 csv2timeseries_graphs(directories_dic={dir_name:[dir_name,0, 0]}, 94 97 base_name='gauge_', 95 98 plot_numbers='',#['1'], 96 quantities=['stage' ],99 quantities=['stage','xmomentum','ymomentum'], 97 100 extra_plot_name='', 98 101 assess_all_csv_files=True, … … 100 103 verbose=True) 101 104 102 #sww_filename='most/good_simulationsourcefriction20080529_160144' 103 #sww_filename='cloud/good_simulation_ursfriction20080529_160110' 104 sww_filename='polyline/good_simulation_polylinefriction20080529_115703' 105 #sww_filename='old_cloud/good_simulation_ursfriction20080509_143712' 106 get_sww_gauge_data(sww_filename,gauge_filename) 107 #for i in range(len(points)): 108 # quantity=quantities['stage'][:,i] 109 # plot_gauge(time,quantity,'stage',elevation[i],i,points[i]) 105 sww_filename='polyline/good_polyline20080607_205151' 106 #sww_filename='most/good_most20080607_205101' 107 #sww_filename='cloud/good_urs20080607_205134' 108 #get_sww_gauge_data(sww_filename,gauge_filename) 109 for i in range(len(points)): 110 quantity=quantities['stage'][:,i] 111 plot_gauge(time,quantity,'stage',elevation[i],i,points[i]) 112 quantity=quantities['xmomentum'][:,i] 113 plot_gauge(time,quantity,'xmomentum',elevation[i],i,points[i]) 114 #quantity=quantities['ymomentum'][:,i] 115 #plot_gauge(time,quantity,'ymomentum',elevation[i],i,points[i]) -
anuga_work/development/boxingday08/project.py
r5393 r5402 1 """Common filenames and locations for topographic data, meshes and outputs. 2 Author: John Jakeman, The Australian National University (2008) 3 """ 1 from anuga.utilities.numerical_tools import ensure_numeric 2 from Scientific.IO.NetCDF import NetCDFFile 3 from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ 4 compress,zeros,fabs 5 from anuga.utilities.polygon import inside_polygon,read_polygon 6 from os import sep 7 from time import localtime, strftime, gmtime 8 def create_sts_boundary(filename,verbose=False): 9 fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read 10 x = fid.variables['x'][:]+fid.xllcorner #x-coordinates of vertices 11 y = fid.variables['y'][:]+fid.yllcorner #y-coordinates of vertices 12 coordinates=transpose(asarray([x.tolist(),y.tolist()])) 13 14 #find indices of the southernmost gauges 15 #if more than one find gauges furthest left 16 #will fail if domain crosses 180 degrees 17 #i.e. if left most most point is not west i.e x=1 and other points 170<x<180 18 #or even in utm just a diffenet zone 19 south=min(y) 20 south_indices=compress(y-south<1e-6,arange(len(y))) 21 west=x[south_indices[0]] 22 west_index=south_indices[0] 23 if south_indices>1: 24 #find most western gauge 25 for i,south_index in enumerate(south_indices[1:]): 26 if x[south_index]<west: 27 west_index=south_index 28 29 n=len(x) 30 boundary=zeros((n,2),Float) 31 indices=arange(n).tolist() 32 boundary[0][0]=x[west_index] 33 boundary[0][1]=y[west_index] 34 indices.remove(west_index) 35 for i in range(n-1): 36 min_dist=1e12 37 for j,k in enumerate(indices): 38 dist=sqrt((boundary[i][0]-x[k])**2+(boundary[i][1]-y[k])**2) 39 if dist<min_dist: 40 min_dist=dist 41 index=k 42 elif fabs(dist-min_dist)<1e-6: 43 msg='Two or more points are equal distance apart. Please resolve' 44 raise msg 45 indices.remove(index) 46 boundary[i+1][0]=x[index] 47 boundary[i+1][1]=y[index] 48 #ensure boundary is a polygon i.e no intesections 49 #e.g. 50 #msg='Boundary is not a regular polygon. Intersections found' 51 #assert boundary.is_polygon, msg 52 53 #So far i cannot tell if the wrong point is closer. 54 return boundary 4 55 5 from os import sep, environ, getenv, getcwd 6 from os.path import expanduser 7 import sys 8 from time import localtime, strftime, gmtime 9 from anuga.utilities.polygon import read_polygon, plot_polygons, \ 10 polygon_area, is_inside_polygon 11 from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm 12 from time import localtime, strftime 56 #read in polyline boundary 57 from anuga.shallow_water.data_manager import urs2sts 58 base_name='polyline' 59 import os 60 if os.path.exists(base_name+'.sts'): 61 print 'sts boundary file already created.' 62 else: 63 urs2sts(base_name,mean_stage=0.0,verbose=False) 64 boundary=create_sts_boundary(base_name,verbose=True) 65 urs_boundary=boundary.tolist() 13 66 14 def convert_arcgis_latlon_list_to_utm(points): 15 #Used because arc gis produced csv files put lat lon in 16 #reverse order to those accpeted by convert_latlon_to_utm() 17 18 from anuga.coordinate_transforms.geo_reference import Geo_reference 19 from anuga.coordinate_transforms.redfearn import redfearn 20 old_geo = Geo_reference() 21 utm_points = [] 22 for point in points: 23 zone, easting, northing = redfearn(float(point[1]),float(point[0])) 24 new_geo = Geo_reference(zone) 25 old_geo.reconcile_zones(new_geo) 26 utm_points.append([easting,northing]) 67 #Read in large domain boundary 68 dir='mesh_polygons' 69 extent = 'extent.csv' 70 bp = read_polygon(dir+sep+extent) 27 71 28 return utm_points, old_geo.get_zone() 72 indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False) 73 bounding_polygon=[] 74 for i in indices[:-2]: 75 #disregard last two entries because they are not on boundary 76 bounding_polygon.append(urs_boundary[i]) 77 78 bounding_polygon.extend(bp[2:len(bp)]) 29 79 30 80 ############################### … … 45 95 dir='mesh_polygons' 46 96 extent = 'extent.csv' 47 #extent = 'boundary_new_pts.csv'#old boundary48 # bounding polygon for study area49 bounding_polygon = read_polygon(dir+sep+extent)50 97 patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv') 51 98 … … 54 101 ############################### 55 102 56 island_south = read_polygon(dir+sep+'south_island.csv') 57 island_south2 = read_polygon(dir+sep+'south_island2.csv') 58 island_north = read_polygon(dir+sep+'north_island.csv') 103 contour20m = read_polygon(dir+sep+'20m.csv') 104 contour20m[0][1]+=3000 105 contour20m[-2][0]+=3000 106 contour20m[-2][1]-=13000 107 contour20m[-1][1]+=9000 108 patong = read_polygon(dir+sep+'patong.csv') 59 109 bay = read_polygon(dir+sep+'bay.csv') 60 patong = read_polygon(dir+sep+'patong.csv')61 contour20m = read_polygon(dir+sep+'20m.csv')62 110 63 64 65 #max_extent_of MOST file 66 #east = 94.93325 67 #west = 80.0003 68 #north = 10.0000 69 #south = -4.9330 70 71 #max_extent_of MOST file should be 72 #east = 104.0 73 #west = 80.0 74 #north = 10.0 75 #south = -5.0 76 77 #max_extent of old simulation is why does this work and not below 78 #7.40911081272 8.71484358635 99.1683687224 97.513856322 79 80 #max extent of clipped MOST file 111 plot=False 112 if plot: 113 from pylab import plot,show,hold 114 bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only 115 plot(bounding_polygon[:,0],bounding_polygon[:,1],'o') 116 plot(bounding_polygon[:,0],bounding_polygon[:,1]) 117 contour20m = ensure_numeric(contour20m) 118 plot(contour20m[:,0],contour20m[:,1]) 119 patong=ensure_numeric(patong) 120 plot(patong[:,0],patong[:,1]) 121 show() 81 122 82 123 east = 97.7 … … 84 125 north = 9.0 85 126 south = 5.9 86 87 127 128 check_ferret_extent=False 129 if check_ferret_extent: 130 east = east 131 west = west 132 north = north 133 south = south 134 135 from anuga.coordinate_transforms.redfearn import redfearn 136 zone0, easting0, northing0 = redfearn(north,east) 137 print zone0,easting0,northing0 138 zone1, easting1, northing1 = redfearn(south,east) 139 print zone1,easting1,northing1 140 zone2, easting2, northing2 = redfearn(south,west) 141 print zone2,easting2,northing2 142 zone3, easting3, northing3 = redfearn(north,west) 143 print zone3,easting3,northing3 144 145 assert zone0==zone1==zone2==zone3 146 147 ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]] 148 points_outside=[] 149 from anuga.utilities.polygon import is_inside_polygon 150 for i,point in enumerate(bounding_polygon): 151 if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False): 152 print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1]) 153 points_outside.append(point) 154 155 from pylab import plot,show 156 ferret_polygon=ensure_numeric(ferret_polygon) 157 points_outside=ensure_numeric(points_outside) 158 plot(ferret_polygon[:,0],ferret_polygon[:,1]) 159 plot(points_outside[:,0],points_outside[:,1],'o') 160 show() 161 print ferret_polygon[0] 162 print ferret_polygon[1] 163 print ferret_polygon[2] 164 print ferret_polygon[3] 165 print points_outside 166
Note: See TracChangeset
for help on using the changeset viewer.