Changeset 5462

Timestamp:

Jul 4, 2008, 11:48:37 AM (15 years ago)

Author:

ole

Message:

John and Ole implemented urs2sts using an ordering file creating smaller and pre-ordered sts files for use with boundary conditions.

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (10 diffs)
• test_data_manager.py (modified) (4 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4903 +4903 @@
 def urs2sts(basename_in, basename_out = None, weights=None,
             verbose = False, origin = None,
-            mean_stage=0.0,zscale=1.0,
+            mean_stage=0.0, zscale=1.0,
+            ordering_filename = None,
             minlat = None, maxlat = None,
             minlon = None, maxlon = None):
@@ -4913 +4914 @@
     origin is a 3-tuple with geo referenced
     UTM coordinates (zone, easting, northing)
+    
+    input:
+    basename_in: list of source file prefixes
+    
+    These are combined with the extensions:
+    WAVEHEIGHT_MUX2_LABEL = '_waveheight-z-mux2' for stage
+    EAST_VELOCITY_MUX2_LABEL =  '_velocity-e-mux2' xmomentum
+    NORTH_VELOCITY_MUX2_LABEL =  '_velocity-n-mux2' and ymomentum
+    
+    to create a 2D list of mux2 file. The rows are associated with each quantity and must have the above extensions
+    the columns are the list of file prefixes.
+    
+    ordering: a .txt file name specifying which mux2 gauge points are sto be stored. This is indicated by the index of
+    the gauge in the ordering file.
+    
+    ordering file format:
+    header
+    index,longitude,latitude
+    
+    header='index,longitude,latitude\n'
+    If ordering is None all points are taken in the order they appear in the mux2 file subject to rectangular clipping box (see below).
+    
+    
+    minlat, minlon, maxlat, maxlon define a rectangular clipping box and can be used to extract gauge information in a rectangular box. This is useful for extracting
+    information at gauge points not on the boundary.
+    if ordering_filename is specified clipping box cannot be used
+    
+    
+    output:
+    baename_out: name of sts file in which mux2 data is stored.
+    
+    
     """
     import os
@@ -4944 +4977 @@
     precision = Float
 
+    if ordering_filename is not None:
+        msg = 'If ordering_filename is specified,'
+        msg += ' rectangular clipping box cannot be specified as well. '
+        assert minlat is None and maxlat is None and\
+            minlon is None and maxlon is None, msg
+    
+    
     msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
-
     if minlat != None:
         assert -90 < minlat < 90 , msg
@@ -4959 +4998 @@
             assert maxlon > minlon
 
+
+    # Check output filename    
     if basename_out is None:
-        stsname = basename_in[0] + '.sts'
-    else:
-        stsname = basename_out + '.sts'
+        msg = 'STS filename must be specified'
+        raise Exception, msg
+    
+    stsname = basename_out + '.sts'
 
     files_in=[[],[],[]]
@@ -4991 +5033 @@
     mux={}
     for i,quantity in enumerate(quantities):
-        times_urs, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights)
+        # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
+        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights)
         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(times,times_old),msg
             assert allclose(latitudes_urs,latitudes_urs_old),msg
             assert allclose(longitudes_urs,longitudes_urs_old),msg
             assert allclose(elevation,elevation_old),msg
             assert allclose(starttime,starttime_old)
-        times_urs_old=times_urs
+        times_old=times
         latitudes_urs_old=latitudes_urs
         longitudes_urs_old=longitudes_urs
         elevation_old=elevation
         starttime_old=starttime
-        
-    if (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
+
+        
+    if ordering_filename is not None:
+        # Read ordering file 
+        
+        try:
+            fid=open(ordering_filename, 'r')
+            file_header=fid.readline().split(',')
+            ordering_lines=fid.readlines()
+            fid.close()
+        except:
+            msg = 'Cannot open %s'%ordering_filename
+            raise Exception, msg
+
+        reference_header = 'index, longitude, latitude\n'.split(',')
+        for i in range(3):
+            if not file_header[i] == reference_header[i]:
+                msg = 'File must contain header\n'+header+'\n'
+                raise Exception, msg
+
+        if len(ordering_lines)<2:
+            msg = 'File must contain at least two points'
+            raise Exception, msg
+
+   
+        
+        permutation = [int(line.split(',')[0]) for line in ordering_lines]
+   
+        latitudes = take(latitudes_urs, permutation)
+        longitudes = take(longitudes_urs, permutation)
+        
+        # Self check - can be removed to improve speed
+        ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines]                
+        ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines]        
+        
+        msg = 'Longitudes specified in ordering file do not match those found in mux files'
+        msg += 'I got %s instead of %s (only beginning shown)'\
+            %(str(longitudes[:10]) + '...',
+              str(ref_longitudes[:10]) + '...')        
+        assert allclose(longitudes, ref_longitudes), msg        
+        
+        msg = 'Latitudes specified in ordering file do not match those found in mux files. '
+        msg += 'I got %s instead of %s (only beginning shown)'\
+            %(str(latitudes[:10]) + '...',
+              str(ref_latitudes[:10]) + '...')                
+        assert allclose(latitudes, ref_latitudes), msg
+        
+    elif (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None):
+        #FIXME - this branch is probably not working    
         if verbose: print 'Cliping urs data'
         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)
+        
+        permutation = range(latitudes.shape[0])        
+        
+        msg = 'Clipping is not done yet'
+        raise Exception, msg
+        
     else:
         latitudes=latitudes_urs
         longitudes=longitudes_urs
-        times=times_urs
+        permutation = range(latitudes.shape[0])                
+        
 
     msg='File is empty and or clipped region not in file region'
@@ -5032 +5128 @@
     #starttime = times[0]
     sts = Write_sts()
-    sts.store_header(outfile, times+starttime,number_of_points, description=description,
-                     verbose=verbose,sts_precision=Float)
+    sts.store_header(outfile, 
+                     times+starttime,
+                     number_of_points, 
+                     description=description,
+                     verbose=verbose,
+                     sts_precision=Float)
     
     # Store
@@ -5067 +5167 @@
     stage = outfile.variables['stage']
     xmomentum = outfile.variables['xmomentum']
-    ymomentum =outfile.variables['ymomentum']
+    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]
+        for i, index in enumerate(permutation):
+            w = zscale*mux['HA'][index,j] + mean_stage
+            h=w-elevation[index]
             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
+
+            xmomentum[j,i] = mux['UA'][index,j]*h
+            ymomentum[j,i] = mux['VA'][index,j]*h
 
     outfile.close()
@@ -5119 +5219 @@
     xllcorner = fid.xllcorner[0]
     yllcorner = fid.yllcorner[0]
-    #Points stored in sts file are normailsed to [xllcorner,yllcorner] but 
+    #Points stored in sts file are normalised to [xllcorner,yllcorner] but 
     #we cannot assume that boundary polygon will be. At least the
     #additional points specified by the user after this function is called
@@ -5132 +5232 @@
     yllcorner = fid.yllcorner[0]
 
+    # Walk through ordering file
     boundary_polygon=[]
     for line in lines:
@@ -5139 +5240 @@
             zone,easting,northing=redfearn(float(fields[1]),float(fields[2]))
             if not zone == fid.zone[0]:
-                msg = 'Inconsitent zones'
+                msg = 'Inconsistent zones'
                 raise Exception, msg
         else:
             easting = float(fields[1])
             northing = float(fields[2])
+            
         if not allclose(array([easting,northing]),sts_points[index]):
             msg = "Points in sts file do not match those in the"+\
-                ".txt file spcifying the order of the boundary points"
+                ".txt file specifying the order of the boundary points"
             raise Exception, msg
         boundary_polygon.append(sts_points[index].tolist())

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -6319 +6319 @@
                                       va=va)
 
-        urs2sts(base_name,mean_stage=tide,verbose=False)
+        urs2sts(base_name,
+                basename_out=base_name, 
+                mean_stage=tide,verbose=False)
 
         # now I want to check the sts file ...
@@ -6441 +6443 @@
 
         # Call urs2sts with multiple mux files
-        urs2sts([base_nameI, base_nameII], weights=[1.0,1.0],
+        urs2sts([base_nameI, base_nameII], 
+                basename_out=base_nameI, 
+                weights=[1.0,1.0],
                 mean_stage=tide,
                 verbose=False)
@@ -6534 +6538 @@
         self.delete_mux(filesII)
         os.remove(sts_file)
+        
+        
+    def test_urs2sts_ordering(self):
+        """Test multiple sources with ordering file
+        """
+        
+        tide=0
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        n=len(lat_long_points)
+        first_tstep=ones(n,Int)
+        first_tstep[0]+=1
+        first_tstep[2]+=1
+        last_tstep=(time_step_count)*ones(n,Int)
+        last_tstep[0]-=1
+
+        gauge_depth=20*ones(n,Float)
+        ha=2*ones((n,time_step_count),Float)
+        ha[0]=arange(0,time_step_count)
+        ha[1]=arange(time_step_count,2*time_step_count)
+        ha[2]=arange(2*time_step_count,3*time_step_count)
+        ha[3]=arange(3*time_step_count,4*time_step_count)
+        ua=5*ones((n,time_step_count),Float)
+        va=-10*ones((n,time_step_count),Float)
+
+        # Create two identical mux files to be combined by urs2sts
+        base_nameI, filesI = self.write_mux2(lat_long_points,
+                                             time_step_count, time_step,
+                                             first_tstep, last_tstep,
+                                             depth=gauge_depth,
+                                             ha=ha,
+                                             ua=ua,
+                                             va=va)
+
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha,
+                                               ua=ua,
+                                               va=va)
+
+                                               
+        # Create ordering file
+        permutation = [3,0,2]
+
+        # Do it wrongly at first and check that exception is being raised
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][0], 
+                                             lat_long_points[index][1]))
+        order_fid.close()
+        
+        try:
+            urs2sts([base_nameI, base_nameII], 
+                    basename_out=base_nameI, 
+                    ordering_filename=ordering_filename,
+                    weights=[1.0,1.0],
+                    mean_stage=tide,
+                    verbose=False)  
+            os.remove(ordering_filename)            
+        except:
+            pass
+        else:
+            msg = 'Should have caught wrong lat longs'
+            raise Exception, msg
+
+
+        
+            
+        # Then do it right
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][1], 
+                                             lat_long_points[index][0]))
+        order_fid.close()
+        
+            
+
+                                               
+        # Call urs2sts with multiple mux files
+        urs2sts([base_nameI, base_nameII], 
+                basename_out=base_nameI, 
+                ordering_filename=ordering_filename,
+                weights=[1.0,1.0],
+                mean_stage=tide,
+                verbose=False)
+
+        # now I want to check the sts file ...
+        sts_file = base_nameI + '.sts'
+
+        #Let's interrogate the sts file
+        # Note, the sts 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]
+
+        for i, index in enumerate(permutation):
+            # Check that STS points are stored in the correct order
+            
+            # Work out the UTM coordinates sts point i
+            zone, e, n = redfearn(lat_long_points[index][0], 
+                                  lat_long_points[index][1])             
+            
+            assert allclose([x[i],y[i]], [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 sts values
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        # Set original data used to write mux file to be zero when gauges are
+        # not recdoring
+        
+        ha[0][0]=0.0
+        ha[0][time_step_count-1]=0.0
+        ha[2][0]=0.0
+        ua[0][0]=0.0
+        ua[0][time_step_count-1]=0.0
+        ua[2][0]=0.0
+        va[0][0]=0.0
+        va[0][time_step_count-1]=0.0
+        va[2][0]=0.0;
+
+        
+        # The stage stored in the .sts file should be the sum of the stage
+        # in the two mux2 files because both have weights = 1. In this case
+        # the mux2 files are the same so stage == 2.0 * ha
+        #print 2.0*transpose(ha) - stage 
+        
+        ha_permutation = take(ha, permutation) 
+        ua_permutation = take(ua, permutation)         
+        va_permutation = take(va, permutation)                 
+        gauge_depth_permutation = take(gauge_depth, permutation)                         
+
+        
+        assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth=zeros((len(lat_long_points),time_step_count),Float)
+        for i in range(len(lat_long_points)):
+            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
+            #2.0*ha necessary because using two files with weights=1 are used
+            
+        depth_permutation = take(depth, permutation)                     
+        
+
+        # The xmomentum stored in the .sts file should be the sum of the ua
+        # in the two mux2 files multiplied by the depth.
+        assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        # The ymomentum stored in the .sts file should be the sum of the va
+        # in the two mux2 files multiplied by the depth.
+        assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert allclose(-gauge_depth_permutation, elevation)  #Meters
+
+        fid.close()
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
+        os.remove(sts_file)
+        
 
     def test_file_boundary_stsI(self):
@@ -6563 +6769 @@
 
         sts_file=base_name
-        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=tide,verbose=False)
         self.delete_mux(files)