Changeset 2648


Ignore:
Timestamp:
Mar 31, 2006, 6:21:26 PM (18 years ago)
Author:
steve
Message:

Looking at island.py example. I changed the order of the implicit and explicit update (in quantity_ext.c) which I think improved the situation a little. But this has left a few fails in the test_all suite which we need to fix up. Mainly small differences in results.

Files:
11 added
16 edited

Legend:

Unmodified
Added
Removed
  • inundation/examples/island.py

    r2638 r2648  
    3333                                           'top': [2],
    3434                                           'left': [3]},
    35                           maximum_triangle_area = 20,
     35                          maximum_triangle_area = 5,
    3636                          filename = 'island.msh'
    37                           , interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
     37              , interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
    3838                          )
    3939
     
    7777    return z
    7878
    79 domain.set_quantity('friction', 0.1)  #Honky dory
    80 #domain.set_quantity('friction', 2)     #Creep
     79#domain.set_quantity('friction', 0.1)  #Honky dory
     80domain.set_quantity('friction', 1000)     #Creep
    8181domain.set_quantity('elevation', island)
    8282domain.set_quantity('stage', 1)
     83domain.max_timestep = 0.01
    8384
    8485
  • inundation/pyvolution/config.py

    r2566 r2648  
    1212
    1313min_timestep = 1.0e-6 #Should be computed based on geometry
    14 max_timestep = 1000
     14max_timestep = 1.0e+3
    1515#This is how:
    1616#Define maximal possible speed in open water v_max, e.g. 500m/s (soundspeed?)
  • inundation/pyvolution/data_manager.py

    r2644 r2648  
    10491049    if zone is None:
    10501050        assert xllcorner is None, 'xllcorner must be None'
    1051         assert yllcorner is None, 'yllcorner must be None'       
     1051        assert yllcorner is None, 'yllcorner must be None'
    10521052        Geo_reference().write_NetCDF(outfile)
    10531053    else:
    1054         Geo_reference(zone, xllcorner, yllcorner).write_NetCDF(outfile)       
    1055        
     1054        Geo_reference(zone, xllcorner, yllcorner).write_NetCDF(outfile)
     1055
    10561056
    10571057
     
    10751075
    10761076
    1077 def dem2pts(basename_in, basename_out=None, 
     1077def dem2pts(basename_in, basename_out=None,
    10781078            easting_min=None, easting_max=None,
    10791079            northing_min=None, northing_max=None,
     
    10991099
    11001100
    1101     kwargs = {'basename_out': basename_out, 
     1101    kwargs = {'basename_out': basename_out,
    11021102              'easting_min': easting_min,
    11031103              'easting_max': easting_max,
     
    11141114    else:
    11151115        result = apply(_dem2pts, [basename_in], kwargs)
    1116        
     1116
    11171117    return result
    11181118
     
    11231123    """Read Digitial Elevation model from the following NetCDF format (.dem)
    11241124
    1125     Internal function. See public function dem2pts for details.   
     1125    Internal function. See public function dem2pts for details.
    11261126    """
    11271127
     
    12191219
    12201220                if k == 0:
    1221                     i1_0 = i 
     1221                    i1_0 = i
    12221222                    j1_0 = j
    1223                 k += 1   
     1223                k += 1
    12241224
    12251225    index1 = j1_0
    12261226    index2 = thisj
    1227    
     1227
    12281228    # dimension definitions
    12291229    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
    12301230    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
    1231    
    1232     clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0) 
     1231
     1232    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
    12331233    nopoints = clippednopoints-nn
    12341234
    12351235    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
    1236    
     1236
    12371237    if verbose and nn > 0:
    12381238        print 'There are %d values in the elevation' %totalnopoints
    12391239        print 'There are %d values in the clipped elevation' %clippednopoints
    12401240        print 'There are %d NODATA_values in the clipped elevation' %nn
    1241        
     1241
    12421242    outfile.createDimension('number_of_points', nopoints)
    12431243    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
     
    12731273
    12741274        local_index = 0
    1275        
     1275
    12761276        y = (nrows-i-1)*cellsize + yllcorner
    12771277        #for j in range(ncols):
     
    12861286                global_index += 1
    12871287                local_index += 1
    1288                
     1288
    12891289        upper_index = global_index
    12901290
    1291         if upper_index == lower_index + newcols: 
     1291        if upper_index == lower_index + newcols:
    12921292            points[lower_index:upper_index, :] = tpoints
    12931293            elevation[lower_index:upper_index] = telev
     
    19141914                                  use_cache = False,
    19151915                                  verbose = False):
    1916     """Read Digitial Elevation model from the following ASCII format (.asc)   
     1916    """Read Digitial Elevation model from the following ASCII format (.asc)
    19171917
    19181918    Example:
     
    19581958    else:
    19591959        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
    1960        
     1960
    19611961    return result
    19621962
    1963    
    1964    
     1963
     1964
    19651965
    19661966
     
    21612161    if maxlon != None:
    21622162        assert -180 < maxlon < 180 , msg
    2163  
     2163
    21642164
    21652165    #Get NetCDF data
     
    21682168    file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
    21692169    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
    2170     file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)
     2170    file_e = NetCDFFile(basename_in + '_e.nc', 'r')  #Elevation (z) (m)
    21712171
    21722172    if basename_out is None:
     
    26052605    from Numeric import array
    26062606    from config import time_format
    2607     from utilities.numerical_tools import ensure_numeric   
     2607    from utilities.numerical_tools import ensure_numeric
    26082608
    26092609
     
    28232823
    28242824
    2825     try:   
     2825    try:
    28262826        domain = Domain(coordinates, volumes, boundary)
    28272827    except AssertionError, e:
     
    28872887            X = resize(X,(X.shape[0]/3,3))
    28882888        domain.set_quantity(quantity,X)
    2889        
     2889
    28902890    fid.close()
    28912891    return domain
     
    35923592    quantities which can be derived from those, i.e.
    35933593    ['depth', 'xmomentum', 'ymomentum', 'momentum', 'velocity', 'bearing'].
    3594    
     3594
    35953595    Momentum is the absolute momentum, sqrt(xmomentum^2 + ymomentum^2).
    35963596    Note, units of momentum are m^2/s and depth is m.
     
    36113611
    36123612    # extract gauge locations from gauge file
    3613    
     3613
    36143614    # extract all quantities from swwfile (f = file_function)
    36153615    if time_min is None: time_min = min(f.T)
    36163616    if time_max is None: time_max = max(f.T)
    3617    
     3617
    36183618    # loop through appropriate range of time
    36193619    # plot prescribed quantities and export data if requested
    3620    
     3620
    36213621    #if gauge_data_outname is None:
    36223622
  • inundation/pyvolution/domain.py

    r2633 r2648  
    7878        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
    7979        #Or maybe get rid of order altogether and use beta_w and beta_h
    80         self.set_default_order(1)       
     80        self.set_default_order(1)
    8181        #self.default_order = 1
    8282        #self.order = self.default_order
     
    115115        if mesh_filename is not None:
    116116            # If the mesh file passed any quantity values
    117             # , initialise with these values. 
     117            # , initialise with these values.
    118118            self.set_quantity_vertices_dict(vertex_quantity_dict)
    119119
    120        
     120
    121121
    122122    def set_default_order(self, n):
     
    126126        msg = 'Default order must be either 1 or 2. I got %s' %n
    127127        assert n in [1,2], msg
    128        
     128
    129129        self.default_order = n
    130130        self.order = self.default_order
     
    412412        Input:
    413413          quantities: either None, a string or a list of strings naming the quantities to be reported
    414           tags:       either None, a string or a list of strings naming the tags to be reported         
     414          tags:       either None, a string or a list of strings naming the tags to be reported
    415415
    416416
     
    514514               yieldstep = None,
    515515               finaltime = None,
    516                duration = None,             
     516               duration = None,
    517517               skip_initial_step = False):
    518518        """Evolve model through time starting from self.starttime.
    519        
    520        
     519
     520
    521521        yieldstep: Interval between yields where results are stored,
    522522                   statistics written and domain inspected or
    523523                   possibly modified. If omitted the internal predefined
    524524                   max timestep is used.
    525                    Internally, smaller timesteps may be taken.               
     525                   Internally, smaller timesteps may be taken.
    526526
    527527        duration: Duration of simulation
     
    534534        skip_initial_step: Boolean flag that decides whether the first
    535535        yield step is skipped or not. This is useful for example to avoid
    536         duplicate steps when multiple evolve processes are dove tailed. 
     536        duplicate steps when multiple evolve processes are dove tailed.
    537537
    538538
     
    543543
    544544
    545         All times are given in seconds       
     545        All times are given in seconds
    546546
    547547        """
     
    575575            if duration is not None:
    576576                self.finaltime = self.starttime + float(duration)
    577                
    578 
    579        
     577
     578
     579
    580580
    581581        self.yieldtime = 0.0 #Time between 'yields'
     
    678678            if B is None:
    679679                print 'WARNING: Ignored boundary segment %d (None)'
    680             else:               
     680            else:
    681681                q = B.evaluate(vol_id, edge_id)
    682682
     
    693693    def update_timestep(self, yieldstep, finaltime):
    694694
    695         from config import min_timestep
     695        from config import min_timestep, max_timestep
    696696
    697697        # self.timestep is calculated from speed of characteristics
    698698        # Apply CFL condition here
    699         timestep = self.CFL*self.timestep
     699        timestep = min(self.CFL*self.timestep,max_timestep)
    700700
    701701        #Record maximal and minimal values of timestep for reporting
  • inundation/pyvolution/quantity.py

    r2526 r2648  
    197197          Arbitrary geo spatial dataset in the form of the class
    198198          Geospatial_data. Mesh points are populated using least squares
    199           fitting           
     199          fitting
    200200
    201201        points:
     
    210210
    211211        data_georef:
    212           If points is specified, geo_reference applies to each point.       
    213          
     212          If points is specified, geo_reference applies to each point.
     213
    214214        filename:
    215215          Name of a .pts file containing data points and attributes for
    216216          use with least squares.
    217          
     217
    218218        attribute_name:
    219219          If specified, any array matching that name
     
    289289        if not(points is None and values is None and data_georef is None):
    290290            from warnings import warn
    291            
     291
    292292            msg = 'Using points, values or data_georef with set_quantity '
    293293            msg += 'is obsolete. Please use a Geospatial_data object instead.'
     
    312312                                              location, indices, verbose)
    313313            elif isinstance(numeric, Geospatial_data):
    314                 self.set_values_from_geospatial_data(numeric,               
     314                self.set_values_from_geospatial_data(numeric,
    315315                                                     alpha,
    316316                                                     location, indices,
     
    338338            print 'The usage of points in set_values will be deprecated.' +\
    339339                  'Please use the geospatial_data object.'
    340            
     340
    341341            msg = 'When points are specified, associated values must also be.'
    342342            assert values is not None, msg
     
    350350                                      location, indices,
    351351                                      verbose = verbose,
    352                                       use_cache = use_cache) 
     352                                      use_cache = use_cache)
    353353        else:
    354354            raise 'This can\'t happen :-)'
     
    586586
    587587
    588        
     588
    589589    def set_values_from_geospatial_data(self, geospatial_data, alpha,
    590590                                        location, indices,
     
    604604                                    data_georef = data_georef,
    605605                                    verbose = verbose,
    606                                     use_cache = use_cache)                                   
    607 
    608                                
     606                                    use_cache = use_cache)
     607
     608
    609609
    610610    def set_values_from_points(self, points, values, alpha,
     
    621621        from least_squares import fit_to_mesh
    622622        from coordinate_transforms.geo_reference import Geo_reference
    623        
     623
    624624
    625625        points = ensure_numeric(points, Float)
     
    645645        #print data_georef
    646646        #print points
    647            
    648 
    649         #Call least squares method           
     647
     648
     649        #Call least squares method
    650650        args = (coordinates, triangles, points, values)
    651651        kwargs = {'data_origin': data_georef.get_origin(),
     
    656656
    657657        #print kwargs
    658        
     658
    659659        if use_cache is True:
    660660            try:
     
    672672            vertex_attributes = apply(fit_to_mesh,
    673673                                      args, kwargs)
    674            
    675         #Call underlying method using array values   
     674
     675        #Call underlying method using array values
    676676        self.set_values_from_array(vertex_attributes,
    677677                                   location, indices, verbose)
     
    730730            data_georef = None
    731731
    732        
     732
    733733
    734734        #Call underlying method for geospatial data
    735735        geospatial_data = points_dictionary2geospatial_data(points_dict)
    736736        geospatial_data.set_default_attribute_name(attribute_name)
    737            
     737
    738738        self.set_values_from_geospatial_data(geospatial_data,
    739739                                             alpha,
     
    741741                                             verbose = verbose,
    742742                                             use_cache = use_cache)
    743        
     743
    744744        #Call underlying method for points
    745745        #self.set_values_from_points(points, z, alpha,
     
    770770
    771771        The values will be stored in elements following their
    772         internal ordering.       
     772        internal ordering.
    773773
    774774        """
     
    11321132    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
    11331133
     1134    if sum(less(denominator, 1.0)) > 0.0:
     1135        msg = 'denominator < 1.0 in semi implicit update. Call Stephen :-)'
     1136        raise msg
     1137
    11341138    if sum(equal(denominator, 0.0)) > 0.0:
    11351139        msg = 'Zero division in semi implicit update. Call Stephen :-)'
  • inundation/pyvolution/quantity_ext.c

    r1506 r2648  
    8181
    8282      //Two point gradient
    83       _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);     
    84      
    85      
     83      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);
     84
     85
    8686      //Old (wrong code)
    8787      //det = x0*y1 - x1*y0;
     
    187187
    188188
    189   //Explicit updates
    190   for (k=0; k<N; k++) {
    191     centroid_values[k] += timestep*explicit_update[k];
    192   }
    193 
    194189  //Semi implicit updates
    195190  for (k=0; k<N; k++) {
     
    202197    }
    203198  }
     199
     200
     201  //Explicit updates
     202  for (k=0; k<N; k++) {
     203    centroid_values[k] += timestep*explicit_update[k];
     204  }
     205
     206
    204207  //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
    205   for (k=0;k<N;k++)
     208  for (k=0;k<N;k++){
    206209    semi_implicit_update[k]=0.0;
     210  }
     211
    207212  return 0;
    208213}
  • inundation/pyvolution/shallow_water.py

    r2620 r2648  
    8383        from config import minimum_allowed_height, maximum_allowed_speed, g
    8484        self.minimum_allowed_height = minimum_allowed_height
    85         self.maximum_allowed_speed = maximum_allowed_speed     
     85        self.maximum_allowed_speed = maximum_allowed_speed
    8686        self.g = g
    8787
     88
     89        self.forcing_terms.append(manning_friction)
    8890        self.forcing_terms.append(gravity)
    89         self.forcing_terms.append(manning_friction)
    9091
    9192        #Realtime visualisation
     
    108109        self.quantities_to_be_stored = ['stage','xmomentum','ymomentum']
    109110
    110        
     111
    111112    def set_quantities_to_be_stored(self, q):
    112113        """Specify which quantities will be stored in the sww file.
     
    118119
    119120        In the two first cases, the named quantities will be stored at each yieldstep
    120         (This is in addition to the quantities elevation and friction) 
     121        (This is in addition to the quantities elevation and friction)
    121122        If q is None, storage will be switched off altogether.
    122123        """
     
    124125
    125126        if q is None:
    126             self.quantities_to_be_stored = []           
     127            self.quantities_to_be_stored = []
    127128            self.store = False
    128129            return
     
    131132            q = [q] # Turn argument into a list
    132133
    133         #Check correcness   
     134        #Check correcness
    134135        for quantity_name in q:
    135136            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
    136             assert quantity_name in self.conserved_quantities, msg 
    137        
     137            assert quantity_name in self.conserved_quantities, msg
     138
    138139        self.quantities_to_be_stored = q
    139        
     140
    140141
    141142    def initialise_visualiser(self,scale_z=1.0,rect=None):
     
    224225               yieldstep = None,
    225226               finaltime = None,
    226                duration = None,               
     227               duration = None,
    227228               skip_initial_step = False):
    228229        """Specialisation of basic evolve method from parent class
     
    670671    hc = wc - zc  #Water depths at centroids
    671672
    672     #Update 
     673    #Update
    673674    #FIXME: Modify accroditg to c-version - or discard altogether.
    674675    for k in range(domain.number_of_elements):
     
    695696    from shallow_water_ext import protect
    696697
    697     protect(domain.minimum_allowed_height, domain.maximum_allowed_speed, 
     698    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
    698699            domain.epsilon, wc, zc, xmomc, ymomc)
    699700
     
    11831184
    11841185
    1185 def manning_friction_c(domain):
     1186def manning_friction_implicit_c(domain):
    11861187    """Wrapper for c version
    11871188    """
     
    12001201    xmom_update = xmom.semi_implicit_update
    12011202    ymom_update = ymom.semi_implicit_update
     1203
     1204    N = domain.number_of_elements
     1205    eps = domain.minimum_allowed_height
     1206    g = domain.g
     1207
     1208    from shallow_water_ext import manning_friction
     1209    manning_friction(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
     1210
     1211
     1212def manning_friction_explicit_c(domain):
     1213    """Wrapper for c version
     1214    """
     1215
     1216    xmom = domain.quantities['xmomentum']
     1217    ymom = domain.quantities['ymomentum']
     1218
     1219    w = domain.quantities['stage'].centroid_values
     1220    z = domain.quantities['elevation'].centroid_values
     1221
     1222    uh = xmom.centroid_values
     1223    vh = ymom.centroid_values
     1224    eta = domain.quantities['friction'].centroid_values
     1225
     1226    xmom_update = xmom.explicit_update
     1227    ymom_update = ymom.explicit_update
    12021228
    12031229    N = domain.number_of_elements
     
    17911817    extrapolate_second_order_sw=extrapolate_second_order_sw_c
    17921818    gravity = gravity_c
    1793     manning_friction = manning_friction_c
     1819    manning_friction = manning_friction_implicit_c
    17941820    h_limiter = h_limiter_c
    17951821    balance_deep_and_shallow = balance_deep_and_shallow_c
  • inundation/pyvolution/shallow_water_ext.c

    r2636 r2648  
    260260
    261261
     262void _manning_friction_explicit(double g, double eps, int N,
     263                       double* w, double* z,
     264                       double* uh, double* vh,
     265                       double* eta, double* xmom, double* ymom) {
     266
     267  int k;
     268  double S, h;
     269
     270  for (k=0; k<N; k++) {
     271    if (eta[k] > eps) {
     272      h = w[k]-z[k];
     273      if (h >= eps) {
     274        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
     275        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
     276        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     277        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
     278
     279
     280        //Update momentum
     281        xmom[k] += S*uh[k];
     282        ymom[k] += S*vh[k];
     283      }
     284    }
     285  }
     286}
     287
    262288
    263289int _balance_deep_and_shallow(int N,
     
    356382int _protect(int N,
    357383             double minimum_allowed_height,
    358              double maximum_allowed_speed,           
     384             double maximum_allowed_speed,
    359385             double epsilon,
    360386             double* wc,
     
    365391  int k;
    366392  double hc;
    367   double u, v, reduced_speed; 
     393  double u, v, reduced_speed;
    368394
    369395  //Protect against initesimal and negative heights
     
    372398
    373399    if (hc < minimum_allowed_height) {
    374    
     400
    375401      //Old code: Set momentum to zero and ensure h is non negative
    376402      //xmomc[k] = 0.0;
    377       //ymomc[k] = 0.0;     
    378       //if (hc <= 0.0) wc[k] = zc[k]; 
    379      
     403      //ymomc[k] = 0.0;
     404      //if (hc <= 0.0) wc[k] = zc[k];
     405
    380406
    381407      //New code: Adjust momentum to guarantee speeds are physical
    382       //          ensure h is non negative     
    383       //FIXME (Ole): This is only implemented in this C extension and 
    384       //             has no Python equivalent     
     408      //          ensure h is non negative
     409      //FIXME (Ole): This is only implemented in this C extension and
     410      //             has no Python equivalent
    385411      if (hc <= 0.0) {
    386         wc[k] = zc[k]; 
     412        wc[k] = zc[k];
    387413        xmomc[k] = 0.0;
    388414        ymomc[k] = 0.0;
    389415      } else {
    390         //Reduce excessive speeds derived from division by small hc 
     416        //Reduce excessive speeds derived from division by small hc
    391417        u = xmomc[k]/hc;
    392418        if (fabs(u) > maximum_allowed_speed) {
    393419          reduced_speed = maximum_allowed_speed * u/fabs(u);
    394           //printf("Speed (u) has been reduced from %.3f to %.3f\n", 
    395           //     u, reduced_speed);     
     420          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     421          //     u, reduced_speed);
    396422          xmomc[k] = reduced_speed * hc;
    397423        }
    398        
    399         v = ymomc[k]/hc;       
     424
     425        v = ymomc[k]/hc;
    400426        if (fabs(v) > maximum_allowed_speed) {
    401           reduced_speed = maximum_allowed_speed * v/fabs(v);   
    402           //printf("Speed (v) has been reduced from %.3f to %.3f\n",   
    403           //     v, reduced_speed);     
    404           ymomc[k] = reduced_speed * hc;         
    405         }       
     427          reduced_speed = maximum_allowed_speed * v/fabs(v);
     428          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     429          //     v, reduced_speed);
     430          ymomc[k] = reduced_speed * hc;
     431        }
    406432      }
    407433    }
     
    530556}
    531557
     558
     559PyObject *manning_friction_explicit(PyObject *self, PyObject *args) {
     560  //
     561  // manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
     562  //
     563
     564
     565  PyArrayObject *w, *z, *uh, *vh, *eta, *xmom, *ymom;
     566  int N;
     567  double g, eps;
     568
     569  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
     570                        &g, &eps, &w, &z, &uh, &vh, &eta,
     571                        &xmom, &ymom))
     572    return NULL;
     573
     574  N = w -> dimensions[0];
     575  _manning_friction_explicit(g, eps, N,
     576                    (double*) w -> data,
     577                    (double*) z -> data,
     578                    (double*) uh -> data,
     579                    (double*) vh -> data,
     580                    (double*) eta -> data,
     581                    (double*) xmom -> data,
     582                    (double*) ymom -> data);
     583
     584  return Py_BuildValue("");
     585}
    532586PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
    533587  /*Compute the vertex values based on a linear reconstruction on each triangle
     
    10671121  if (!PyArg_ParseTuple(args, "dddOOOO",
    10681122                        &minimum_allowed_height,
    1069                         &maximum_allowed_speed,                 
     1123                        &maximum_allowed_speed,
    10701124                        &epsilon,
    10711125                        &wc, &zc, &xmomc, &ymomc))
  • inundation/pyvolution/test_data_manager.py

    r2639 r2648  
    7979        latitudes = [-34.5, -34.33333, -34.16667, -34]
    8080
    81         long_name = 'LON'
     81        long_name = 'LON'
    8282        lat_name = 'LAT'
    8383
     
    8989        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
    9090            fid = NetCDFFile(self.test_MOST_file + ext, 'w')
    91        
     91
    9292            fid.createDimension(long_name,nx)
    9393            fid.createVariable(long_name,'d',(long_name,))
     
    107107            fid.variables['TIME'].units='seconds'
    108108            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
    109            
     109
    110110
    111111            name = ext[1:3].upper()
     
    113113            fid.createVariable(name,'d',('TIME', lat_name, long_name))
    114114            fid.variables[name].units='CENTIMETERS'
    115             fid.variables[name].missing_value=-1.e+034       
    116            
     115            fid.variables[name].missing_value=-1.e+034
     116
    117117            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
    118118                                              [-0.1214216, 0, 0, 0],
     
    139139                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
    140140                                              [0, 0.0004811212, 0.0004811212, 0]]])
    141            
    142        
     141
     142
    143143            fid.close()
    144144
    145        
     145
    146146
    147147
     
    150150        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
    151151            #print 'Trying to remove', self.test_MOST_file + ext
    152             os.remove(self.test_MOST_file + ext)       
     152            os.remove(self.test_MOST_file + ext)
    153153
    154154    def test_sww_constant(self):
     
    520520    #def test_write_pts(self):
    521521    #    #Obsolete
    522     #   
     522    #
    523523    #    #Get (enough) datapoints
    524524    #
     
    546546    #    #Check contents
    547547    #    #Get NetCDF
    548     #    from Scientific.IO.NetCDF import NetCDFFile       
     548    #    from Scientific.IO.NetCDF import NetCDFFile
    549549    #    fid = NetCDFFile(ptsfile, 'r')
    550550    #
     
    555555    #
    556556    #    #Check values#
    557     # 
     557    #
    558558    #    #print points[:]
    559559    #    #print ref_points
    560     #    assert allclose(points, points1) 
     560    #    assert allclose(points, points1)
    561561    #
    562562    #    #print attributes[:]
     
    569569    #    import os
    570570    #    os.remove(ptsfile)
    571        
     571
    572572
    573573    def test_dem2pts_bounding_box_v2(self):
     
    591591NODATA_value  -9999
    592592""")
    593         #Create linear function       
     593        #Create linear function
    594594        ref_points = []
    595595        ref_elevation = []
     
    610610
    611611        fid.close()
    612            
     612
    613613        #print 'sending pts', ref_points
    614614        #print 'sending elev', ref_elevation
     
    650650
    651651        #create new reference points
    652         newz = []       
     652        newz = []
    653653        newz[0:5] = ref_elevation[32:38]
    654654        newz[6:11] = ref_elevation[42:48]
     
    668668                x = x0 + xvec[j]
    669669                xnew = x - 2002.0
    670                 ref_points.append ([xnew,ynew]) #Relative point values       
     670                ref_points.append ([xnew,ynew]) #Relative point values
    671671
    672672        assert allclose(points, ref_points)
     
    704704NODATA_value  -9999
    705705""")
    706         #Create linear function       
     706        #Create linear function
    707707        ref_points = []
    708708        ref_elevation = []
     
    721721                ref_points.append ([x,y])
    722722                count += 1
    723                 z[count] = (4*i - 3*j)%13 
     723                z[count] = (4*i - 3*j)%13
    724724                if j == 4: z[count] = NODATA_value #column inside clipping region
    725725                if j == 8: z[count] = NODATA_value #column outside clipping region
    726726                if i == 9: z[count] = NODATA_value #row outside clipping region
    727                 if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region               
     727                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
    728728                ref_elevation.append( z[count] )
    729729                fid.write('%f ' %z[count])
    730730            fid.write('\n')
    731731
    732         fid.close()             
    733        
     732        fid.close()
     733
    734734        #print 'sending elev', ref_elevation
    735735
     
    781781        newz[16:19] = ref_elevation[65:68]
    782782
    783        
     783
    784784        ref_elevation = newz
    785785        ref_points = []
     
    795795                x = x0 + xvec[j]
    796796                xnew = x - 2002.0
    797                 if j <> 2 and (i<>1 or j<>4):   
     797                if j <> 2 and (i<>1 or j<>4):
    798798                    ref_points.append([x,y])
    799799                    new_ref_points.append ([xnew,ynew])
    800800
    801                
     801
    802802        assert allclose(points, new_ref_points)
    803803        assert allclose(elevation, ref_elevation)
     
    834834NODATA_value  -9999
    835835""")
    836         #Create linear function       
     836        #Create linear function
    837837        ref_points = []
    838838        ref_elevation = []
     
    851851                ref_points.append ([x,y])
    852852                count += 1
    853                 z[count] = (4*i - 3*j)%13 
     853                z[count] = (4*i - 3*j)%13
    854854                if j == 4: z[count] = NODATA_value #column inside clipping region
    855855                if j == 8: z[count] = NODATA_value #column outside clipping region
    856856                if i == 6: z[count] = NODATA_value #row on clipping boundary
    857                 if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region               
     857                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
    858858                ref_elevation.append( z[count] )
    859859                fid.write('%f ' %z[count])
    860860            fid.write('\n')
    861861
    862         fid.close()             
    863        
     862        fid.close()
     863
    864864        #print 'sending elev', ref_elevation
    865865
     
    910910
    911911
    912        
     912
    913913        ref_elevation = newz
    914914        ref_points = []
     
    924924                x = x0 + xvec[j]
    925925                xnew = x - 2002.0
    926                 if j <> 2 and (i<>1 or j<>4) and i<>3:   
     926                if j <> 2 and (i<>1 or j<>4) and i<>3:
    927927                    ref_points.append([x,y])
    928928                    new_ref_points.append ([xnew,ynew])
     
    932932        #print new_ref_points, len(new_ref_points)
    933933
    934         assert allclose(elevation, ref_elevation)                   
     934        assert allclose(elevation, ref_elevation)
    935935        assert allclose(points, new_ref_points)
    936936
     
    973973  PROJECTION ZONE: 50
    974974  DATUM: AGD66
    975   VERTICAL DATUM: 
    976   NUMBER OF REACHES:  19 
    977   NUMBER OF CROSS-SECTIONS:  2 
     975  VERTICAL DATUM:
     976  NUMBER OF REACHES:  19
     977  NUMBER OF CROSS-SECTIONS:  2
    978978END HEADER:
    979979
     
    984984    STREAM ID:Southern-Wungong
    985985    REACH ID:Southern-Wungong
    986     STATION:21410   
     986    STATION:21410
    987987    CUT LINE:
    988       407546.08 , 6437277.542 
    989       407329.32 , 6437489.482 
    990       407283.11 , 6437541.232 
     988      407546.08 , 6437277.542
     989      407329.32 , 6437489.482
     990      407283.11 , 6437541.232
    991991    SURFACE LINE:
    992992     407546.08,   6437277.54,   52.14
     
    999999
    10001000  CROSS-SECTION:
    1001     STREAM ID:Swan River     
    1002     REACH ID:Swan Mouth     
    1003     STATION:840.*   
     1001    STREAM ID:Swan River
     1002    REACH ID:Swan Mouth
     1003    STATION:840.*
    10041004    CUT LINE:
    1005       381178.0855 , 6452559.0685 
    1006       380485.4755 , 6453169.272 
     1005      381178.0855 , 6452559.0685
     1006      380485.4755 , 6453169.272
    10071007    SURFACE LINE:
    10081008     381178.09,   6452559.07,   4.17
     
    10171017     381063.46,   6452660.06,   3.67
    10181018     381054.41,   6452668.03,   3.67
    1019   END: 
     1019  END:
    10201020END CROSS-SECTIONS:
    10211021""")
     
    10441044                      [407510.08, 6437312.74]]
    10451045
    1046         ref_points += [[381178.09, 6452559.07], 
    1047                        [381169.49, 6452566.64], 
    1048                        [381157.78, 6452576.96], 
    1049                        [381155.97, 6452578.56], 
    1050                        [381143.72, 6452589.35], 
    1051                        [381136.69, 6452595.54], 
    1052                        [381114.74, 6452614.88], 
    1053                        [381075.53, 6452649.43], 
    1054                        [381071.47, 6452653.00], 
    1055                        [381063.46, 6452660.06], 
    1056                        [381054.41, 6452668.03]] 
    1057        
    1058                      
     1046        ref_points += [[381178.09, 6452559.07],
     1047                       [381169.49, 6452566.64],
     1048                       [381157.78, 6452576.96],
     1049                       [381155.97, 6452578.56],
     1050                       [381143.72, 6452589.35],
     1051                       [381136.69, 6452595.54],
     1052                       [381114.74, 6452614.88],
     1053                       [381075.53, 6452649.43],
     1054                       [381071.47, 6452653.00],
     1055                       [381063.46, 6452660.06],
     1056                       [381054.41, 6452668.03]]
     1057
     1058
    10591059        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
    10601060        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
     
    11271127                verbose = False,
    11281128                format = 'asc')
    1129                
     1129
    11301130        #Check prj (meta data)
    11311131        prjid = open(prjfile)
     
    12251225        cellsize      10.000000
    12261226        NODATA_value  -9999
    1227         -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200 
    1228          -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 
    1229          -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 
    1230          -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 
    1231          -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 
    1232          -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 
     1227        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
     1228         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
     1229         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
     1230         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
     1231         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
     1232         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
    12331233         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
    12341234         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
     
    12361236         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
    12371237           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
    1238        
     1238
    12391239        """
    12401240
     
    12671267        #
    12681268        domain.set_quantity('elevation', lambda x,y: -x-y)
    1269         domain.set_quantity('stage', 0)       
     1269        domain.set_quantity('stage', 0)
    12701270
    12711271        B = Transmissive_boundary(domain)
     
    12831283        cellsize = 10  #10m grid
    12841284
    1285        
     1285
    12861286        #Check contents
    12871287        #Get NetCDF
     
    13031303                verbose = False,
    13041304                format = 'asc')
    1305        
    1306                
     1305
     1306
    13071307        #Check prj (meta data)
    13081308        prjid = open(prjfile)
     
    13801380                #assert allclose(float(value), -(10-i+j)*cellsize)
    13811381                assert float(value) == -(10-i+j)*cellsize
    1382        
     1382
    13831383
    13841384        fid.close()
     
    13981398        Original extent is 100m x 100m:
    13991399
    1400         Eastings:   308500 -  308600 
     1400        Eastings:   308500 -  308600
    14011401        Northings: 6189000 - 6189100
    14021402
    14031403        Bounding box changes this to the 50m x 50m square defined by
    14041404
    1405         Eastings:   308530 -  308570 
     1405        Eastings:   308530 -  308570
    14061406        Northings: 6189050 - 6189100
    14071407
    14081408        The cropped values should be
    14091409
    1410          -130 -140 -150 -160 -170 
    1411          -120 -130 -140 -150 -160 
    1412          -110 -120 -130 -140 -150 
    1413          -100 -110 -120 -130 -140 
    1414           -90 -100 -110 -120 -130 
     1410         -130 -140 -150 -160 -170
     1411         -120 -130 -140 -150 -160
     1412         -110 -120 -130 -140 -150
     1413         -100 -110 -120 -130 -140
     1414          -90 -100 -110 -120 -130
    14151415          -80  -90 -100 -110 -120
    14161416
    1417         and the new lower reference point should be   
    1418         Eastings:   308530 
     1417        and the new lower reference point should be
     1418        Eastings:   308530
    14191419        Northings: 6189050
    1420        
     1420
    14211421        Original dataset is the same as in test_sww2dem_larger()
    1422        
     1422
    14231423        """
    14241424
     
    14511451        #
    14521452        domain.set_quantity('elevation', lambda x,y: -x-y)
    1453         domain.set_quantity('stage', 0)       
     1453        domain.set_quantity('stage', 0)
    14541454
    14551455        B = Transmissive_boundary(domain)
     
    14671467        cellsize = 10  #10m grid
    14681468
    1469        
     1469
    14701470        #Check contents
    14711471        #Get NetCDF
     
    14871487                easting_min = 308530,
    14881488                easting_max = 308570,
    1489                 northing_min = 6189050, 
    1490                 northing_max = 6189100,                               
     1489                northing_min = 6189050,
     1490                northing_max = 6189100,
    14911491                verbose = False,
    14921492                format = 'asc')
    1493        
    1494         fid.close()
    1495 
    1496        
     1493
     1494        fid.close()
     1495
     1496
    14971497        #Check prj (meta data)
    14981498        prjid = open(prjfile)
     
    15681568        for i, line in enumerate(lines[6:]):
    15691569            for j, value in enumerate( line.split() ):
    1570                 #assert float(value) == -(10-i+j)*cellsize               
     1570                #assert float(value) == -(10-i+j)*cellsize
    15711571                assert float(value) == -(10-i+j+3)*cellsize
    1572        
     1572
    15731573
    15741574
     
    20002000        #        cellsize = cellsize,
    20012001        #        verbose = False)
    2002                
     2002
    20032003        sww2dem(self.domain.filename,
    20042004                quantity = 'elevation',
     
    20102010        #Check header data
    20112011        from ermapper_grids import read_ermapper_header, read_ermapper_data
    2012        
     2012
    20132013        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
    20142014        #print header
    20152015        assert header['projection'].lower() == '"utm-56"'
    20162016        assert header['datum'].lower() == '"wgs84"'
    2017         assert header['units'].lower() == '"meters"'   
    2018         assert header['value'].lower() == '"elevation"'         
     2017        assert header['units'].lower() == '"meters"'
     2018        assert header['value'].lower() == '"elevation"'
    20192019        assert header['xdimension'] == '0.25'
    2020         assert header['ydimension'] == '0.25'   
     2020        assert header['ydimension'] == '0.25'
    20212021        assert float(header['eastings']) == 308500.0   #xllcorner
    2022         assert float(header['northings']) == 6189000.0 #yllcorner       
     2022        assert float(header['northings']) == 6189000.0 #yllcorner
    20232023        assert int(header['nroflines']) == 5
    2024         assert int(header['nrofcellsperline']) == 5     
     2024        assert int(header['nrofcellsperline']) == 5
    20252025        assert int(header['nullcellvalue']) == NODATA_value
    2026         #FIXME - there is more in the header                   
    2027 
    2028                
    2029         #Check grid data               
    2030         grid = read_ermapper_data(self.domain.filename + '_elevation') 
    2031        
     2026        #FIXME - there is more in the header
     2027
     2028
     2029        #Check grid data
     2030        grid = read_ermapper_data(self.domain.filename + '_elevation')
     2031
    20322032        #FIXME (Ole): Why is this the desired reference grid for -x-y?
    20332033        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
    20342034                    -1,    -1.25, -1.5,  -1.75, -2.0,
    2035                     -0.75, -1.0,  -1.25, -1.5,  -1.75,             
     2035                    -0.75, -1.0,  -1.25, -1.5,  -1.75,
    20362036                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
    2037                     -0.25, -0.5,  -0.75, -1.0,  -1.25]             
    2038                                          
     2037                    -0.25, -0.5,  -0.75, -1.0,  -1.25]
     2038
    20392039
    20402040        #print grid
     
    20422042
    20432043        fid.close()
    2044        
     2044
    20452045        #Cleanup
    20462046        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
    2047         #Done (Ole) - it was because sww2ers didn't close it's sww file 
     2047        #Done (Ole) - it was because sww2ers didn't close it's sww file
    20482048        os.remove(sww.filename)
    20492049        os.remove(self.domain.filename + '_elevation')
    2050         os.remove(self.domain.filename + '_elevation.ers')                 
     2050        os.remove(self.domain.filename + '_elevation.ers')
    20512051
    20522052
     
    20672067        # Fourth value (index==3) is -6.50198 cm
    20682068
    2069        
     2069
    20702070
    20712071        #Read
    20722072        from coordinate_transforms.redfearn import redfearn
    20732073        #fid = NetCDFFile(self.test_MOST_file)
    2074         fid = NetCDFFile(self.test_MOST_file + '_ha.nc')       
     2074        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
    20752075        first_value = fid.variables['HA'][:][0,0,0]
    20762076        fourth_value = fid.variables['HA'][:][0,0,3]
     
    20902090        #Read output file 'small.sww'
    20912091        #fid = NetCDFFile('small.sww')
    2092         fid = NetCDFFile(self.test_MOST_file + '.sww')       
     2092        fid = NetCDFFile(self.test_MOST_file + '.sww')
    20932093
    20942094        x = fid.variables['x'][:]
     
    22112211        h2_list = [-34.5,-34.33333]
    22122212
    2213         long_name = 'LON'
     2213        long_name = 'LON'
    22142214        lat_name = 'LAT'
    2215         time_name = 'TIME'       
     2215        time_name = 'TIME'
    22162216
    22172217        nx = 3
     
    22532253
    22542254        name = {}
    2255         name[fid1]='HA'
    2256         name[fid2]='UA'
    2257         name[fid3]='VA'
    2258         name[fid4]='ELEVATION'
     2255        name[fid1]='HA'
     2256        name[fid2]='UA'
     2257        name[fid3]='VA'
     2258        name[fid4]='ELEVATION'
    22592259
    22602260        units = {}
    2261         units[fid1]='cm'
    2262         units[fid2]='cm/s'
    2263         units[fid3]='cm/s'
    2264         units[fid4]='m'
    2265 
    2266         values = {}
    2267         values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
    2268         values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
    2269         values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
    2270         values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
     2261        units[fid1]='cm'
     2262        units[fid2]='cm/s'
     2263        units[fid3]='cm/s'
     2264        units[fid4]='m'
     2265
     2266        values = {}
     2267        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
     2268        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
     2269        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
     2270        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
    22712271
    22722272        for fid in [fid1,fid2,fid3]:
     
    22762276          fid.variables[name[fid]].assignValue(values[fid])
    22772277          fid.variables[name[fid]].missing_value = -99999999.
    2278           if fid == fid3: break
     2278          if fid == fid3: break
    22792279
    22802280        for fid in [fid4]:
    2281             fid.createVariable(name[fid],'d',(lat_name,long_name))
     2281            fid.createVariable(name[fid],'d',(lat_name,long_name))
    22822282            fid.variables[name[fid]].point_spacing='uneven'
    22832283            fid.variables[name[fid]].units=units[fid]
     
    23232323        elevation = fid.variables['elevation'][:]
    23242324        stage = fid.variables['stage'][:]
    2325         xmomentum = fid.variables['xmomentum'][:]
    2326         ymomentum = fid.variables['ymomentum'][:]
    2327 
    2328         #print ymomentum
     2325        xmomentum = fid.variables['xmomentum'][:]
     2326        ymomentum = fid.variables['ymomentum'][:]
     2327
     2328        #print ymomentum
    23292329        first_height = first_amp/100 - first_elevation
    23302330        third_height = third_amp/100 - third_elevation
     
    23712371        h2_list = [-34.5,-34.33333]
    23722372
    2373         long_name = 'LON961_1261'
    2374         lat_name = 'LAT481_841'
    2375         time_name = 'TIME1'
     2373#        long_name = 'LON961_1261'
     2374#        lat_name = 'LAT481_841'
     2375#        time_name = 'TIME1'
     2376
     2377        long_name = 'LON'
     2378        lat_name = 'LAT'
     2379        time_name = 'TIME'
    23762380
    23772381        nx = 3
     
    24132417
    24142418        name = {}
    2415         name[fid1]='HA'
    2416         name[fid2]='UA'
    2417         name[fid3]='VA'
    2418         name[fid4]='ELEVATION'
     2419        name[fid1]='HA'
     2420        name[fid2]='UA'
     2421        name[fid3]='VA'
     2422        name[fid4]='ELEVATION'
    24192423
    24202424        units = {}
    2421         units[fid1]='cm'
    2422         units[fid2]='cm/s'
    2423         units[fid3]='cm/s'
    2424         units[fid4]='m'
    2425 
    2426         values = {}
    2427         values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
    2428         values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
    2429         values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
    2430         values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
     2425        units[fid1]='cm'
     2426        units[fid2]='cm/s'
     2427        units[fid3]='cm/s'
     2428        units[fid4]='m'
     2429
     2430        values = {}
     2431        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
     2432        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
     2433        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
     2434        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
    24312435
    24322436        for fid in [fid1,fid2,fid3]:
     
    24362440          fid.variables[name[fid]].assignValue(values[fid])
    24372441          fid.variables[name[fid]].missing_value = -99999999.
    2438           if fid == fid3: break
     2442          if fid == fid3: break
    24392443
    24402444        for fid in [fid4]:
     
    24682472
    24692473        #Call conversion (with zero origin)
    2470         ferret2sww('test', verbose=False,
    2471                    origin = (56, 0, 0))
     2474        ferret2sww('test', verbose=False, origin = (56, 0, 0))
    24722475
    24732476        os.remove('test_va.nc')
     
    24832486        elevation = fid.variables['elevation'][:]
    24842487        stage = fid.variables['stage'][:]
    2485         xmomentum = fid.variables['xmomentum'][:]
    2486         ymomentum = fid.variables['ymomentum'][:]
    2487 
    2488         #print ymomentum
     2488        xmomentum = fid.variables['xmomentum'][:]
     2489        ymomentum = fid.variables['ymomentum'][:]
     2490
     2491        #print ymomentum
    24892492        first_height = first_amp/100 - first_elevation
    24902493        third_height = third_amp/100 - third_elevation
     
    25182521        #fid = NetCDFFile('small.sww', 'r')
    25192522        fid = NetCDFFile(self.test_MOST_file + '.sww')
    2520        
     2523
    25212524        x = fid.variables['x'][:]
    25222525        y = fid.variables['y'][:]
     
    26422645
    26432646        #os.remove(domain.filename + '.sww')
    2644         os.remove(filename)       
     2647        os.remove(filename)
    26452648
    26462649        bits = ['vertex_coordinates']
     
    27822785        #Clean up
    27832786        os.remove(filename)
    2784            
     2787
    27852788
    27862789        bits = [ 'geo_reference.get_xllcorner()',
     
    32043207
    32053208        os.remove(filename)
    3206        
     3209
    32073210    def test_asc_csiro2sww(self):
    32083211        import tempfile
    32093212
    32103213        bath_dir = tempfile.mkdtemp()
    3211         bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
     3214        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
    32123215        #bath_dir = 'bath_data_manager_test'
    3213         #print "os.getcwd( )",os.getcwd( ) 
     3216        #print "os.getcwd( )",os.getcwd( )
    32143217        elevation_dir =  tempfile.mkdtemp()
    32153218        #elevation_dir = 'elev_expanded'
    3216         elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
    3217         elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
    3218        
     3219        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
     3220        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
     3221
    32193222        fid = open(bath_dir_filename, 'w')
    32203223        fid.write(""" ncols             3
     
    32253228 nodata_value   -9999.0
    32263229    9000.000    -1000.000    3000.0
    3227    -1000.000    9000.000  -1000.000 
     3230   -1000.000    9000.000  -1000.000
    32283231""")
    32293232        fid.close()
    3230        
     3233
    32313234        fid = open(elevation_dir_filename1, 'w')
    32323235        fid.write(""" ncols             3
     
    32373240 nodata_value   -9999.0
    32383241    9000.000    0.000    3000.0
    3239      0.000     9000.000     0.000 
     3242     0.000     9000.000     0.000
    32403243""")
    32413244        fid.close()
     
    32493252 nodata_value   -9999.0
    32503253    9000.000    4000.000    4000.0
    3251     4000.000    9000.000    4000.000 
     3254    4000.000    9000.000    4000.000
    32523255""")
    32533256        fid.close()
    3254          
     3257
    32553258        ucur_dir =  tempfile.mkdtemp()
    32563259        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
    32573260        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
    3258        
     3261
    32593262        fid = open(ucur_dir_filename1, 'w')
    32603263        fid.write(""" ncols             3
     
    32653268 nodata_value   -9999.0
    32663269    90.000    60.000    30.0
    3267     10.000    10.000    10.000 
     3270    10.000    10.000    10.000
    32683271""")
    32693272        fid.close()
     
    32763279 nodata_value   -9999.0
    32773280    90.000    60.000    30.0
    3278     10.000    10.000    10.000 
     3281    10.000    10.000    10.000
    32793282""")
    32803283        fid.close()
    3281        
     3284
    32823285        vcur_dir =  tempfile.mkdtemp()
    32833286        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
    32843287        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
    3285        
     3288
    32863289        fid = open(vcur_dir_filename1, 'w')
    32873290        fid.write(""" ncols             3
     
    32923295 nodata_value   -9999.0
    32933296    90.000    60.000    30.0
    3294     10.000    10.000    10.000 
     3297    10.000    10.000    10.000
    32953298""")
    32963299        fid.close()
     
    33033306 nodata_value   -9999.0
    33043307    90.000    60.000    30.0
    3305     10.000    10.000    10.000 
     3308    10.000    10.000    10.000
    33063309""")
    33073310        fid.close()
    3308        
     3311
    33093312        sww_file = 'a_test.sww'
    33103313        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
    33113314
    33123315        # check the sww file
    3313        
     3316
    33143317        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
    33153318        x = fid.variables['x'][:]
     
    33193322        xmomentum = fid.variables['xmomentum'][:]
    33203323        geo_ref = Geo_reference(NetCDFObject=fid)
    3321         #print "geo_ref",geo_ref 
     3324        #print "geo_ref",geo_ref
    33223325        x_ref = geo_ref.get_xllcorner()
    33233326        y_ref = geo_ref.get_yllcorner()
     
    33253328        assert allclose(x_ref, 587798.418) # (-38, 148)
    33263329        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
    3327        
    3328         #Zone:   55   
    3329         #Easting:  588095.674  Northing: 5821451.722 
     3330
     3331        #Zone:   55
     3332        #Easting:  588095.674  Northing: 5821451.722
    33303333        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
    33313334        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
    33323335
    3333         #Zone:   55   
    3334         #Easting:  632145.632  Northing: 5820863.269 
    3335         #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
     3336        #Zone:   55
     3337        #Easting:  632145.632  Northing: 5820863.269
     3338        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
    33363339        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
    33373340
    3338         #Zone:   55   
    3339         #Easting:  609748.788  Northing: 5793447.860 
    3340         #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 '' 
     3341        #Zone:   55
     3342        #Easting:  609748.788  Northing: 5793447.860
     3343        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
    33413344        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
    33423345
     
    33463349        #(4000+1000)*60
    33473350        assert allclose(xmomentum[1][1],300000.0 )
    3348        
    3349        
    3350         fid.close()
    3351    
     3351
     3352
     3353        fid.close()
     3354
    33523355        #tidy up
    33533356        os.remove(bath_dir_filename)
     
    33573360        os.remove(elevation_dir_filename2)
    33583361        os.rmdir(elevation_dir)
    3359        
     3362
    33603363        os.remove(ucur_dir_filename1)
    33613364        os.remove(ucur_dir_filename2)
    33623365        os.rmdir(ucur_dir)
    3363        
     3366
    33643367        os.remove(vcur_dir_filename1)
    33653368        os.remove(vcur_dir_filename2)
     
    33693372        # remove sww file
    33703373        os.remove(sww_file)
    3371        
     3374
    33723375    def test_asc_csiro2sww2(self):
    33733376        import tempfile
    33743377
    33753378        bath_dir = tempfile.mkdtemp()
    3376         bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
     3379        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
    33773380        #bath_dir = 'bath_data_manager_test'
    3378         #print "os.getcwd( )",os.getcwd( ) 
     3381        #print "os.getcwd( )",os.getcwd( )
    33793382        elevation_dir =  tempfile.mkdtemp()
    33803383        #elevation_dir = 'elev_expanded'
    3381         elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
    3382         elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
    3383        
     3384        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
     3385        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
     3386
    33843387        fid = open(bath_dir_filename, 'w')
    33853388        fid.write(""" ncols             3
     
    33903393 nodata_value   -9999.0
    33913394    9000.000    -1000.000    3000.0
    3392    -1000.000    9000.000  -1000.000 
     3395   -1000.000    9000.000  -1000.000
    33933396""")
    33943397        fid.close()
    3395        
     3398
    33963399        fid = open(elevation_dir_filename1, 'w')
    33973400        fid.write(""" ncols             3
     
    34023405 nodata_value   -9999.0
    34033406    9000.000    0.000    3000.0
    3404      0.000     -9999.000     -9999.000 
     3407     0.000     -9999.000     -9999.000
    34053408""")
    34063409        fid.close()
     
    34143417 nodata_value   -9999.0
    34153418    9000.000    4000.000    4000.0
    3416     4000.000    9000.000    4000.000 
     3419    4000.000    9000.000    4000.000
    34173420""")
    34183421        fid.close()
    3419          
     3422
    34203423        ucur_dir =  tempfile.mkdtemp()
    34213424        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
    34223425        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
    3423        
     3426
    34243427        fid = open(ucur_dir_filename1, 'w')
    34253428        fid.write(""" ncols             3
     
    34303433 nodata_value   -9999.0
    34313434    90.000    60.000    30.0
    3432     10.000    10.000    10.000 
     3435    10.000    10.000    10.000
    34333436""")
    34343437        fid.close()
     
    34413444 nodata_value   -9999.0
    34423445    90.000    60.000    30.0
    3443     10.000    10.000    10.000 
     3446    10.000    10.000    10.000
    34443447""")
    34453448        fid.close()
    3446        
     3449
    34473450        vcur_dir =  tempfile.mkdtemp()
    34483451        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
    34493452        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
    3450        
     3453
    34513454        fid = open(vcur_dir_filename1, 'w')
    34523455        fid.write(""" ncols             3
     
    34573460 nodata_value   -9999.0
    34583461    90.000    60.000    30.0
    3459     10.000    10.000    10.000 
     3462    10.000    10.000    10.000
    34603463""")
    34613464        fid.close()
     
    34683471 nodata_value   -9999.0
    34693472    90.000    60.000    30.0
    3470     10.000    10.000    10.000 
     3473    10.000    10.000    10.000
    34713474""")
    34723475        fid.close()
    3473        
     3476
    34743477        try:
    34753478            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
     
    34833486            os.remove(elevation_dir_filename2)
    34843487            os.rmdir(elevation_dir)
    3485        
     3488
    34863489            os.remove(ucur_dir_filename1)
    34873490            os.remove(ucur_dir_filename2)
    34883491            os.rmdir(ucur_dir)
    3489        
     3492
    34903493            os.remove(vcur_dir_filename1)
    34913494            os.remove(vcur_dir_filename2)
     
    35003503            os.rmdir(elevation_dir)
    35013504            raise 'Should raise exception'
    3502        
     3505
    35033506            os.remove(ucur_dir_filename1)
    35043507            os.remove(ucur_dir_filename2)
    35053508            os.rmdir(ucur_dir)
    3506        
     3509
    35073510            os.remove(vcur_dir_filename1)
    35083511            os.remove(vcur_dir_filename2)
    35093512            os.rmdir(vcur_dir)
    35103513
    3511        
    3512      
     3514
     3515
    35133516    def test_asc_csiro2sww3(self):
    35143517        import tempfile
    35153518
    35163519        bath_dir = tempfile.mkdtemp()
    3517         bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
     3520        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
    35183521        #bath_dir = 'bath_data_manager_test'
    3519         #print "os.getcwd( )",os.getcwd( ) 
     3522        #print "os.getcwd( )",os.getcwd( )
    35203523        elevation_dir =  tempfile.mkdtemp()
    35213524        #elevation_dir = 'elev_expanded'
    3522         elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
    3523         elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
    3524        
     3525        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
     3526        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
     3527
    35253528        fid = open(bath_dir_filename, 'w')
    35263529        fid.write(""" ncols             3
     
    35313534 nodata_value   -9999.0
    35323535    9000.000    -1000.000    3000.0
    3533    -1000.000    9000.000  -1000.000 
     3536   -1000.000    9000.000  -1000.000
    35343537""")
    35353538        fid.close()
    3536        
     3539
    35373540        fid = open(elevation_dir_filename1, 'w')
    35383541        fid.write(""" ncols             3
     
    35433546 nodata_value   -9999.0
    35443547    9000.000    0.000    3000.0
    3545      0.000     -9999.000     -9999.000 
     3548     0.000     -9999.000     -9999.000
    35463549""")
    35473550        fid.close()
     
    35553558 nodata_value   -9999.0
    35563559    9000.000    4000.000    4000.0
    3557     4000.000    9000.000    4000.000 
     3560    4000.000    9000.000    4000.000
    35583561""")
    35593562        fid.close()
    3560          
     3563
    35613564        ucur_dir =  tempfile.mkdtemp()
    35623565        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
    35633566        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
    3564        
     3567
    35653568        fid = open(ucur_dir_filename1, 'w')
    35663569        fid.write(""" ncols             3
     
    35713574 nodata_value   -9999.0
    35723575    90.000    60.000    30.0
    3573     10.000    10.000    10.000 
     3576    10.000    10.000    10.000
    35743577""")
    35753578        fid.close()
     
    35823585 nodata_value   -9999.0
    35833586    90.000    60.000    30.0
    3584     10.000    10.000    10.000 
     3587    10.000    10.000    10.000
    35853588""")
    35863589        fid.close()
    3587        
     3590
    35883591        vcur_dir =  tempfile.mkdtemp()
    35893592        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
    35903593        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
    3591        
     3594
    35923595        fid = open(vcur_dir_filename1, 'w')
    35933596        fid.write(""" ncols             3
     
    35983601 nodata_value   -9999.0
    35993602    90.000    60.000    30.0
    3600     10.000    10.000    10.000 
     3603    10.000    10.000    10.000
    36013604""")
    36023605        fid.close()
     
    36093612 nodata_value   -9999.0
    36103613    90.000    60.000    30.0
    3611     10.000    10.000    10.000 
     3614    10.000    10.000    10.000
    36123615""")
    36133616        fid.close()
    3614        
     3617
    36153618        sww_file = 'a_test.sww'
    36163619        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
     
    36193622
    36203623        # check the sww file
    3621        
     3624
    36223625        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
    36233626        x = fid.variables['x'][:]
     
    36273630        xmomentum = fid.variables['xmomentum'][:]
    36283631        geo_ref = Geo_reference(NetCDFObject=fid)
    3629         #print "geo_ref",geo_ref 
     3632        #print "geo_ref",geo_ref
    36303633        x_ref = geo_ref.get_xllcorner()
    36313634        y_ref = geo_ref.get_yllcorner()
     
    36333636        assert allclose(x_ref, 587798.418) # (-38, 148)
    36343637        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
    3635        
    3636         #Zone:   55   
    3637         #Easting:  588095.674  Northing: 5821451.722 
     3638
     3639        #Zone:   55
     3640        #Easting:  588095.674  Northing: 5821451.722
    36383641        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
    36393642        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
    36403643
    3641         #Zone:   55   
    3642         #Easting:  632145.632  Northing: 5820863.269 
    3643         #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
     3644        #Zone:   55
     3645        #Easting:  632145.632  Northing: 5820863.269
     3646        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
    36443647        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
    36453648
    3646         #Zone:   55   
    3647         #Easting:  609748.788  Northing: 5793447.860 
    3648         #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 '' 
     3649        #Zone:   55
     3650        #Easting:  609748.788  Northing: 5793447.860
     3651        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
    36493652        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
    36503653
     
    36553658        #(100.0 - 9000)*10
    36563659        assert allclose(xmomentum[0][4], -89000.0 )
    3657        
     3660
    36583661        #(100.0 - -1000.000)*10
    36593662        assert allclose(xmomentum[0][5], 11000.0 )
    3660        
    3661         fid.close()
    3662    
     3663
     3664        fid.close()
     3665
    36633666        #tidy up
    36643667        os.remove(bath_dir_filename)
     
    36683671        os.remove(elevation_dir_filename2)
    36693672        os.rmdir(elevation_dir)
    3670        
     3673
    36713674        os.remove(ucur_dir_filename1)
    36723675        os.remove(ucur_dir_filename2)
    36733676        os.rmdir(ucur_dir)
    3674        
     3677
    36753678        os.remove(vcur_dir_filename1)
    36763679        os.remove(vcur_dir_filename2)
     
    36793682        # remove sww file
    36803683        os.remove(sww_file)
    3681    
    3682      
     3684
     3685
    36833686    def test_asc_csiro2sww4(self):
    36843687        """
    36853688        Test specifying the extent
    36863689        """
    3687        
     3690
    36883691        import tempfile
    36893692
    36903693        bath_dir = tempfile.mkdtemp()
    3691         bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
     3694        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
    36923695        #bath_dir = 'bath_data_manager_test'
    3693         #print "os.getcwd( )",os.getcwd( ) 
     3696        #print "os.getcwd( )",os.getcwd( )
    36943697        elevation_dir =  tempfile.mkdtemp()
    36953698        #elevation_dir = 'elev_expanded'
    3696         elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
    3697         elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
    3698        
     3699        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
     3700        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
     3701
    36993702        fid = open(bath_dir_filename, 'w')
    37003703        fid.write(""" ncols             4
     
    37103713""")
    37113714        fid.close()
    3712        
     3715
    37133716        fid = open(elevation_dir_filename1, 'w')
    37143717        fid.write(""" ncols             4
     
    37383741""")
    37393742        fid.close()
    3740          
     3743
    37413744        ucur_dir =  tempfile.mkdtemp()
    37423745        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
    37433746        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
    3744        
     3747
    37453748        fid = open(ucur_dir_filename1, 'w')
    37463749        fid.write(""" ncols             4
     
    37693772""")
    37703773        fid.close()
    3771        
     3774
    37723775        vcur_dir =  tempfile.mkdtemp()
    37733776        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
    37743777        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
    3775        
     3778
    37763779        fid = open(vcur_dir_filename1, 'w')
    37773780        fid.write(""" ncols             4
     
    38003803""")
    38013804        fid.close()
    3802        
     3805
    38033806        sww_file = tempfile.mktemp(".sww")
    38043807        #sww_file = 'a_test.sww'
     
    38123815
    38133816        # check the sww file
    3814        
     3817
    38153818        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
    38163819        x = fid.variables['x'][:]
     
    38213824        ymomentum = fid.variables['ymomentum'][:]
    38223825        geo_ref = Geo_reference(NetCDFObject=fid)
    3823         #print "geo_ref",geo_ref 
     3826        #print "geo_ref",geo_ref
    38243827        x_ref = geo_ref.get_xllcorner()
    38253828        y_ref = geo_ref.get_yllcorner()
     
    38303833        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
    38313834
    3832         #Easting:  632145.632  Northing: 5820863.269 
    3833         #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
    3834        
     3835        #Easting:  632145.632  Northing: 5820863.269
     3836        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
     3837
    38353838        #print "x",x
    3836         #print "y",y 
     3839        #print "y",y
    38373840        self.failUnless(len(x) == 4,'failed') # 2*2
    38383841        self.failUnless(len(x) == 4,'failed') # 2*2
    3839        
     3842
    38403843        #Zone:   55
    3841         #Easting:  632145.632  Northing: 5820863.269 
     3844        #Easting:  632145.632  Northing: 5820863.269
    38423845        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
    3843         # magic number - y is close enough for me. 
     3846        # magic number - y is close enough for me.
    38443847        assert allclose(x[3], 632145.63 - x_ref)
    38453848        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
    38463849
    38473850        assert allclose(z[0],9000.0 ) #z is elevation info
    3848         #print "z",z 
     3851        #print "z",z
    38493852        # 2 time steps, 4 points
    3850         self.failUnless(xmomentum.shape == (2,4), 'failed') 
    3851         self.failUnless(ymomentum.shape == (2,4), 'failed') 
    3852        
     3853        self.failUnless(xmomentum.shape == (2,4), 'failed')
     3854        self.failUnless(ymomentum.shape == (2,4), 'failed')
     3855
    38533856        #(100.0 - -1000.000)*10
    38543857        #assert allclose(xmomentum[0][5], 11000.0 )
    3855        
    3856         fid.close()
    3857        
     3858
     3859        fid.close()
     3860
    38583861        # is the sww file readable?
    38593862        #Lets see if we can convert it to a dem!
    3860         #print "sww_file",sww_file 
     3863        #print "sww_file",sww_file
    38613864        #dem_file = tempfile.mktemp(".dem")
    38623865        domain = sww2domain(sww_file) ###, dem_file)
     
    38703873        os.remove(elevation_dir_filename2)
    38713874        os.rmdir(elevation_dir)
    3872        
     3875
    38733876        os.remove(ucur_dir_filename1)
    38743877        os.remove(ucur_dir_filename2)
    38753878        os.rmdir(ucur_dir)
    3876        
     3879
    38773880        os.remove(vcur_dir_filename1)
    38783881        os.remove(vcur_dir_filename2)
     
    38803883
    38813884
    3882      
    3883        
     3885
     3886
    38843887        # remove sww file
    38853888        os.remove(sww_file)
    3886        
     3889
    38873890        # remove dem file
    38883891        #os.remove(dem_file)
     
    39073910                        longitudes == longitudes_news,
    39083911                         'failed')
    3909        
     3912
    39103913        ## 2nd test
    39113914        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
     
    39183921        #print "latitudes_new", latitudes_new
    39193922        #print "longitudes_news",longitudes_news
    3920        
     3923
    39213924        self.failUnless(latitudes == latitudes_new and \
    39223925                        longitudes == longitudes_news,
    39233926                         'failed')
    3924        
     3927
    39253928        ## 3rd test
    39263929        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
     
    39333936        #print "latitudes_new", latitudes_new
    39343937        #print "longitudes_news",longitudes_news
    3935        
     3938
    39363939        self.failUnless(latitudes_new == [2, 1] and \
    39373940                        longitudes_news == [10, 20],
    39383941                         'failed')
    39393942
    3940        
     3943
    39413944        ## 4th test
    39423945        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
     
    39493952        #print "latitudes_new", latitudes_new
    39503953        #print "longitudes_news",longitudes_news
    3951        
     3954
    39523955        self.failUnless(latitudes_new == [2, 1, 0] and \
    39533956                        longitudes_news == [0, 10, 20],
     
    39633966        #print "latitudes_new", latitudes_new
    39643967        #print "longitudes_news",longitudes_news
    3965        
     3968
    39663969        self.failUnless(latitudes_new == [2, 1, 0] and \
    39673970                        longitudes_news == [0, 10, 20],
    39683971                         'failed')
    3969        
     3972
    39703973        ## 6th test
    3971        
     3974
    39723975        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
    39733976            latitudes,longitudes,
     
    39793982        #print "latitudes_new", latitudes_new
    39803983        #print "longitudes_news",longitudes_news
    3981        
     3984
    39823985        self.failUnless(latitudes_new == [3, 2, 1] and \
    39833986                        longitudes_news == [10, 20, 30],
    39843987                         'failed')
    3985        
    3986        
     3988
     3989
    39873990        ## 7th test
    39883991        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
     
    39953998        longitudes_news = longitudes[lmin:lmax]
    39963999        m2d = m2d[kmin:kmax,lmin:lmax]
    3997         #print "m2d", m2d 
     4000        #print "m2d", m2d
    39984001        #print "latitudes_new", latitudes_new
    39994002        #print "longitudes_news",longitudes_news
    4000        
     4003
    40014004        self.failUnless(latitudes_new == [2, 1] and \
    40024005                        longitudes_news == [10, 20],
    40034006                         'failed')
    4004        
     4007
    40054008        self.failUnless(m2d == [[5,6],[9,10]],
    40064009                         'failed')
    4007        
     4010
    40084011    def test_get_min_max_indexes2(self):
    40094012        latitudes = [-30,-35,-40,-45]
     
    40114014
    40124015        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
    4013        
     4016
    40144017        # k - lat
    40154018        # l - lon
     
    40204023        #print "kmin",kmin;print "kmax",kmax
    40214024        #print "lmin",lmin;print "lmax",lmax
    4022         #print "m2d", m2d 
     4025        #print "m2d", m2d
    40234026        #print "latitudes", latitudes
    40244027        #print "longitudes",longitudes
    4025         #print "latitudes[kmax]", latitudes[kmax] 
     4028        #print "latitudes[kmax]", latitudes[kmax]
    40264029        latitudes_new = latitudes[kmin:kmax]
    40274030        longitudes_new = longitudes[lmin:lmax]
     
    40304033        #print "latitudes_new", latitudes_new
    40314034        #print "longitudes_new",longitudes_new
    4032        
     4035
    40334036        self.failUnless(latitudes_new == [-30, -35, -40] and \
    40344037                        longitudes_new == [148, 149,150],
     
    40364039        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
    40374040                         'failed')
    4038        
     4041
    40394042    def test_get_min_max_indexes3(self):
    40404043        latitudes = [-30,-35,-40,-45,-50,-55,-60]
     
    40474050            -43,-37,148.5,149.5)
    40484051
    4049        
     4052
    40504053        #print "kmin",kmin;print "kmax",kmax
    40514054        #print "lmin",lmin;print "lmax",lmax
     
    40564059        #print "latitudes_new", latitudes_new
    40574060        #print "longitudes_news",longitudes_news
    4058        
     4061
    40594062        self.failUnless(latitudes_new == [-35, -40, -45] and \
    40604063                        longitudes_news == [148, 149,150],
    40614064                         'failed')
    4062        
     4065
    40634066    def test_get_min_max_indexes4(self):
    40644067        latitudes = [-30,-35,-40,-45,-50,-55,-60]
     
    40694072        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
    40704073            latitudes,longitudes)
    4071        
     4074
    40724075        #print "kmin",kmin;print "kmax",kmax
    40734076        #print "lmin",lmin;print "lmax",lmax
     
    40784081        #print "latitudes_new", latitudes_new
    40794082        #print "longitudes_news",longitudes_news
    4080        
     4083
    40814084        self.failUnless(latitudes_new == latitudes  and \
    40824085                        longitudes_news == longitudes,
    40834086                         'failed')
    4084        
     4087
    40854088    def test_tsh2sww(self):
    40864089        import os
    40874090        import tempfile
    4088        
     4091
    40894092        tsh_file = tempfile.mktemp(".tsh")
    40904093        file = open(tsh_file,"w")
     
    41244127
    41254128        #sww_file = tempfile.mktemp(".sww")
    4126         #print "sww_file",sww_file 
     4129        #print "sww_file",sww_file
    41274130        #print "sww_file",tsh_file
    41284131        tsh2sww(tsh_file)
     
    41354138    suite = unittest.makeSuite(Test_Data_Manager,'test')
    41364139    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
    4137     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')   
     4140    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')
    41384141    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
    41394142    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
     
    41414144    runner = unittest.TextTestRunner()
    41424145    runner.run(suite)
    4143  
     4146
    41444147
    41454148
     
    43784381                if i == 3 and j == 3: z = NODATA_value
    43794382
    4380                
     4383
    43814384                if z <> NODATA_value:
    43824385                    new_ref_elev.append(z)
    43834386                    new_ref_pts.append( [x,y] )
    4384                
     4387
    43854388                ref_points.append( [x,y] )
    43864389                ref_elevation.append(z)
     
    43914394        fid.close()
    43924395
    4393        
     4396
    43944397        #Write prj file with metadata
    43954398        metafilename = root+'.prj'
     
    44824485                    new_ref_elev1.append(z)
    44834486                    new_ref_pts1.append( [x,y] )
    4484                    
     4487
    44854488                ref_points.append( [x,y] )
    44864489                ref_elevation.append(z)
     
    45064509""")
    45074510        fid.close()
    4508        
     4511
    45094512        #Convert to NetCDF pts
    45104513        convert_dem_from_ascii2netcdf(root)
     
    45494552                    new_ref_pts2.append( [x_new,y_new] )
    45504553
    4551                    
     4554
    45524555                ref_points.append( [x_new,y_new] )
    45534556                ref_elevation.append(z)
     
    45614564        #assert allclose(elevation, ref_elevation)
    45624565
    4563        
     4566
    45644567        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
    45654568        assert allclose(points, new_ref_pts2), 'points do not align'
     
    45804583
    45814584
    4582 
  • inundation/zeus/Validation.zpi

    r2255 r2648  
    7373    <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
    7474    <ReleaseProjectReload>Off</ReleaseProjectReload>
    75     <file>..\validation\completed\lwru2\create_mesh.py</file>
    76     <file>..\validation\completed\lwru2\extract_timeseries.py</file>
    77     <file>..\validation\completed\lwru2\lwru2.py</file>
    78     <file>..\validation\completed\lwru2\project.py</file>
    7975    <folder name="Header Files" />
    8076    <folder name="Resource Files" />
  • inundation/zeus/analytic_solutions.zpi

    r2255 r2648  
    7373    <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
    7474    <ReleaseProjectReload>Off</ReleaseProjectReload>
    75     <file>..\analytical solutions\Analytical solution_oblique_shock.py</file>
    76     <file>..\analytical solutions\Analytical solution_oblique_shock_01.py</file>
    77     <file>..\analytical solutions\Analytical_solution_circular_hydraulic_jump.py</file>
    78     <file>..\analytical solutions\Non_symmetrical_dam_break.py</file>
     75    <file>..\..\development\analytical solutions\Analytical_solution_circular_hydraulic_jump.py</file>
     76    <file>..\..\development\analytical solutions\Analytical_solution_Sampson_cross_mesh.py</file>
     77    <file>..\..\development\analytical solutions\Analytical_solution_Yoon_cross_mesh.py</file>
    7978</project>
  • inundation/zeus/anuga-workspace.ini

    r2255 r2648  
    11[Files MRU]
    2 0=c:\home\projects\anuga\inundation\pyvolution\shallow_water.py
    3 1=c:\home\projects\anuga\inundation\pyvolution\general_mesh.py
    4 2=c:\apps\lyx\bin\reLyX.bat
    5 3=c:\program files\zfw\readme.txt
     20=c:\home\projects\anuga\production\merimbula_2005\run_merimbula_lake.py
     31=c:\home\projects\anuga\production\merimbula_2005\run_new_meribula.py
     42=c:\home\projects\anuga\production\merimbula_2005\prepare.py
     53=c:\home\projects\anuga\inundation\pyvolution\shallow_water.py
     64=c:\home\projects\anuga\inundation\pyvolution\general_mesh.py
     75=c:\apps\lyx\bin\reLyX.bat
     86=c:\program files\zfw\readme.txt
  • inundation/zeus/anuga-workspace.zwi

    r2622 r2648  
    22<workspace name="anuga-workspace">
    33    <mode>Debug</mode>
    4     <active>Merimbula</active>
     4    <active>analytic_solutions</active>
    55    <project name="analytic_solutions">analytic_solutions.zpi</project>
    66    <project name="euler">euler.zpi</project>
  • inundation/zeus/parallel.zpi

    r2255 r2648  
    8484    <file>..\parallel\run_parallel_advection.py</file>
    8585    <file>..\parallel\run_parallel_merimbula.py</file>
    86     <file>..\parallel\run_parallel_mesh.py</file>
    8786    <file>..\parallel\run_parallel_sw_merimbula.py</file>
    8887    <file>..\parallel\run_parallel_sw_merimbula_metis.py</file>
  • production/merimbula_2005/project.py

    r2622 r2648  
    77bathymetry_filename = 'merimbula_bathymetry.xya'
    88mesh_filename = 'merimbula_10785.tsh'
     9sww_file = 'merimbula_weed.sww'
  • production/merimbula_2005/run_new_meribula_weed.py

    r2622 r2648  
    156156# close to bearing = 0 and 360!
    157157#---------------------------------------------
    158 filename = project.original_wind_filename[:-4]+'.tms'
    159 print 'Wind field from ',filename
    160 wind = file_function(filename, domain=domain, quantities=['speed','bearing'])
    161 domain.forcing_terms.append(Wind_stress(wind))
     158#filename = project.original_wind_filename[:-4]+'.tms'
     159#print 'Wind field from ',filename
     160#wind = file_function(filename, domain=domain, quantities=['speed','bearing'])
     161#domain.forcing_terms.append(Wind_stress(wind))
    162162
    163163
Note: See TracChangeset for help on using the changeset viewer.