Changeset 5465 for anuga_work/development/boxingday08/project.py

Timestamp:

Jul 4, 2008, 1:56:15 PM (16 years ago)

Author:

ole

Message:

updated boxingday scenario to use ordered sts file

File:

• anuga_work/development/boxingday08/project.py (modified) (7 diffs)

anuga_work/development/boxingday08/project.py

@@ -2 +2 @@
 from Scientific.IO.NetCDF import NetCDFFile
 from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
-    compress,zeros,fabs,allclose
+    compress,zeros,fabs,allclose,ones
 from anuga.utilities.polygon import inside_polygon,read_polygon
 from os import sep
 from time import localtime, strftime, gmtime
-
-def create_sts_boundary_order_file(order_filename,stsfilename,verbose=False):
+from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
+import os
+
+def create_sts_boundary_order_file(urs_filename,order_filename,sts_filename,verbose=False):
     """
     Note This is not a robust algorithm. Be Careful when applying that 
@@ -13 +15 @@
     follows the path it is supposed too.
     """
-    fid = NetCDFFile(stsfilename+'.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()]))
+    from anuga.shallow_water.data_manager import read_mux2_py
+    times, y, x, elevation, quantity, starttime = read_mux2_py([urs_filename],weights=ones(1,Float))
     
     #find indices of the southernmost gauges
@@ -60 +60 @@
               
     #So far i cannot tell if the wrong point is closer.
+    
+    boundary=ensure_numeric(boundary)
+    from anuga.coordinate_transforms import redfearn,convert_from_latlon_to_utm
+    boundary_utm,zone=convert_from_latlon_to_utm(longitudes=list(boundary[:,0]),latitudes=list(boundary[:,1]))
+
     d=","
     fid=open(order_filename,'w')
-    header="index,longitude,latitude\n"
+    header='index, longitude, latitude\n'
     fid.write(header)
     line=str(west_index)+d+str(x[west_index])+d+\
@@ -72 +77 @@
             fid.write(line)
     fid.close()
-    return boundary
-            
-
-#read in polyline boundary
-from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
-urs_name='polyline'
-#base_name=urs_name
-tide=0.35
-base_name='tide_polyline'
-#base_name='most_polyline'
-import os
-if os.path.exists(base_name+'.sts'):
-    print 'sts boundary file already created.'
-else:
-    print 'creatin sts file'
-    urs2sts(urs_name,base_name,mean_stage=tide,verbose=False)
-
-order_filename='sts_order.txt'
-boundary=create_sts_boundary_order_file(order_filename,base_name,verbose=False)
-urs_boundary=create_sts_boundary(order_filename,base_name,lat_long=False)
-assert allclose(boundary.tolist(),urs_boundary)
-
-#Read in large domain boundary
-dir='mesh_polygons'
-extent = 'extent.csv'
-bp = read_polygon(dir+sep+extent)
-
-#Now order of boundary is correct clip sts boundary to small region for
-#this particular scenario
-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])
-
-#Append clipped sts boundary to large domain boundary
-bounding_polygon.extend(bp[2:len(bp)])
-
-###############################
-# Domain definitions
-###############################
+    return boundary,boundary_utm
+
+######################
+# Domain definitions #
+######################
 
 # Bathymetry and topography filenames
@@ -129 +98 @@
 extent = 'extent.csv'
 patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
+    
+tide=0.35
+
+######################################
+# Create urs boundary from mux2files #
+######################################
+urs_filename='polyline'
+base_name='tide_polyline'
+order_filename='boxingday_boundary_order.txt'
+
+# Create ordering file
+# Usually this would be provided before hand
+boundary,boundary_utm=create_sts_boundary_order_file(urs_filename+'_waveheight-z-mux2',order_filename,base_name,verbose=False)
+
+# Create ordered sts file
+if not os.path.exists(base_name+'.sts'):
+    print 'creating sts file'
+    urs2sts(urs_filename,basename_out=base_name,
+    ordering_filename=order_filename,
+    mean_stage=tide,
+    verbose=False)
+else:
+    print 'sts file exists'
+
+# Read in boundary from ordered sts file
+urs_boundary=create_sts_boundary(base_name)
+
+# Checks that the boundary read from the sts file is the same as the one
+# defined by the function create_sts_boundary_order_file
+assert allclose(boundary_utm,urs_boundary)
+
+#############################################################
+# Add urs boundary segment to user defined boundary polygon #
+#############################################################
+
+# Read in specific user defined boundary
+dir='mesh_polygons'
+extent = 'extent.csv'
+user_boundary_polygon = read_polygon(dir+sep+extent)
+
+# The following is specific to this scenario
+indices=inside_polygon(urs_boundary, user_boundary_polygon, 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])
+
+# Append clipped sts boundary to user defined domain boundary
+bounding_polygon.extend(user_boundary_polygon[2:len(user_boundary_polygon)])
 
 ###############################
-# Interior region definitions
+# Interior region definitions #
 ###############################
 
@@ -141 +159 @@
 patong = read_polygon(dir+sep+'patong.csv')
 bay = read_polygon(dir+sep+'bay.csv')
+
+##########################################
+# Plot urs boundary and internal regions #
+##########################################
 
 plot=False
@@ -156 +178 @@
     show()
 
+    
+#######################
+# For MOST SIMULATION #
+#######################
 east  = 97.7
 west  = 96.7