Changeset 7794 for branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py

Timestamp:

Jun 7, 2010, 9:38:23 AM (13 years ago)

Author:

sexton

Message:

updates to bunbury scripts

File:

• branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py (modified) (6 diffs)

branches/anuga_1_1/anuga_work/production/bunbury_storm_surge_2009/build_gems_boundary.py

@@ -13 +13 @@
 import sys
 # sys.path.append('W:\\sandpits\\lfountain\\anuga_work\\production\\mandurah_storm_surge_2009\\gems_comparison')
-sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
-import project_GEMS_grid
+#sys.path.append('/nas/gemd/georisk_models/inundation/sandpits/lfountain/anuga_work/production/mandurah_storm_surge_2009/gems_comparison')
+import project
 from Scientific.IO.NetCDF import NetCDFFile
 from anuga.coordinate_transforms.geo_reference import Geo_reference, write_NetCDF_georeference
@@ -23 +23 @@
 state = 'western_australia'
 scenario_folder = 'bunbury_storm_surge_scenario_2009'
-event = 'case_1' # Baseline TC Alby event
+event = '20100527_gcom_12min' # Baseline TC Alby event
 
 ENV_INUNDATIONHOME = 'INUNDATIONHOME'
@@ -35 +35 @@
 # Define boundary points to extract from master STS file
 #--------------------------------------------------------
-gems_boundary = numpy.array(read_polygon(project_GEMS_grid.gems_order))
+gems_boundary = numpy.array(read_polygon(project.gems_order))
+#gems_boundary = numpy.array(read_polygon(project.gauges))
 number_of_points = gems_boundary.shape[0]
+print 'hello', number_of_points
+print 'points', gems_boundary
 
 #---------------------------------------------------------
@@ -42 +45 @@
 #---------------------------------------------------------
 print "Reading master NetCDF file"
-infile = NetCDFFile(os.path.join(event_folder, event + '_master.sts'), 'r')
+infile = NetCDFFile(os.path.join(event_folder, event + '_master_2_1.sts'), 'rl') #rl for large netcdf files
 
 grid_size = infile.grid_size
@@ -69 +72 @@
 origin = numpy.array([x_origin[0], y_origin[0]])
 index_array_2D = ((gems_boundary - origin)/grid_size).astype(int)
-assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols
+assert index_array_2D[:,0].min() >= 0 and index_array_2D[:,0].max() <= ncols # These test that the boundary points lie within the data area
 assert index_array_2D[:,1].min() >= 0 and index_array_2D[:,1].max() <= nrows
 
@@ -135 +138 @@
     s[i] = stage[i,:]
     xm[i] = xmomentum[i,:]
-    ym[i] = ymomentum[i,:]
- 
+    ym[i] = ymomentum[i,:] 
+
+#---------------------------------------------------------------------------
+# Get timeseries values for wave height and components of momentum
+#---------------------------------------------------------------------------
+
+basename = 'sts_gauge_tide'
+for i in xrange(number_of_points):
+    out_file = os.path.join(event_folder,
+                            basename+'_'+str(i)+'.csv')
+    fid_sts = open(out_file, 'w')
+    fid_sts.write('time, stage, xmomentum, ymomentum, elevation \n')
+    for j in xrange(number_of_timesteps):
+        fid_sts.write('%.6f, %.6f, %.6f, %.6f, %.6f\n'
+                    % (time[j], stage[j,i], xmomentum[j,i], ymomentum[j,i], elevation[i]))
+
+    fid_sts.close() 
+
 origin = Geo_reference(refzone, x_origin, y_origin)
 geo_ref = write_NetCDF_georeference(origin, fid)
 
 fid.variables['x'][:] = x 
-fid.variables['y'][:] = y 
+fid.variables['y'][:] = y  
 
 fid.close()
+
+sys.exit()
+
+print 'get timeseries'
+#-------------------------------------------------------------------------------
+# Get gauges (timeseries of index points)
+#-------------------------------------------------------------------------------
+
+basename = 'sts_gauge'
+quantity_names = ['stage', 'xmomentum', 'ymomentum']
+quantities = {}
+for i, name in enumerate(quantity_names):
+    quantities[name] = fid.variables[name][:]
+
+
+for j in range(len(x)):
+    index = j # permutation[j]
+    stage = quantities['stage'][:,j]
+    #print 'hello', j, stage
+    xmomentum = quantities['xmomentum'][:,j]
+    ymomentum = quantities['ymomentum'][:,j]
+
+    out_file = os.path.join(event_folder,
+                            basename+'_'+str(index)+'.csv')
+    fid_sts = open(out_file, 'w')
+    fid_sts.write('time, stage, xmomentum, ymomentum \n')
+
+    #-----------------------------------------------------------------------
+    # End of the get gauges
+    #-----------------------------------------------------------------------
+    for k in range(len(time)-1):
+        fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
+                      % (time[k], stage[k], xmomentum[k], ymomentum[k]))
+
+