Changeset 3944

Timestamp:

Nov 8, 2006, 5:16:44 PM (16 years ago)

Author:

sexton

Message:

(i) update sww2timeseries so can handle gauges which don't fall within sww domain (ii) script for generating report for dampier based on sww files from parallel setup

Files:

• anuga_core/documentation/user_manual/anuga_user_manual.tex (modified) (1 diff)
• anuga_core/documentation/user_manual/examples/project.py (modified) (1 diff)
• anuga_core/documentation/user_manual/examples/runsydney.py (modified) (1 diff)
• anuga_core/source/anuga/abstract_2d_finite_volumes/util.py (modified) (14 diffs)
• anuga_work/development/cairns_2006/run_cairns.py (modified) (2 diffs)
• anuga_work/production/broome_2006/export_results.py (modified) (3 diffs)
• anuga_work/production/broome_2006/project.py (modified) (1 diff)
• anuga_work/production/dampier_2006/make_report_cipma.py (added)
• anuga_work/production/hobart_2006/export_results.py (modified) (2 diffs)
• anuga_work/production/hobart_2006/make_report.py (modified) (1 diff)
• anuga_work/production/hobart_2006/report/hobart_2006_report_MRT.tex (modified) (1 diff)

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -3219 +3219 @@
 http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
 
-
+\bibitem[grid250]{grid250}
+Australian Bathymetry and Topography Grid, June 2005.
+Webster, M.A. and Petkovic, P.
+Geoscience Australia Record 2005/12. ISBN: 1-920871-46-2.\\
+http://www.ga.gov.au/meta/ANZCW0703008022.html
 
 

anuga_core/documentation/user_manual/examples/project.py

@@ -45 +45 @@
 # demo poly 
 j0 = [385000, 6280000]
-j1 = [360000, 6272500]
-j2 = [335000, 6272500]
+j1 = [360000, 6273000]
+j2 = [335000, 6273000]
 j3 = [330000, 6265000]
 j31 = [325000, 6260000]
 j4 = [316000, 6260000]
-j5 = [316000, 6247000]
-j6 = [350000, 6247000]
+j5 = [316000, 6246750]
+j6 = [350000, 6246750]
 j7 = [385000, 6238000]
 
 demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
 
-from utilities.polygon import read_polygon
-polygonptsfile4 = polygondir + 'poly1'
-polygonptsfile0 = polygondir + 'poly2'
-polygonptsfile1 = polygondir + 'poly3'
-polygonptsfile2 = polygondir + 'poly4'
-polygonptsfile3 = polygondir + 'poly5'
-northern_polygon = read_polygon(polygonptsfile0 + '.csv')
-manly_polygon = read_polygon(polygonptsfile1 + '.csv')
-harbour_polygon = read_polygon(polygonptsfile2 + '.csv')
-southern_polygon = read_polygon(polygonptsfile3 + '.csv')
-top_polygon = read_polygon(polygonptsfile4 + '.csv')
+from anuga.utilities.polygon import read_polygon, plot_polygons
+#polygonptsfile4 = polygondir + 'poly1'
+#polygonptsfile0 = polygondir + 'poly2'
+#polygonptsfile1 = polygondir + 'poly3'
+#polygonptsfile2 = polygondir + 'poly4'
+#polygonptsfile3 = polygondir + 'poly5'
+#northern_polygon = read_polygon(polygonptsfile0 + '.csv')
+#manly_polygon = read_polygon(polygonptsfile1 + '.csv')
+#harbour_polygon = read_polygon(polygonptsfile2 + '.csv')
+#southern_polygon = read_polygon(polygonptsfile3 + '.csv')
+#top_polygon = read_polygon(polygonptsfile4 + '.csv')
+coastal_polygon = read_polygon(polygondir+'coastal'+'.csv')
+shallow_polygon = read_polygon(polygondir+'shallow'+'.csv')
+
+#plot_polygons([demopoly,northern_polygon,manly_polygon,harbour_polygon,southern_polygon,top_polygon],'model_setup',verbose=False)
+plot_polygons([demopoly,coastal_polygon,shallow_polygon],'new_model_setup',verbose=False)
 
 slump_origin = [372500.0, 6255000.0] #Absolute UTM

