Changeset 6343

Timestamp:

Feb 16, 2009, 9:00:59 AM (16 years ago)

Author:

kristy

Message:

Added build_boundary

File:

• anuga_work/production/busselton/standardised_version/run_model.py (modified) (3 diffs)

anuga_work/production/busselton/standardised_version/run_model.py

@@ -22 +22 @@
 # Standard modules
 import os
+import os.path
 import time
+from time import localtime, strftime, gmtime
 
 # Related major packages
+from Scientific.IO.NetCDF import NetCDFFile
+import Numeric as num
+
 from anuga.interface import create_domain_from_regions
 from anuga.interface import Transmissive_stage_zero_momentum_boundary
@@ -32 +37 @@
 from anuga.interface import create_sts_boundary
 from anuga.interface import csv2building_polygons
+<<<<<<< .mine
+=======
 from file_length import file_length
 
+>>>>>>> .r6341
 from anuga.shallow_water.data_manager import start_screen_catcher
 from anuga.shallow_water.data_manager import copy_code_files
+from anuga.shallow_water.data_manager import urs2sts
 from anuga.utilities.polygon import read_polygon, Polygon_function
-    
+
 # Application specific imports
 from setup_model import project
-import build_urs_boundary as bub
-
+
+
+#-------------------------------------------------------------------------------
+# Get gauges (timeseries of index points)
+#-------------------------------------------------------------------------------
+def get_sts_gauge_data(filename, verbose=False):
+    fid = NetCDFFile(filename+'.sts', 'r')      #Open existing file for read
+    permutation = fid.variables['permutation'][:]
+    x = fid.variables['x'][:] + fid.xllcorner   #x-coordinates of vertices
+    y = fid.variables['y'][:] + fid.yllcorner   #y-coordinates of vertices
+    points = num.transpose(num.asarray([x.tolist(), y.tolist()]))
+    time = fid.variables['time'][:] + fid.starttime
+    elevation = fid.variables['elevation'][:]
+       
+    basename = 'sts_gauge'
+    quantity_names = ['stage', 'xmomentum', 'ymomentum']
+    quantities = {}
+    for i, name in enumerate(quantity_names):
+        quantities[name] = fid.variables[name][:]
+
+    #---------------------------------------------------------------------------
+    # Get maximum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+
+    maxname = 'max_sts_stage.csv'
+    fid_max = open(os.path.join(project.event_sts, maxname), 'w')
+    fid_max.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage)))
+      
+    #---------------------------------------------------------------------------
+    # Get minimum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+
+    minname = 'min_sts_stage.csv'
+    fid_min = open(os.path.join(project.event_sts, minname), 'w')
+    fid_min.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)))
+
+        out_file = os.path.join(project.event_sts,
+                                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]))
+
+        fid_sts.close()      
+
+    fid.close()
+
+    return quantities,elevation,time
+
+
+##
+# @brief Build boundary STS files from one or more MUX files.
+# @param mux_dir Directory containing the MUX files.
+# @param event_file Path to meta-file containing mux files+weight data.
+# @param output_dir Directory to write STS data to.
+# @note 'event_file' is produced by EventSelection.
+def build_urs_boundary(mux_dir, event_file, output_dir):
+    '''Build a boundary STS file from a set of MUX files.'''
+
+    print 'build_urs_boundary: mux_dir=%s' % mux_dir
+    print 'build_urs_boundary: event_file=%s' % event_file
+    print 'build_urs_boundary: output_dir=%s' % output_dir
+
+    # if we are using an EventSelection multi-mux file
+    if project.multi_mux:    
+        # get the mux+weight data from the event file
+        mux_event_file = os.path.join(mux_dir, event_file)
+        try:
+            fd = open(mux_event_file, 'r')
+            mux_data = fd.readlines()
+            fd.close()
+        except IOError, e:
+            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
+            raise Exception, msg
+        except:
+            raise
+
+        # first line of file is # filenames+weight in rest of file
+        num_lines = int(mux_data[0].strip())
+        mux_data = mux_data[1:]
+        print 'number of sources %d' % num_lines
+
+        # quick sanity check on input mux meta-file
+        if num_lines != len(mux_data):
+            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
+                   % (event_file, len(mux_data), num_lines))
+            raise Exception, msg
+
+        # Create filename and weights lists.
+        # Must chop GRD filename just after '*.grd'.
+        mux_filenames = []
+        for line in mux_data:
+            muxname = line.strip().split()[0]
+            split_index = muxname.index('.grd')
+            muxname = muxname[:split_index+len('.grd')]
+            muxname = os.path.join(mux_dir, muxname)
+            mux_filenames.append(muxname)
+        
+        mux_weights = [float(line.strip().split()[1]) for line in mux_data]
+
+        # Call legacy function to create STS file.
+        # This should be replaced in the future.
+        print 'reading', project.urs_order
+        print 'creating STS file'
+        print 'mux_filenames=%s' % str(mux_filenames)
+        print 'basename_out=%s' % str(output_dir)
+        print 'ordering_filename=%s' % str(project.urs_order)
+        print 'weights=%s' % str(mux_weights)
+        print 'mean_stage=%s' % str(project.tide)
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=project.urs_order,
+                weights=mux_weights,
+                mean_stage=project.tide,
+                verbose=True)
+    else:                           # a single mux stem file
+        urs_filenames = [os.path.join(mux_dir, event_file)]
+
+        weight_factor = 1.0
+        weights = weight_factor*num.ones(len(urs_filenames), num.float)
+            
+        order_filename = os.path.join(project.order_filename_dir)
+
+        print 'reading', order_filename
+        # Create ordered sts file
+        print 'creating sts file'
+        urs2sts(urs_filenames,
+                basename_out=os.path.join(project.boundaries_dir,
+                                          project.scenario_name),
+                ordering_filename=order_filename,
+                weights=weights,
+                mean_stage=project.tide,
+                verbose=True)
+
+    # report on stuff so far
+    quantities, elevation, time = get_sts_gauge_data(project.event_folder,
+                                                     verbose=False)
+    print len(elevation), len(quantities['stage'][0,:])
 
 #-------------------------------------------------------------------------------
@@ -62 +226 @@
 
 # Create the STS file
+<<<<<<< .mine
+build_urs_boundary(project.mux_data_folder,
+                   project.mux_input,
+                   os.path.join(project.event_folder,
+                                project.scenario_name))
+=======
 print 'project.mux_data_folder=%s' % project.mux_data_folder
 bub.build_urs_boundary(project.mux_input_filename,
                        os.path.join(project.event_folder,
                                     project.scenario_name))
+>>>>>>> .r6341
 
 # Read in boundary from ordered sts file