Changeset 5402

Timestamp:

Jun 13, 2008, 10:37:14 AM (15 years ago)

Author:

ole

Message:

Project file and get_gauges to support boxing day validation scripts committed in changeset:5401

Location:

anuga_work/development/boxingday08

Files:

• get_gauges.py (modified) (8 diffs)
• project.py (modified) (4 diffs)

anuga_work/development/boxingday08/get_gauges.py

@@ -1 +1 @@
 import os
 from os import sep
-import create_boundary
-from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
+#import create_boundary
+import project
+try:
+    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
+except:
+    print 'Cannot import pylab plotting will not work. Csv files are still created'
 from anuga.utilities.polygon import inside_polygon
 from anuga.utilities.numerical_tools import ensure_numeric
@@ -16 +20 @@
     elevation=fid.variables['elevation'][:]
     indices=inside_polygon(points, polygon, closed=True, verbose=False)
-    inside_points=[]
-    for i in indices:
-        inside_points.append(points[i])
-
+    inside_points=take(points,indices,0)
+    elevation=take(elevation,indices,0)
     if verbose:
         inside_points_array=ensure_numeric(inside_points)#for plotting only
-        #plot(points[:,0],points[:,1],'o')
+        plot(points[:,0],points[:,1],'o')
         plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o')
+        polygon=ensure_numeric(polygon)
         plot(polygon[:,0],polygon[:,1])
         show()
-
+        
     quantity_names=['stage','xmomentum','ymomentum']
     quantities = {}
@@ -38 +41 @@
     fid.close()
 
-    return points,quantities,elevation,time
+    return inside_points,quantities,elevation,time
 
 from anuga.shallow_water.data_manager import urs2sts
@@ -48 +51 @@
     urs2sts(base_name,mean_stage=0.0,verbose=False)
 
-polygon=create_boundary.bounding_polygon
-points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=True)
-
+#polygon=create_boundary.bounding_polygon
+polygon=project.bounding_polygon
+points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=False)
+assert len(points)==len(elevation)==len(quantities['stage'][0,:])
 
 def plot_gauge(time,quantity,quantity_name,elevation,gauge_id,pos):
-    figname=figname = '%s_%s%s.png' %('urs_gauge',
+    figname=figname = '%s_%s%s.png' %('tgs/urs_gauge',
                                       quantity_name,
                                       str(gauge_id))
@@ -60 +64 @@
     title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1]))
     xlabel('Time (minutes)')
-    ylabel('Stage (m)')
+    ylabel(quantity_name)
     legend()
     print 'saving figure here %s' %figname
@@ -67 +71 @@
 points=points[::5]
 elevation=elevation[::5]
-print points.shape
 gauge_filename='urs_gauges'
 d=','
@@ -91 +94 @@
         
     print 'generating graphs'