anuga_core/documentation/user_manual/examples/runsydney.py

@@ -59 +59 @@
 # Interior regions and mesh resolutions
 interior_res = 5000
-interior_regions = [[project.northern_polygon, interior_res],
-                    [project.harbour_polygon, interior_res],
-                    [project.southern_polygon, interior_res],
-                    [project.manly_polygon, interior_res],
-                    [project.top_polygon, interior_res]]
-
+shallow_res = 15000
+##interior_regions = [[project.northern_polygon, interior_res],
+##                    [project.harbour_polygon, interior_res],
+##                    [project.southern_polygon, interior_res],
+##                    [project.manly_polygon, interior_res],
+##                    [project.top_polygon, interior_res]]
+interior_regions = [[project.coastal_polygon, interior_res],
+                    [project.shallow_polygon, shallow_res]]
 
 create_mesh_from_regions(project.demopoly,

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -681 +681 @@
         
     assert type(gauge_filename) == type(''),\
-               'Gauge filename must be a string'
+           'Gauge filename must be a string'
     
     try:
@@ -693 +693 @@
         report = False
         
-    
     if plot_quantity is None:
         plot_quantity = ['depth', 'speed', 'bearing']
@@ -733 +732 @@
                           verbose = True,
                           use_cache = True)
- 
+
+        # determine which gauges are contained in sww file
+        count = 0
+        gauge_index = []
+        print 'swwfile', swwfile
+        for k, g in enumerate(gauges):
+            if f(0.0, point_id = k)[2] > 1.0e6:
+                count += 1
+                if count == 1: print 'Gauges not contained here:'
+                print locations[k]
+            else:
+                gauge_index.append(k)
+
+        if len(gauge_index) > 0:
+            print 'Gauges contained here: \n',
+        else:
+            print 'No gauges contained here. \n'
+        for i in range(len(gauge_index)):
+             print locations[gauge_index[i]]
+             
         index = swwfile.rfind(sep)
         file_loc.append(swwfile[:index+1])
@@ -759 +777 @@
             raise Exception, msg
 
-    if verbose: print 'Inputs OK - going to generate figures'
-
-    return generate_figures(plot_quantity, file_loc, report, reportname, surface,
-                            leg_label, f_list, gauges, locations, elev, production_dirs,
-                            time_min, time_max, title_on, label_id, verbose)
+    if verbose and len(gauge_index) > 0: print 'Inputs OK - going to generate figures'
+
+    if len(gauge_index) <> 0:
+        texfile, elev_output = generate_figures(plot_quantity, file_loc, report, reportname, surface,
+                                                leg_label, f_list, gauges, locations, elev, gauge_index,
+                                                production_dirs, time_min, time_max, title_on, label_id, verbose)
+    else:
+        texfile = ''
+        elev_output = []
+
+    return texfile, elev_output
                          
 #Fixme - Use geospatial to read this file - it's an xya file
 #Need to include other information into this filename, so xya + Name - required for report
 def get_gauges_from_file(filename):
+    """ Read in gauge information from file
+    """
     from os import sep, getcwd, access, F_OK, mkdir
     fid = open(filename)
@@ -803 +829 @@
 
 def check_list(quantity):
-
+    """ Check that input quantities in quantity list are possible
+    """
     all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
                     'ymomentum', 'speed', 'bearing', 'elevation']
@@ -817 +844 @@
 
 def calc_bearing(uh, vh):
-
+    """ Calculate velocity bearing from North
+    """
     from math import atan, degrees
     
@@ -836 +864 @@
 
 def generate_figures(plot_quantity, file_loc, report, reportname, surface,
-                     leg_label, f_list,
-                     gauges, locations, elev, production_dirs, time_min, time_max,
-                     title_on, label_id, verbose):
-
+                     leg_label, f_list, gauges, locations, elev, gauge_index,
+                     production_dirs, time_min, time_max, title_on, label_id,
+                     verbose):
+    """ Generate figures based on required quantities and gauges for each sww file
+    """
     from math import sqrt, atan, degrees
     from Numeric import ones, allclose, zeros, Float, ravel
