Changeset 5347

Timestamp:

May 21, 2008, 9:52:16 AM (16 years ago)

Author:

jakeman

Message:

added sts format to datamanager and mux2sts

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (5 diffs)
• test_data_manager.py (modified) (3 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -2697 +2697 @@
     latitudes = file_h.variables[dim_h_latitude]
     longitudes = file_h.variables[dim_h_longitude]
-    
+
+    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes[:],
+                                                  longitudes[:],
+                                                  minlat, maxlat,
+                                                  minlon, maxlon)
     # get dimensions for file_e
     for dimension in file_e.dimensions.keys():
@@ -3912 +3916 @@
     latitudes = ensure_numeric(latitudes)
     longitudes = ensure_numeric(longitudes)
-    
+
     assert allclose(sort(longitudes), longitudes)
+
+    print latitudes[0],longitudes[0]
+    print len(latitudes),len(longitudes)
+    print latitudes[len(latitudes)-1],longitudes[len(longitudes)-1]
     
     lat_ascending = True
@@ -4816 +4824 @@
     # Do some conversions while writing the sww file
 
-    
+    ##################################
+    # READ MUX2 FILES line of points #
+    ##################################
+
+WAVEHEIGHT_MUX2_LABEL = '_waveheight-z-mux2'
+EAST_VELOCITY_MUX2_LABEL =  '_velocity-e-mux2'
+NORTH_VELOCITY_MUX2_LABEL =  '_velocity-n-mux2'    
+
+def read_mux2_py(filenames,weights=None):
+
+    from Numeric import ones,Float,compress,zeros,arange
+    from urs_ext import read_mux2
+
+    if weights is None:
+        weights=ones(len(filenames),Float)/len(filenames) #default 1/numSrc
+
+    file_params=-1*ones(3,Float)#[nsta,dt,nt]
+    write=1 #write txt files to current directory as well
+    data=read_mux2(1,filenames,weights,file_params,write)
+
+    msg='File parameter values were not read in correctly from c file'
+    assert len(compress(file_params>0,file_params))!=0,msg
+    msg='The number of stations specifed in the c array and in the file are inconsitent'
+    assert file_params[0]==data.shape[0],msg
+    
+    nsta=int(file_params[0])
+    msg='Must have at least one station'
+    assert nsta>0,msg
+    dt=file_params[1]
+    msg='Must have a postive timestep'
+    assert dt>0,msg
+    nt=int(file_params[2])
+    msg='Must have at least one gauge value'
+    assert nt>0,msg
+    
+    OFFSET=5 #number of site parameters p passed back with data
+    #p=[geolat,geolon,depth,start_tstep,finish_tstep]
+
+    times=dt*arange(1,(data.shape[1]-OFFSET)+1)
+    latitudes=zeros(data.shape[0],Float)
+    longitudes=zeros(data.shape[0],Float)
+    elevation=zeros(data.shape[0],Float)
+    stage=zeros((data.shape[0],data.shape[1]-OFFSET),Float)
+    for i in range(0,data.shape[0]):
+        latitudes[i]=data[i][data.shape[1]-OFFSET]
+        longitudes[i]=data[i][data.shape[1]-OFFSET+1]
+        elevation[i]=-data[i][data.shape[1]-OFFSET+2]
+        stage[i]=data[i][:-OFFSET]
+
+    return times, latitudes, longitudes, elevation, stage
+
 def mux2sww_time(mux_times, mint, maxt):
     """
@@ -4834 +4892 @@
 
     return mux_times_start_i, mux_times_fin_i
+
+
+def urs2sts(basename_in, basename_out = None, verbose = False, origin = None,
+            mean_stage=0.0,zscale=1.0,
+            minlat = None, maxlat = None,
+            minlon = None, maxlon = None):
+    """Convert URS mux2 format for wave propagation to sts format
+
+    Specify only basename_in and read files of the form
+    out_waveheight-z-mux2
+
+    Also convert latitude and longitude to UTM. All coordinates are
+    assumed to be given in the GDA94 datum
+
+    origin is a 3-tuple with geo referenced
+    UTM coordinates (zone, easting, northing)
+    """
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, Int, Int32, searchsorted, zeros, array
+    from Numeric import allclose, around
+
+    precision = Float
+
+    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
+
+    if minlat != None:
+        assert -90 < minlat < 90 , msg
+    if maxlat != None:
+        assert -90 < maxlat < 90 , msg
+        if minlat != None:
+            assert maxlat > minlat
+    if minlon != None:
+        assert -180 < minlon < 180 , msg
+    if maxlon != None:
+        assert -180 < maxlon < 180 , msg
+        if minlon != None:
+            assert maxlon > minlon
+
+    if basename_out is None:
+        stsname = basename_in + '.sts'
+    else:
+        stsname = basename_out + '.sts'
+
+    files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL,
+                basename_in + EAST_VELOCITY_MUX2_LABEL,
+                basename_in + NORTH_VELOCITY_MUX2_LABEL]
+    quantities = ['HA','UA','VA']
+
+    #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
+    #for now set x_momentum and y_moentum quantities to zero
+    if (verbose): print 'reading mux2 file'
+    mux={}
+    for quantity, file in map(None, quantities, files_in):
+        times_urs, latitudes_urs, longitudes_urs, elevation, mux[quantity] = read_mux2_py([file])
+        if quantity!=quantities[0]:
+            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
+            assert allclose(times_urs,times_urs_old),msg
+            assert allclose(latitudes_urs,latitudes_urs_old),msg
+            assert allclose(longitudes_urs,longitudes_urs_old),msg
+            assert allclose(elevation,elevation_old),msg
+        times_urs_old=times_urs
+        latitudes_urs_old=latitudes_urs
+        longitudes_urs_old=longitudes_urs
+        elevation_old=elevation
+        
+    if (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
+        latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs)
+        longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs)
+        times = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),times_urs)
+    else:
+        latitudes=latitudes_urs
+        longitudes=longitudes_urs
+        times=times_urs
+
+    number_of_points = latitudes.shape[0]
+    number_of_times = times.shape[0]
+    number_of_latitudes = latitudes.shape[0]
+    number_of_longitudes = longitudes.shape[0]
+
+    # NetCDF file definition
+    outfile = NetCDFFile(stsname, 'w')
+
+    description = 'Converted from URS mux2 files: %s'\
+                  %(basename_in)
+    
+    # Create new file
+    starttime = times[0]
+    sts = Write_sts()
+    sts.store_header(outfile, times,number_of_points, description=description,
+                     verbose=verbose,sts_precision=Float)
+    
+    # Store
+    from anuga.coordinate_transforms.redfearn import redfearn
+    x = zeros(number_of_points, Float)  #Easting
+    y = zeros(number_of_points, Float)  #Northing
+
+    # Check zone boundaries
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0])  
+
+    for i in range(number_of_points):
+        zone, easting, northing = redfearn(latitudes[i],longitudes[i])
+        x[i] = easting
+        y[i] = northing
+        #print zone,easting,northing
+
+    if origin is None:
+        origin = Geo_reference(refzone,min(x),min(y))
+    geo_ref = write_NetCDF_georeference(origin, outfile)
+
+    z = elevation
+    
+    #print geo_ref.get_xllcorner()
+    #print geo_ref.get_yllcorner()
+
+    z = resize(z,outfile.variables['z'][:].shape)
+    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
+    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
+    outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
+    outfile.variables['elevation'][:] = z
+
+    stage = outfile.variables['stage']
+    xmomentum = outfile.variables['xmomentum']
+    ymomentum =outfile.variables['ymomentum']
+
+    if verbose: print 'Converting quantities'
+    for j in range(len(times)):
+        for i in range(number_of_points):
+            w = zscale*mux['HA'][i,j] + mean_stage
+            h=w-elevation[i]
+            stage[j,i] = w
+            #delete following two lines once velcotu files are read in.
+            xmomentum[j,i] = mux['UA'][i,j]*h
+            ymomentum[j,i] = mux['VA'][i,j]*h
+
+    outfile.close()
 
 
@@ -5226 +5420 @@
             
             i +=1
+
+class Write_sts:
+
+    sts_quantities = ['stage','xmomentum','ymomentum']
+
+
+    RANGE = '_range'
+    EXTREMA = ':extrema'
+    
+    def __init__(self):
+        pass
+
+    def store_header(self,
+                     outfile,
+                     times,
+                     number_of_points,
+                     description='Converted from URS mux2 format',
+                     sts_precision=Float32,
+                     verbose=False):
+        """
+        outfile - the name of the file that will be written
+        times - A list of the time slice times OR a start time
+        Note, if a list is given the info will be made relative.
+        number_of_points - the number of urs gauges sites
+        """
+
+        outfile.institution = 'Geoscience Australia'
+        outfile.description = description
+        
+        try:
+            revision_number = get_revision_number()
+        except:
+            revision_number = None
+        # Allow None to be stored as a string                
+        outfile.revision_number = str(revision_number) 
+        
+        # times - A list or array of the time slice times OR a start time
+        # Start time in seconds since the epoch (midnight 1/1/1970)
+
+        # This is being used to seperate one number from a list.
+        # what it is actually doing is sorting lists from numeric arrays.
+        if type(times) is list or type(times) is ArrayType:  
+            number_of_times = len(times)
+            times = ensure_numeric(times)  
+            if number_of_times == 0:
+                starttime = 0
+            else:
+                starttime = times[0]
+                times = times - starttime  #Store relative times
+        else:
+            number_of_times = 0
+            starttime = times
+
+        outfile.starttime = starttime
+
+        # Dimension definitions
+        outfile.createDimension('number_of_points', number_of_points)
+        outfile.createDimension('number_of_timesteps', number_of_times)
+        outfile.createDimension('numbers_in_range', 2)
+
+        # Variable definitions
+        outfile.createVariable('x', sts_precision, ('number_of_points',))
+        outfile.createVariable('y', sts_precision, ('number_of_points',))
+        outfile.createVariable('elevation', sts_precision,('number_of_points',))
+
+        q = 'elevation'
+        outfile.createVariable(q+Write_sts.RANGE, sts_precision,
+                               ('numbers_in_range',))
+
+        # Initialise ranges with small and large sentinels.
+        # If this was in pure Python we could have used None sensibly
+        outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min 
+        outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
+
+        outfile.createVariable('z', sts_precision, ('number_of_points',))
+        # Doing sts_precision instead of Float gives cast errors.
+        outfile.createVariable('time', Float, ('number_of_timesteps',))
+
+        for q in Write_sts.sts_quantities:
+            outfile.createVariable(q, sts_precision,
+                                   ('number_of_timesteps',
+                                    'number_of_points'))
+            outfile.createVariable(q+Write_sts.RANGE, sts_precision,
+                                   ('numbers_in_range',))
+            # Initialise ranges with small and large sentinels.
+            # If this was in pure Python we could have used None sensibly
+            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
+            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
+
+        if type(times) is list or type(times) is ArrayType:  
+            outfile.variables['time'][:] = times    #Store time relative
+
+        if verbose:
+            print '------------------------------------------------'
+            print 'Statistics:'
+            print '    t in [%f, %f], len(t) == %d'\
+                  %(min(times.flat), max(times.flat), len(times.flat))
+
+    def store_points(self,
+                     outfile,
+                     points_utm,
+                     elevation, zone=None, new_origin=None, 
+                     points_georeference=None, verbose=False):
+
+        """
+        points_utm - currently a list or array of the points in UTM.
+        points_georeference - the georeference of the points_utm
+
+        How about passing new_origin and current_origin.
+        If you get both, do a convertion from the old to the new.
+        
+        If you only get new_origin, the points are absolute,
+        convert to relative
+        
+        if you only get the current_origin the points are relative, store
+        as relative.
+        
+        if you get no georefs create a new georef based on the minimums of
+        points_utm.  (Another option would be to default to absolute)
+        
+        Yes, and this is done in another part of the code.
+        Probably geospatial.
+        
+        If you don't supply either geo_refs, then supply a zone. If not
+        the default zone will be used.
+
+        precondition:
+             header has been called.
+        """
+
+        number_of_points = len(points_utm)
+        points_utm = array(points_utm)
+
+        # given the two geo_refs and the points, do the stuff
+        # described in the method header
+        points_georeference = ensure_geo_reference(points_georeference)
+        new_origin = ensure_geo_reference(new_origin)
+        
+        if new_origin is None and points_georeference is not None:
+            points = points_utm
+            geo_ref = points_georeference
+        else:
+            if new_origin is None:
+                new_origin = Geo_reference(zone,min(points_utm[:,0]),
+                                           min(points_utm[:,1]))
+            points = new_origin.change_points_geo_ref(points_utm,
+                                                      points_georeference)
+            geo_ref = new_origin
+
+        # At this stage I need a georef and points
+        # the points are relative to the georef
+        geo_ref.write_NetCDF(outfile)
+
+        x =  points[:,0]
+        y =  points[:,1]
+        z = outfile.variables['z'][:]
+    
+        if verbose:
+            print '------------------------------------------------'
+            print 'More Statistics:'
+            print '  Extent (/lon):'
+            print '    x in [%f, %f], len(lat) == %d'\
+                  %(min(x), max(x),
+                    len(x))
+            print '    y in [%f, %f], len(lon) == %d'\
+                  %(min(y), max(y),
+                    len(y))
+            print '    z in [%f, %f], len(z) == %d'\
+                  %(min(elevation), max(elevation),
+                    len(elevation))
+            print 'geo_ref: ',geo_ref
+            print '------------------------------------------------'
+            
+        #z = resize(bath_grid,outfile.variables['z'][:].shape)
+        #print "points[:,0]", points[:,0]
+        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
+        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
+        outfile.variables['z'][:] = elevation
+        outfile.variables['elevation'][:] = elevation  #FIXME HACK4
+
+        q = 'elevation'
+        # This updates the _range values
+        outfile.variables[q+Write_sts.RANGE][0] = min(elevation)
+        outfile.variables[q+Write_sts.RANGE][1] = max(elevation)
+
+    def store_quantities(self, outfile, sts_precision=Float32,
+                         slice_index=None, time=None,
+                         verbose=False, **quant):
+        
+        """
+        Write the quantity info.
+
+        **quant is extra keyword arguments passed in. These must be
+          the sts quantities, currently; stage.
+        
+        if the time array is already been built, use the slice_index
+        to specify the index.
+        
+        Otherwise, use time to increase the time dimension
+
+        Maybe make this general, but the viewer assumes these quantities,
+        so maybe we don't want it general - unless the viewer is general
+        
+        precondition:
+            triangulation and
+            header have been called.
+        """
+        if time is not None:
+            file_time = outfile.variables['time']
+            slice_index = len(file_time)
+            file_time[slice_index] = time    
+
+        # Write the conserved quantities from Domain.
+        # Typically stage,  xmomentum, ymomentum
+        # other quantities will be ignored, silently.
+        # Also write the ranges: stage_range
+        for q in Write_sts.sts_quantities:
+            if not quant.has_key(q):
+                msg = 'STS file can not write quantity %s' %q
+                raise NewQuantity, msg
+            else:
+                q_values = quant[q]
+                outfile.variables[q][slice_index] = \
+                                q_values.astype(sts_precision)
+
+                # This updates the _range values
+                q_range = outfile.variables[q+Write_sts.RANGE][:]
+                q_values_min = min(q_values)
+                if q_values_min < q_range[0]:
+                    outfile.variables[q+Write_sts.RANGE][0] = q_values_min
+                q_values_max = max(q_values)
+                if q_values_max > q_range[1]:
+                    outfile.variables[q+Write_sts.RANGE][1] = q_values_max
+
+
     
 class Urs_points:

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -5977 +5977 @@
         os.remove(sww_file)
         
+
+    def write_mux2(self,lat_long_points, time_step_count, time_step,
+                  depth=None, ha=None, ua=None, va=None
+                  ):
+        """
+        This will write 3 non-gridded mux files, for testing.
+        If no quantities are passed in,
+        na and va quantities will be the Easting values.
+        Depth and ua will be the Northing value.
+        """
+        #print "lat_long_points", lat_long_points
+        #print "time_step_count",time_step_count
+        #print "time_step",
+
+        #irrelevant header information
+        ig=ilon=ilat=0
+        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
+
+        points_num = len(lat_long_points)
+        latlondeps = []
+        quantities = ['HA','UA','VA']
+
+        mux_names = [WAVEHEIGHT_MUX2_LABEL,
+                     EAST_VELOCITY_MUX2_LABEL,
+                     NORTH_VELOCITY_MUX2_LABEL]
+        quantities_init = [[],[],[]]
+        # urs binary is latitude fastest
+        for point in lat_long_points:
+            lat = point[0]
+            lon = point[1]
+            _ , e, n = redfearn(lat, lon)
+            if depth is None:
+                this_depth = n
+            else:
+                this_depth = depth
+            if ha is None:
+                this_ha = e
+            else:
+                this_ha = ha
+            if ua is None:
+                this_ua = n
+            else:
+                this_ua = ua
+            if va is None:
+                this_va = e   
+            else:
+                this_va = va         
+            latlondeps.append([lat, lon, this_depth])
+            quantities_init[0].append(this_ha) # HA
+            quantities_init[1].append(this_ua) # UA
+            quantities_init[2].append(this_va) # VA
+
+        file_handle, base_name = tempfile.mkstemp("")
+        os.close(file_handle)
+        os.remove(base_name)
+
+        files = []        
+        for i,q in enumerate(quantities): 
+            quantities_init[i] = ensure_numeric(quantities_init[i])
+            #print "HA_init", HA_init
+            q_time = zeros((time_step_count, points_num), Float)
+            for time in range(time_step_count):
+                q_time[time,:] = quantities_init[i] #* time * 4
+
+            #Write C files
+            columns = 3 # long, lat , depth
+            file = base_name + mux_names[i]
+            #print "base_name file",file 
+            f = open(file, 'wb')
+            files.append(file)
+
+            f.write(pack('i',points_num))
+            #write mux 2 header
+            for latlondep in latlondeps:
+                f.write(pack('f',latlondep[0]))
+                f.write(pack('f',latlondep[1]))
+                f.write(pack('f',mcolat))
+                f.write(pack('f',mcolon))
+                f.write(pack('i',ig))
+                f.write(pack('i',ilon))
+                f.write(pack('i',ilat))
+                f.write(pack('f',latlondep[2]))
+                f.write(pack('f',centerlat))
+                f.write(pack('f',centerlon))
+                f.write(pack('f',offset))
+                f.write(pack('f',az))
+                f.write(pack('f',baz))
+                f.write(pack('f',time_step))
+                f.write(pack('i',time_step_count))
+                for j in range(4): # identifier
+                    f.write(pack('f',id))    
+
+            first_tstep=1
+            last_tstep=time_step_count
+            for latlondep in latlondeps:
+                f.write(pack('i',first_tstep))
+            for latlondep in latlondeps:
+                f.write(pack('i',last_tstep))
+
+            # Write quantity info
+
+            for time in  range(time_step_count):
+                #first timestep always assumed to be zero
+                f.write(pack('f',0.0))
+                for point_i in range(points_num):
+                    f.write(pack('f',q_time[time,point_i]))
+
+            f.close()
+        return base_name, files
+
+    def test_read_mux2_py(self):
+        from Numeric import ones,Float
+        tide = 1
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        depth=20
+        ha=2
+        ua=5
+        va=-10 #-ve added to take into account mux file format where south
+               # is positive.
+        base_name, files = self.write_mux2(lat_long_points,
+                                     time_step_count, time_step,
+                                     depth=depth,
+                                     ha=ha,
+                                     ua=ua,
+                                     va=va)
+
+        weights=ones(1,Float)
+        #ensure that files are indeed mux2 files
+        times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
+        self.delete_mux(files)
+
+        msg='time array has incorrect length'
+        assert times.shape[0]==time_step_count,msg
+        msg = 'time array is incorrect'
+        assert allclose(times,time_step*arange(1,time_step_count+1)),msg
+        msg='Incorrect gauge positions returned'
+        for i,point in enumerate(lat_long_points):
+            assert allclose(latitudes[i],point[0]) and allclose(longitudes[i],point[1]),msg
+        msg='Incorrect gauge depths returned'
+        assert allclose(elevation,-depth*ones(len(lat_long_points),Float)),msg
+        msg='incorrect gauge time series returned'
+        assert allclose(stage,ha*ones((len(lat_long_points),time_step_count),Float))
+
+    def test_urs2sts(self):
+        from Numeric import ones,Float
+        tide = 1
+        time_step_count = 3
+        time_step = 2
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        depth=20
+        ha=2
+        ua=5
+        va=-10 #-ve added to take into account mux file format where south
+               # is positive.
+        base_name, files = self.write_mux2(lat_long,
+                                     time_step_count, time_step,
+                                     depth=depth,
+                                     ha=ha,
+                                     ua=ua,
+                                     va=va)
+
+        urs2sts(base_name,mean_stage=tide,verbose=False)
+
+        # now I want to check the sts file ...
+        sts_file = base_name + '.sts'
+
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sts_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], ha+tide)  #Meters
+
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        answer_x = ua*(ha+tide+depth)
+        actual_x = xmomentum[0,0]
+        #print "answer_x",answer_x
+        #print "actual_x",actual_x 
+        assert allclose(answer_x, actual_x)  #Meters
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        answer_y = va*(ha+tide+depth)
+        actual_y = ymomentum[0,0]
+        #print "answer_y",answer_y
+        #print "actual_y",actual_y 
+        assert allclose(answer_y, actual_y)  #Meters
+
+        # check the stage values, first time step.
+        assert allclose(stage[0], ha +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, depth)  #Meters
+
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sts_file)
+
     def test_lon_lat2grid(self):
         lonlatdep = [
@@ -7555 +7795 @@
 
 
-
         
         domain.set_quantity('elevation', 0.0)
@@ -7575 +7814 @@
                           interpolation_points=I,
                           verbose=False)
-
         for t in range(t_end+1):
             for i in range(3):