-    csv2timeseries_graphs(directories_dic={dir_name:['dir_name',0, 0]},
+    csv2timeseries_graphs(directories_dic={dir_name:[dir_name,0, 0]},
                           base_name='gauge_',
                           plot_numbers='',#['1'],
-                          quantities=['stage'],
+                          quantities=['stage','xmomentum','ymomentum'],
                           extra_plot_name='',
                           assess_all_csv_files=True,
@@ -100 +103 @@
                           verbose=True)
 
-#sww_filename='most/good_simulationsourcefriction20080529_160144'
-#sww_filename='cloud/good_simulation_ursfriction20080529_160110'
-sww_filename='polyline/good_simulation_polylinefriction20080529_115703'
-#sww_filename='old_cloud/good_simulation_ursfriction20080509_143712'
-get_sww_gauge_data(sww_filename,gauge_filename)
-#for i in range(len(points)):
-#    quantity=quantities['stage'][:,i]
-#    plot_gauge(time,quantity,'stage',elevation[i],i,points[i])
+sww_filename='polyline/good_polyline20080607_205151'
+#sww_filename='most/good_most20080607_205101'
+#sww_filename='cloud/good_urs20080607_205134'
+#get_sww_gauge_data(sww_filename,gauge_filename)
+for i in range(len(points)):
+    quantity=quantities['stage'][:,i]
+    plot_gauge(time,quantity,'stage',elevation[i],i,points[i])
+    quantity=quantities['xmomentum'][:,i]
+    plot_gauge(time,quantity,'xmomentum',elevation[i],i,points[i])
+    #quantity=quantities['ymomentum'][:,i]
+    #plot_gauge(time,quantity,'ymomentum',elevation[i],i,points[i])

anuga_work/development/boxingday08/project.py

@@ -1 @@
-"""Common filenames and locations for topographic data, meshes and outputs.
-Author: John Jakeman, The Australian National University (2008)
-"""
+from anuga.utilities.numerical_tools import ensure_numeric
+from Scientific.IO.NetCDF import NetCDFFile
+from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
+     compress,zeros,fabs
+from anuga.utilities.polygon import inside_polygon,read_polygon
+from os import sep
+from time import localtime, strftime, gmtime
+def create_sts_boundary(filename,verbose=False):
+    fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
+    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
+    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
+    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
+    
+    #find indices of the southernmost gauges
+    #if more than one find gauges furthest left
+    #will fail if domain crosses 180 degrees
+    #i.e. if left most most point is not west i.e x=1 and other points 170<x<180
+    #or even in utm just a diffenet zone
+    south=min(y)
+    south_indices=compress(y-south<1e-6,arange(len(y)))
+    west=x[south_indices[0]]
+    west_index=south_indices[0]
+    if south_indices>1:
+         #find most western gauge
+         for i,south_index in enumerate(south_indices[1:]):
+              if x[south_index]<west:
+                   west_index=south_index
+    
+    n=len(x)
+    boundary=zeros((n,2),Float)
+    indices=arange(n).tolist()
+    boundary[0][0]=x[west_index]
+    boundary[0][1]=y[west_index]
+    indices.remove(west_index)
+    for i in range(n-1):
+         min_dist=1e12
+         for j,k in enumerate(indices):
+              dist=sqrt((boundary[i][0]-x[k])**2+(boundary[i][1]-y[k])**2)
+              if dist<min_dist:
+                   min_dist=dist
+                   index=k
+              elif fabs(dist-min_dist)<1e-6:
+                   msg='Two or more points are equal distance apart. Please resolve'
+                   raise msg
+         indices.remove(index)
+         boundary[i+1][0]=x[index]
+         boundary[i+1][1]=y[index]
+    #ensure boundary is a polygon i.e no intesections
+    #e.g.
+    #msg='Boundary is not a regular polygon. Intersections found'
+    #assert boundary.is_polygon, msg
+              
+    #So far i cannot tell if the wrong point is closer.
+    return boundary
 
-from os import sep, environ, getenv, getcwd
-from os.path import expanduser
-import sys
-from time import localtime, strftime, gmtime
-from anuga.utilities.polygon import read_polygon, plot_polygons, \
-                                    polygon_area, is_inside_polygon
-from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
-from time import localtime, strftime
+#read in polyline boundary
+from anuga.shallow_water.data_manager import urs2sts
+base_name='polyline'
+import os
+if os.path.exists(base_name+'.sts'):
+        print 'sts boundary file already created.'
+else:
+    urs2sts(base_name,mean_stage=0.0,verbose=False)
+boundary=create_sts_boundary(base_name,verbose=True)
+urs_boundary=boundary.tolist()
 
-def convert_arcgis_latlon_list_to_utm(points):
-     #Used because arc gis produced csv files put lat lon in
-     #reverse order to those accpeted by convert_latlon_to_utm()
-     
-     from anuga.coordinate_transforms.geo_reference import Geo_reference
-     from anuga.coordinate_transforms.redfearn import redfearn
-     old_geo = Geo_reference() 
-     utm_points = []
-     for point in points:
-         zone, easting, northing = redfearn(float(point[1]),float(point[0]))
-         new_geo = Geo_reference(zone)
-         old_geo.reconcile_zones(new_geo) 
-         utm_points.append([easting,northing])
+#Read in large domain boundary
+dir='mesh_polygons'
+extent = 'extent.csv'
+bp = read_polygon(dir+sep+extent)
 
-     return utm_points, old_geo.get_zone()
+indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
+bounding_polygon=[]
+for i in indices[:-2]:
+    #disregard last two entries because they are not on boundary
+    bounding_polygon.append(urs_boundary[i])
+
+bounding_polygon.extend(bp[2:len(bp)])
 
 ###############################
@@ -45 +95 @@
 dir='mesh_polygons'
 extent = 'extent.csv'
-#extent = 'boundary_new_pts.csv'#old boundary
-# bounding polygon for study area
-bounding_polygon = read_polygon(dir+sep+extent)
 patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
 
@@ -54 +101 @@
 ###############################
 
-island_south = read_polygon(dir+sep+'south_island.csv')
-island_south2 = read_polygon(dir+sep+'south_island2.csv')
-island_north = read_polygon(dir+sep+'north_island.csv')
+contour20m = read_polygon(dir+sep+'20m.csv')
+contour20m[0][1]+=3000
+contour20m[-2][0]+=3000
+contour20m[-2][1]-=13000
+contour20m[-1][1]+=9000
+patong = read_polygon(dir+sep+'patong.csv')
 bay = read_polygon(dir+sep+'bay.csv')
-patong = read_polygon(dir+sep+'patong.csv')
-contour20m = read_polygon(dir+sep+'20m.csv')
 
-
-
-#max_extent_of MOST file
-#east  = 94.93325
-#west  = 80.0003
-#north = 10.0000
-#south = -4.9330
-
-#max_extent_of MOST file should be
-#east  = 104.0 
-#west  = 80.0
-#north = 10.0
-#south = -5.0
-
-#max_extent of old simulation is why does this work and not below
-#7.40911081272 8.71484358635 99.1683687224 97.513856322
-
-#max extent of clipped MOST file
+plot=False
+if plot:
+    from pylab import plot,show,hold
+    bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only
+    plot(bounding_polygon[:,0],bounding_polygon[:,1],'o')
+    plot(bounding_polygon[:,0],bounding_polygon[:,1])
+    contour20m = ensure_numeric(contour20m)
+    plot(contour20m[:,0],contour20m[:,1])
+    patong=ensure_numeric(patong)
+    plot(patong[:,0],patong[:,1])
+    show()
 
 east  = 97.7
@@ -84 +125 @@
 north = 9.0
 south = 5.9
-     
-     
+
+check_ferret_extent=False
+if check_ferret_extent:
+    east  = east
+    west  = west
+    north = north
+    south = south
+    
+    from anuga.coordinate_transforms.redfearn import redfearn
+    zone0, easting0, northing0 = redfearn(north,east)
+    print zone0,easting0,northing0
+    zone1, easting1, northing1 = redfearn(south,east)
+    print zone1,easting1,northing1
+    zone2, easting2, northing2 = redfearn(south,west)
+    print zone2,easting2,northing2
+    zone3, easting3, northing3 = redfearn(north,west)
+    print zone3,easting3,northing3
+    
+    assert zone0==zone1==zone2==zone3
+    
+    ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]]
+    points_outside=[]
+    from anuga.utilities.polygon import is_inside_polygon
+    for i,point in enumerate(bounding_polygon):
+        if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False):
+            print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1])
+            points_outside.append(point)
+
+    from pylab import plot,show
+    ferret_polygon=ensure_numeric(ferret_polygon)
+    points_outside=ensure_numeric(points_outside)
+    plot(ferret_polygon[:,0],ferret_polygon[:,1])
+    plot(points_outside[:,0],points_outside[:,1],'o')
+    show()
+    print ferret_polygon[0]
+    print ferret_polygon[1]
+    print ferret_polygon[2]
+    print ferret_polygon[3]
+    print points_outside
+