@@ -849 +878 @@
     if surface is True:
         import mpl3d.mplot3d as p3
-
+        
     if report == True:    
         texdir = getcwd()+sep+'report'+sep
@@ -904 +933 @@
         fid_compare = open(comparefile, 'w')
         ##### loop over each gauge #####
-        for k, g in enumerate(gauges):
+        for k in gauge_index:
+            g = gauges[k]
             if verbose: print 'Gauge %d of %d' %(k, len(gauges))
             min_stage = 10
@@ -992 +1022 @@
         savefig('profilefig')
                 
-    #stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
-    stage_axis = axis([time_min/60.0, time_max/60.0, -3.0, 3.0])    
+    stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
+    #stage_axis = axis([time_min/60.0, time_max/60.0, -3.0, 3.0])    
     vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
     mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1])  
@@ -1002 +1032 @@
     elev_output = []
     if len(label_id) > 1: graphname_report = []
-    for k, g in enumerate(gauges):
+    for k in gauge_index:
+        g = gauges[k]
         count1 = 0
         if report == True and len(label_id) > 1:
@@ -1008 +1039 @@
             fid.write(s)
         if len(label_id) > 1: graphname_report = []
-        #### generate figures for each gauge ###
+        #### generate figures for each gauge ####
         for j, f in enumerate(f_list):
             ion()
@@ -1109 +1140 @@
                         word_quantity += plot_quantity[i]
                     
-                word_quantity += ' and ' + plot_quantity[nn-1]                
+                word_quantity += ' and ' + plot_quantity[nn-1]
                 caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' '))
                 if elev[k] == 0.0:
@@ -1164 +1195 @@
     return texfile2, elev_output
 
+
+def gauge_in_sww(swwfile,
+                 gauge_filename,
+                 verbose = False):
+    
+    """ Determine which gauges fall in sww file.
+
+    Input variables:
+
+    swwfiles        - sww files 
+    
+    gauge_filename  - name of file containing gauge data
+                        - easting, northing, name , elevation
+                    - OR (this is not yet done)
+                        - structure which can be converted to a Numeric array,
+                          such as a geospatial data object
+                      
+    Output:
+    
+    - Modified gauge_filename containing gauges which fall is given sww file.
+      Name = gauges_filename_mod 
+    Other important information:
+                          
+    """
+
+    
+    k = _gauge_in_sww(swwfile,
+                      gauge_filename,
+                      verbose)
+
+    return k
+
+def _gauge_in_sww(swwfile, 
+                  gauge_filename,
+                  verbose = False):   
+
+    try:
+        fid = open(swwfile)
+    except Exception, e:
+        msg = 'File "%s" could not be opened: Error="%s"'\
+              %(swwfile, e)
+        raise msg
+        
+    assert type(gauge_filename) == type(''),\
+           'Gauge filename must be a string'
+    
+    try:
+        fid = open(gauge_filename)
+    except Exception, e:
+        msg = 'File "%s" could not be opened: Error="%s"'\
+                  %(gauge_filename, e)
+        raise msg
+
+    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+    f = file_function(swwfile,
+                      quantities = sww_quantity,
+                      #interpolation_points = gauges,
+                      verbose = True,
+                      use_cache = True)
+    print dir(f)
+    
+    #from bstract_2d_finite_volumes.neighbour_mesh import get_boundary_polygon
+    from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
+    #from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
+    #from anuga.fit_interpolate.interpolate import Interpolation_function
+    from anuga.fit_interpolate.interpolate import Interpolate
+
+    test = Mesh(f.vertex_coordinates)
+    print test.mesh.get_boundary_polygon()
+    
+    return check_gauge(gauge_filename,bounding_polygon)
+
+def check_gauge(gauge_filename,bounding_polygon):
+
+    from anuga.utilities.polygon import is_inside_polygon 
+    
+    filename_out = gauge_filename + '_mod'
+    fid_out = open(filename_out, 'w')
+    fid_out.write(s)
+
+    gauges, locations, elev = get_gauges_from_file(gauge_filename)
+    
+    for i, gauge in enumerate(gauges):
+        v = is_inside_polygon(gauge,bounding_polygon,verbose=False)
+        if v == True:
+            s = '%.2f, %.2f, %.2f, %s\n' %(gauges[0], gauges[1], elev[i], locations[i])
+    
+    return fid_out
+
 from os.path import basename
 

