Changeset 5402


Ignore:
Timestamp:
Jun 13, 2008, 10:37:14 AM (14 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:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/boxingday08/get_gauges.py

    r5393 r5402  
    11import os
    22from os import sep
    3 import create_boundary
    4 from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
     3#import create_boundary
     4import project
     5try:
     6    from pylab import plot,show,xlabel,ylabel,legend,title,savefig,hold
     7except:
     8    print 'Cannot import pylab plotting will not work. Csv files are still created'
    59from anuga.utilities.polygon import inside_polygon
    610from anuga.utilities.numerical_tools import ensure_numeric
     
    1620    elevation=fid.variables['elevation'][:]
    1721    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)
    2224    if verbose:
    2325        inside_points_array=ensure_numeric(inside_points)#for plotting only
    24         #plot(points[:,0],points[:,1],'o')
     26        plot(points[:,0],points[:,1],'o')
    2527        plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o')
     28        polygon=ensure_numeric(polygon)
    2629        plot(polygon[:,0],polygon[:,1])
    2730        show()
    28 
     31       
    2932    quantity_names=['stage','xmomentum','ymomentum']
    3033    quantities = {}
     
    3841    fid.close()
    3942
    40     return points,quantities,elevation,time
     43    return inside_points,quantities,elevation,time
    4144
    4245from anuga.shallow_water.data_manager import urs2sts
     
    4851    urs2sts(base_name,mean_stage=0.0,verbose=False)
    4952
    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
     54polygon=project.bounding_polygon
     55points,quantities,elevation,time=get_sts_gauge_data(base_name,polygon,verbose=False)
     56assert len(points)==len(elevation)==len(quantities['stage'][0,:])
    5357
    5458def 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',
    5660                                      quantity_name,
    5761                                      str(gauge_id))
     
    6064    title('Gauge %d at (%f,%f)' %(gauge_id,pos[0],pos[1]))
    6165    xlabel('Time (minutes)')
    62     ylabel('Stage (m)')
     66    ylabel(quantity_name)
    6367    legend()
    6468    print 'saving figure here %s' %figname
     
    6771points=points[::5]
    6872elevation=elevation[::5]
    69 print points.shape
    7073gauge_filename='urs_gauges'
    7174d=','
     
    9194       
    9295    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]},
    9497                          base_name='gauge_',
    9598                          plot_numbers='',#['1'],
    96                           quantities=['stage'],
     99                          quantities=['stage','xmomentum','ymomentum'],
    97100                          extra_plot_name='',
    98101                          assess_all_csv_files=True,
     
    100103                          verbose=True)
    101104
    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])
     105sww_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)
     109for 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 """
     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
     5from anuga.utilities.polygon import inside_polygon,read_polygon
     6from os import sep
     7from time import localtime, strftime, gmtime
     8def 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
    455
    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
     57from anuga.shallow_water.data_manager import urs2sts
     58base_name='polyline'
     59import os
     60if os.path.exists(base_name+'.sts'):
     61        print 'sts boundary file already created.'
     62else:
     63    urs2sts(base_name,mean_stage=0.0,verbose=False)
     64boundary=create_sts_boundary(base_name,verbose=True)
     65urs_boundary=boundary.tolist()
    1366
    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
     68dir='mesh_polygons'
     69extent = 'extent.csv'
     70bp = read_polygon(dir+sep+extent)
    2771
    28      return utm_points, old_geo.get_zone()
     72indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
     73bounding_polygon=[]
     74for i in indices[:-2]:
     75    #disregard last two entries because they are not on boundary
     76    bounding_polygon.append(urs_boundary[i])
     77
     78bounding_polygon.extend(bp[2:len(bp)])
    2979
    3080###############################
     
    4595dir='mesh_polygons'
    4696extent = 'extent.csv'
    47 #extent = 'boundary_new_pts.csv'#old boundary
    48 # bounding polygon for study area
    49 bounding_polygon = read_polygon(dir+sep+extent)
    5097patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
    5198
     
    54101###############################
    55102
    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')
     103contour20m = read_polygon(dir+sep+'20m.csv')
     104contour20m[0][1]+=3000
     105contour20m[-2][0]+=3000
     106contour20m[-2][1]-=13000
     107contour20m[-1][1]+=9000
     108patong = read_polygon(dir+sep+'patong.csv')
    59109bay = read_polygon(dir+sep+'bay.csv')
    60 patong = read_polygon(dir+sep+'patong.csv')
    61 contour20m = read_polygon(dir+sep+'20m.csv')
    62110
    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
     111plot=False
     112if 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()
    81122
    82123east  = 97.7
     
    84125north = 9.0
    85126south = 5.9
    86      
    87      
     127
     128check_ferret_extent=False
     129if 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.