Changeset 4280

Timestamp:

Feb 26, 2007, 9:50:44 AM (18 years ago)

Author:

duncan

Message:

working on updated urs2sww

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (18 diffs)
• test_data_manager.py (modified) (1 diff)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -78 +78 @@
 from anuga.caching.caching import myhash
 from anuga.utilities.anuga_exceptions import ANUGAError
-from anuga.pmesh.mesh import Mesh
 from anuga.shallow_water import Domain
 from anuga.abstract_2d_finite_volumes.pmesh2domain import \
@@ -4419 +4418 @@
 
     """
+    #FIXME cache this function!
     
     from sets import ImmutableSet
@@ -4449 +4449 @@
     interpolate any point on the segment.
     """
+    from math import sqrt
     #print "zone",zone 
     geo_reference = Geo_reference(zone=zone)
@@ -4456 +4457 @@
     #print "seg_lat_long", seg_lat_long
     # 1.415 = 2^0.5, rounded up....
-    buffer = 1.415 * grid_spacing
+    sqrt_2_rounded_up = 1.415
+    buffer = sqrt_2_rounded_up * grid_spacing
     
     #
@@ -4465 +4467 @@
     min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
 
-    #print "ll_lat", ll_lat
-    #print "ll_long", ll_long
-    #print "max_lat", max_lat
-    #print "max_long", max_long 
-    #print "min_lat", min_lat
-    #print "min_long", min_long
     first_row = (min_long - ll_long)/grid_spacing
     # To round up
@@ -4479 +4475 @@
     last_row_long = int(round(last_row))
     #print "last_row",last_row _long
-
     
     first_row = (min_lat - ll_lat)/grid_spacing
@@ -4490 +4485 @@
     #print "last_row",last_row_lat
 
+    ul_lat = ll_lat + grid_spacing*lat_amount
+    ul_lat_close = ul_lat - grid_spacing
+    ll_long_close = ll_long + grid_spacing    
+    _ , x, y = redfearn(ul_lat, ll_long)    
+    _ , _, y_close = redfearn(ul_lat_close, ll_long)    
+    _ , x_close, _ = redfearn(ul_lat, ll_long_close)
+
+    
+    max_distance = sqrt_2_rounded_up * max(x-x_close,y - y_close)
     points_lat_long = []
     # Create a list of the lat long points to include.
@@ -4496 +4500 @@
             lat = ll_lat + index_lat*grid_spacing
             long = ll_long + index_long*grid_spacing
-            points_lat_long.append((lat, long)) #must be hashable
+
+            #filter here to keep good points
+            if keep_point(lat, long, seg, max_distance):
+                points_lat_long.append((lat, long)) #must be hashable
     #print "points_lat_long", points_lat_long