anuga_work/development/cairns_2006/run_cairns.py

@@ -9 +9 @@
 #-------------------------------
 import sys, os
-from anuga.pyvolution.shallow_water import Domain, Reflective_boundary,\
+from anuga.shallow_water import Domain, Reflective_boundary,\
      File_boundary, Transmissive_Momentum_Set_Stage_boundary,\
      Transmissive_boundary, Dirichlet_boundary
-from anuga.pyvolution.mesh_factory import rectangular_cross
-from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
+#from anuga.pyvolution.mesh_factory import rectangular_cross
+#from anuga.pmesh2domain import pmesh_to_domain_instance
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from Numeric import array, zeros, Float, allclose
 import cairns_project as project
@@ -41 +42 @@
 """
 
-domain = cache(pmesh_to_domain_instance,
-               (project.mesh_filename, Domain),
-               dependencies = [project.mesh_filename])
+#cache(pmesh_to_domain_instance,
+#      (project.mesh_filename, Domain),
+#      dependencies = [project.mesh_filename])
+domain = Domain(project.mesh_filename, use_cache = True, verbose = True)
 
 domain.check_integrity()

anuga_work/production/broome_2006/export_results.py

@@ -6 +6 @@
 from os import sep
 
-time_dir = '20061026_220716' 
+time_dir = '20061101_024119' 
 
 directory = project.outputdir
@@ -22 +22 @@
 
 if which_var == 2:  # Depth
-    outname = name + '_depth_bruny'
+    outname = name + '_depth'
     quantityname = 'stage-elevation'  #Depth
 
@@ -37 +37 @@
 sww2dem(name, basename_out = outname,
             quantity = quantityname,
-            cellsize = 25,      # would prefer this at 25
+            cellsize = 100,      # would prefer this at 25
             # define region for viz purposes
             easting_min = project.eastingmin,

anuga_work/production/broome_2006/project.py

@@ -138 +138 @@
 
 # exporting asc grid
-eastingmin = 406215.87
-eastingmax = 440208.78
-northingmin = 7983427.73
-northingmax = 8032834.52
+eastingmin = 412036.43
+eastingmax = 427184.37
+northingmin = 8007984
+northingmax = 8029830
 
         

anuga_work/production/hobart_2006/export_results.py

@@ -49 +49 @@
 ##            format = 'asc')
 
-# for bruny
+### for bruny
+##sww2dem(name, basename_out = outname,
+##            quantity = quantityname,
+##            cellsize = 25,      # would prefer this at 25
+##            # define region for viz purposes
+##            easting_min = project.eastingmin25_2,
+##            easting_max = project.eastingmax25_2,
+##            northing_min = project.northingmin25_2,
+##            northing_max = project.northingmax25_2,        
+##            reduction = max, #this is because we want max quantityname
+##            verbose = True,
+##            format = 'asc')
+
+# for fixed timestep
 sww2dem(name, basename_out = outname,
             quantity = quantityname,
-            cellsize = 25,      # would prefer this at 25
+            cellsize = 10,      # would prefer this at 25
             # define region for viz purposes
             easting_min = project.eastingmin25_2,
@@ -58 +71 @@
             northing_min = project.northingmin25_2,
             northing_max = project.northingmax25_2,        
-            reduction = max, #this is because we want max quantityname
+            #reduction = max, #this is because we want max quantityname
             verbose = True,
             format = 'asc')
+

anuga_work/production/hobart_2006/make_report.py

@@ -99 +99 @@
                                       #reportname = 'latexoutput',
                                       plot_quantity = ['stage', 'speed'],
-                                      surface = False,
+                                      surface = True,
                                       time_min = None,
                                       time_max = None,

anuga_work/production/hobart_2006/report/hobart_2006_report_MRT.tex

@@ -172 +172 @@
      \input{metadata}
 
-\pagebreak
+\clearpage
 
    \section{Time series}