-    return points_lat_long    
-
+
+    # Now that we have these points, lets throw ones out that are too far away
+    return points_lat_long
+
+def keep_point(lat, long, seg, max_distance):
+    """
+    seg is two points, UTM
+    """
+    from math import sqrt
+    _ , x0, y0 = redfearn(lat, long)
+    x1 = seg[0][0]
+    y1 = seg[0][1]
+    x2 = seg[1][0]
+    y2 = seg[1][1]
+
+    x2_1 = x2-x1
+    y2_1 = y2-y1
+    d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1))
+    if d <= max_distance:
+        return True
+    else:
+        return False
+    
     #### CONVERTING UNGRIDDED URS DATA TO AN SWW FILE ####
 def urs_ungridded2sww_link(basename_in='o', basename_out=None, verbose=False,
@@ -4522 +4550 @@
             origin = None,
             zscale=1,
-            fail_on_NaN=True,
-            NaN_filler=0):
+            fail_on_NaN=True):
     """
     parameters not used!
-    NaN_filler
     origin
     #mint=None, maxt=None,
-    """
-    # Convert to utm
-    
+    """    
+    from anuga.pmesh.mesh import Mesh
+
     files_in = [basename_in+'-z-mux',
                 basename_in+'-e-mux',
@@ -4539 +4565 @@
     # instanciate urs_points of the three mux files.
     mux = {}
-    quality_slices = {}
+    #quality_slices = {}
     for quantity, file in map(None, quantities, files_in):
         mux[quantity] = Urs_points(file)
-        quality_slices[quantity] = []
-        for slice in mux[quantity]:
-            quality_slices[quantity].append(slice)
-        
+        #quality_slices[quantity] = []
+        #for slice in mux[quantity]:
+        #    quality_slices[quantity].append(slice)
         
     # FIXME: check that the depth is the same. (hashing)
@@ -4560 +4585 @@
     #print "zone", zone
 
-    elevation = a_mux.lonlatdep[:,0] * -1 #
+    elevation = a_mux.lonlatdep[:,2] * -1 #
     
     # grid ( create a mesh from the selected points)
+    # This mesh has a problem.  Triangles are streched over ungridded areas.
+    #  If these areas could be described as holes in pmesh, that would be great
+    
     mesh = Mesh()
     mesh.add_vertices(points_utm)
@@ -4570 +4598 @@
 
     #mesh.export_mesh_file(basename_in + '.tsh')
-
+    
     times = []
     for i in range(a_mux.time_step_count):
@@ -4593 +4621 @@
     write_sww_header(outfile, times, len(volumes), len(points_utm))
     write_sww_triangulation(outfile, points_utm, volumes,
-                            elevation, zone,  verbose=verbose)
-
-    write_sww_time_slice(outfile, quality_slices['HA'],
-                         quality_slices['UA'], quality_slices['VA'],
-                         elevation,
-                         mean_stage=mean_stage, zscale=zscale,
-                         verbose=verbose)
+                            elevation, zone,  origin=origin, verbose=verbose)
+    j = 0
+    # Read in a time slice from each mux file and write it to the sww file
+    for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):            
+        write_sww_time_slice(outfile, HA, UA, VA,  elevation, j,
+                              mean_stage=mean_stage, zscale=zscale,
+                              verbose=verbose)
+        j += 1
+    outfile.close()
     #
     # Do some conversions while writing the sww file
-
-
 
 def write_sww_header(outfile, times, number_of_volumes, number_of_points ):
@@ -4659 +4687 @@
                            ('number_of_timesteps',
                             'number_of_points'))
-
-class Urs_points:
-    """
-    Read the info in URS mux files.
-
-    for the quantities heres a correlation between the file names and
-    what they mean;
-    z-mux is height above sea level, m
-    e-mux is velocity is Eastern direction, m/s
-    n-mux is velocity is Northern direction, m/s
-
-    
-    """
-    def __init__(self,urs_file):
-        self.iterated = False
-        columns = 3 # long, lat , depth
-        mux_file = open(urs_file, 'rb')
-        
-        # Number of points/stations
-        (self.points_num,)= unpack('i',mux_file.read(4))
-        
-        # nt, int - Number of time steps
-        (self.time_step_count,)= unpack('i',mux_file.read(4))
-
-        #dt, float - time step, seconds
-        (self.time_step,) = unpack('f', mux_file.read(4))
-    
-        msg = "Bad data in the urs file."
-        if self.points_num < 0:
-            mux_file.close()
-            raise ANUGAError, msg
-        if self.time_step_count < 0:
-            mux_file.close()
-            raise ANUGAError, msg
-        if self.time_step < 0:
-            mux_file.close()
-            raise ANUGAError, msg
-
-        # the depth is in meters, and it is the distance from the ocean
-        # to the sea bottom.
-        lonlatdep = p_array.array('f')
-        lonlatdep.read(mux_file, columns * self.points_num)
-        lonlatdep = array(lonlatdep, typecode=Float)    
-        lonlatdep = reshape(lonlatdep, (self.points_num, columns))
-        #print 'lonlatdep',lonlatdep
-        self.lonlatdep = lonlatdep
-        
-        self.mux_file = mux_file
-        # check this array
-
-    def __iter__(self):
-        """
-        iterate over quantity data which is with respect to time.
-
-        Note: You can only interate once over an object
-        
-        returns quantity infomation for each time slice 
-        """
-        msg =  "You can only interate once over a urs file."
-        assert not self.iterated, msg
-        self.time_step = 0
-        self.iterated = True
-        return self
-    
-    def next(self):
-        if self.time_step_count == self.time_step:
-            self.close()
-            raise StopIteration
-        #Read in a time slice  from mux file  
-        hz_p_array = p_array.array('f')
-        hz_p_array.read(self.mux_file, self.points_num)
-        hz_p = array(hz_p_array, typecode=Float)
-        self.time_step += 1
-        
-        return hz_p
-
-    def close(self):
-        self.mux_file.close()
+    
+    outfile.variables['time'][:] = times   
 
 
 def write_sww_triangulation(outfile, points_utm, volumes,
-                            elevation, zone, verbose=False):
-
-    number_of_points = len(points_utm)
-    #Store
-    x = zeros(number_of_points, Float)  #Easting
-    y = zeros(number_of_points, Float)  #Northing
-
+                            elevation, zone, origin=None, verbose=False):
+
+    number_of_points = len(points_utm)   
     volumes = array(volumes)
 
-    geo_ref = Geo_reference(zone,min(x),min(y))
+    if origin is not None:
+        if isinstance(origin, Geo_reference): 
+            geo_ref = origin
+        else:
+            geo_ref = apply(Geo_reference,origin)
+    else:
+        geo_ref = Geo_reference(zone,min(points_utm[:,0]),min(points_utm[:,1]))
+    #geo_ref = Geo_reference(zone,0,0)
     geo_ref.write_NetCDF(outfile)
 
@@ -4773 +4728 @@
 
     #z = resize(bath_grid,outfile.variables['z'][:].shape)
-    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
-    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
+    outfile.variables['x'][:] = points_utm[:,0] - geo_ref.get_xllcorner()
+    outfile.variables['y'][:] = points_utm[:,1] - geo_ref.get_yllcorner()
     outfile.variables['z'][:] = elevation
     outfile.variables['elevation'][:] = elevation  #FIXME HACK
@@ -4782 +4737 @@
     # FIXME: Don't store the whole speeds and stage in  memory.
     # block it?
-def write_sww_time_slice(outfile, has, uas, vas, elevation,
+def write_sww_time_slices(outfile, has, uas, vas, elevation,
                          mean_stage=0, zscale=1,
                          verbose=False):    
@@ -4792 +4747 @@
     if verbose: print 'Converting quantities'
     n = len(has)
-    j = 0
-    
+    j=0
     for ha, ua, va in map(None, has, uas, vas):
         if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
@@ -4802 +4756 @@
         ymomentum[j] = va*h
         j += 1
-    outfile.close()
-    
+
+ 
+def write_sww_time_slice(outfile, ha, ua, va, elevation, slice_index,
+                         mean_stage=0, zscale=1,
+                         verbose=False):    
+    #Time stepping
+    stage = outfile.variables['stage']
+    xmomentum = outfile.variables['xmomentum']
+    ymomentum = outfile.variables['ymomentum']
+
+    if verbose: print 'Converting quantities'
+    
+    w = zscale*ha + mean_stage
+    stage[slice_index] = w
+    h = w - elevation
+    xmomentum[slice_index] = ua*h
+    ymomentum[slice_index] = va*h
+        
+class Urs_points:
+    """
+    Read the info in URS mux files.
+
+    for the quantities heres a correlation between the file names and
+    what they mean;
+    z-mux is height above sea level, m
+    e-mux is velocity is Eastern direction, m/s
+    n-mux is velocity is Northern direction, m/s   
+    """
+    def __init__(self,urs_file):
+        self.iterated = False
+        columns = 3 # long, lat , depth
+        mux_file = open(urs_file, 'rb')
+        
+        # Number of points/stations
+        (self.points_num,)= unpack('i',mux_file.read(4))
+        
+        # nt, int - Number of time steps
+        (self.time_step_count,)= unpack('i',mux_file.read(4))
+        #print "self.time_step_count", self.time_step_count 
+        #dt, float - time step, seconds
+        (self.time_step,) = unpack('f', mux_file.read(4))
+        #print "self.time_step", self.time_step
+        msg = "Bad data in the urs file."
+        if self.points_num < 0:
+            mux_file.close()
+            raise ANUGAError, msg
+        if self.time_step_count < 0:
+            mux_file.close()
+            raise ANUGAError, msg
+        if self.time_step < 0:
+            mux_file.close()
+            raise ANUGAError, msg
+
+        # the depth is in meters, and it is the distance from the ocean
+        # to the sea bottom.
+        lonlatdep = p_array.array('f')
+        lonlatdep.read(mux_file, columns * self.points_num)
+        lonlatdep = array(lonlatdep, typecode=Float)    
+        lonlatdep = reshape(lonlatdep, (self.points_num, columns))
+        #print 'lonlatdep',lonlatdep
+        self.lonlatdep = lonlatdep
+        
+        self.mux_file = mux_file
+        # check this array
+
+    def __iter__(self):
+        """
+        iterate over quantity data which is with respect to time.
+
+        Note: You can only interate once over an object
+        
+        returns quantity infomation for each time slice 
+        """
+        msg =  "You can only interate once over a urs file."
+        assert not self.iterated, msg
+        self.iter_time_step = 0
+        self.iterated = True
+        return self
+    
+    def next(self):
+        if self.time_step_count == self.iter_time_step:
+            self.close()
+            raise StopIteration
+        #Read in a time slice  from mux file  
+        hz_p_array = p_array.array('f')
+        hz_p_array.read(self.mux_file, self.points_num)
+        hz_p = array(hz_p_array, typecode=Float)
+        self.iter_time_step += 1
+        
+        return hz_p
+
+    def close(self):
+        self.mux_file.close()
+
     #### END URS UNGRIDDED 2 SWW ###
 

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -5221 +5221 @@
         time_step_count = 2
         time_step = 400
+        tide = 9000000
         base_name, files = self.write_mux(lat_long,
                                           time_step_count, time_step)
-        urs_ungridded2sww(base_name)
-   
-        
+        urs_ungridded2sww(base_name, mean_stage=tide)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+        
+        assert allclose(ensure_numeric(times),
+                        ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        assert allclose(answer, actual)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+  
+    def test_urs_ungridded2swwII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 2
+        time_step = 400
+        tide = 9000000
+        geo_reference = Geo_reference(50, 3434543,34534543)
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs_ungridded2sww(base_name, mean_stage=tide, origin = geo_reference)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+        
+        assert allclose(ensure_numeric(times),
+                        ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        assert allclose(answer, actual)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        #os.remove(sww_file)
+  
+    def test_urs_ungridded2swwIII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 2
+        time_step = 400
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343))
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        y = points[:,1]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+        
+        assert allclose(ensure_numeric(times),
+                        ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        assert allclose(answer, actual)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_URS_points_needed_and_urs_ungridded2sww(self):
+        # This doesn't actually check anything
+        #  
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+        time_step_count = 2
+        time_step = 400
+        tide = -200000
+
+        boundary_polygon = [[250000,7660000],[280000,7660000],
+                             [280000,7630000],[250000,7630000]]
+        geo=URS_points_needed(boundary_polygon,
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount)
+        lat_long = geo.get_data_points(as_lat_long=True)
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs_ungridded2sww(base_name, mean_stage=tide)
+        self.delete_mux(files)
+        os.remove( base_name + '.sww')
+        
+    def visual_test_URS_points_needed_and_urs_ungridded2sww(self):
+        
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+        time_step_count = 2
+        time_step = 400
+        tide = -200000
+
+        boundary_polygon = [[250000,7660000],[270000,7650000],
+                             [280000,7630000],[250000,7630000]]
+        geo=URS_points_needed(boundary_polygon,
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount)
+        lat_long = geo.get_data_points(as_lat_long=True)
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs_ungridded2sww(base_name, mean_stage=tide)
+        self.delete_mux(files)
+        #os.remove( base_name + '.sww')
 #-------------------------------------------------------------
 if __name__ == "__main__":
     #suite = unittest.makeSuite(Test_Data_Manager,'test_URS_points_needed')
-    #suite = unittest.makeSuite(Test_Data_Manager,'dave_test_URS_poinneeded')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs_ungridded2swwII')
+    #suite = unittest.makeSuite(Test_Data_Manager,'visual_test_URS_points_needed_and_urs_ungridded2sww')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     runner = unittest.TextTestRunner()