Changeset 5412

Timestamp:

Jun 23, 2008, 4:26:45 PM (15 years ago)

Author:

jakeman

Message:

allow urs2sts to accept multiple sources

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (10 diffs)
• test_data_manager.py (modified) (10 diffs)
• urs_ext.c (modified) (6 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4838 +4838 @@
 NORTH_VELOCITY_MUX2_LABEL =  '_velocity-n-mux2'    
 
-def read_mux2_py(filenames,weights=None):
+def read_mux2_py(filenames,weights):
 
     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
+    numSrc=len(filenames)
 
     file_params=-1*ones(3,Float)#[nsta,dt,nt]
     write=0 #if true write txt files to current directory as well
-    data=read_mux2(1,filenames,weights,file_params,write)
+    data=read_mux2(numSrc,filenames,weights,file_params,write)
 
     msg='File parameter values were not read in correctly from c file'
@@ -4873 +4872 @@
     elevation=zeros(data.shape[0],Float)
     quantity=zeros((data.shape[0],data.shape[1]-OFFSET),Float)
+    starttime=1e16
     for i in range(0,data.shape[0]):
         latitudes[i]=data[i][data.shape[1]-OFFSET]
@@ -4878 +4878 @@
         elevation[i]=-data[i][data.shape[1]-OFFSET+2]
         quantity[i]=data[i][:-OFFSET]
-
-    return times, latitudes, longitudes, elevation, quantity
+        starttime=min(dt*data[i][data.shape[1]-OFFSET+3],starttime)
+        
+    return times, latitudes, longitudes, elevation, quantity, starttime
 
 def mux2sww_time(mux_times, mint, maxt):
@@ -4900 +4901 @@
 
 
-def urs2sts(basename_in, basename_out = None, verbose = False, origin = None,
+def urs2sts(basename_in, basename_out = None, weights=None,
+            verbose = False, origin = None,
             mean_stage=0.0,zscale=1.0,
             minlat = None, maxlat = None,
@@ -4906 +4908 @@
     """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
@@ -4918 +4917 @@
     from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import Float, Int, Int32, searchsorted, zeros, array
-    from Numeric import allclose, around
+    from Numeric import allclose, around,ones,Float
+    from types import ListType,StringType
+    from operator import __and__
+    
+    if not isinstance(basename_in, ListType):
+        if verbose: print 'Reading single source'
+        basename_in=[basename_in]
+
+    # Check that basename is a list of strings
+    if not reduce(__and__, map(lambda z:isinstance(z,StringType),basename_in)):
+        msg= 'basename_in must be a string or list of strings'
+        raise Exception, msg
+
+    # Find the number of sources to be used
+    numSrc=len(basename_in)
+
+    # A weight must be specified for each source
+    if weights is None:
+        weights=ones(numSrc,Float)/numSrc
+    else:
+        weights = ensure_numeric(weights)
+        msg = 'When combining multiple sources must specify a weight for '+\
+              'mux2 source file'
+        assert len(weights)== numSrc, msg
 
     precision = Float
@@ -4938 +4960 @@
 
     if basename_out is None:
-        stsname = basename_in + '.sts'
+        stsname = basename_in[0] + '.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]
+    files_in=[[],[],[]]
+    for files in basename_in:
+        files_in[0].append(files + WAVEHEIGHT_MUX2_LABEL),
+        files_in[1].append(files + EAST_VELOCITY_MUX2_LABEL)
+        files_in[2].append(files + NORTH_VELOCITY_MUX2_LABEL)
+        
+    #files_in = [basename_in + WAVEHEIGHT_MUX2_LABEL,
+    #            basename_in + EAST_VELOCITY_MUX2_LABEL,
+    #            basename_in + NORTH_VELOCITY_MUX2_LABEL]
+    
     quantities = ['HA','UA','VA']
 
-    for file_in in files_in:
-        if (os.access(file_in, os.F_OK) == 0):
-            msg = 'File %s does not exist or is not accessible' %file_in
-            raise IOError, msg
+    # For each source file check that there exists three files ending with:
+    # WAVEHEIGHT_MUX2_LABEL,
+    # EAST_VELOCITY_MUX2_LABEL, and 
+    # NORTH_VELOCITY_MUX2_LABEL
+    for i in range(len(quantities)): 
+        for file_in in files_in[i]:
+            if (os.access(file_in, os.F_OK) == 0):
+                msg = 'File %s does not exist or is not accessible' %file_in
+                raise IOError, msg
 
     #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
@@ -4956 +4990 @@
     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])
+    for i,quantity in enumerate(quantities):
+        times_urs, 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])
@@ -4964 +4998 @@
             assert allclose(longitudes_urs,longitudes_urs_old),msg
             assert allclose(elevation,elevation_old),msg
+            assert allclose(starttime,starttime_old)
         times_urs_old=times_urs
         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 verbose: print 'Cliiping urs data'
+        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)
@@ -4994 +5030 @@
     
     # Create new file
-    starttime = times[0]
+    #starttime = times[0]
     sts = Write_sts()
-    sts.store_header(outfile, times,number_of_points, description=description,
+    sts.store_header(outfile, times+starttime,number_of_points, description=description,
                      verbose=verbose,sts_precision=Float)
     

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -6106 +6106 @@
         return base_name, files
 
-    def test_read_mux2_py(self):
+    def test_read_mux2_pyI(self):
         """Constant stage,momentum at each gauge
         """
@@ -6131 +6131 @@
         weights=ones(1,Float)
         #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity=read_mux2_py([files[1]],weights)
+        times, latitudes, longitudes, elevation, stage, starttime =read_mux2_py([files[0]],weights)
+        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
         msg='ha and ua have different gauge meta data'
-        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation),msg
-        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity=read_mux2_py([files[2]],weights)
+        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
+        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity, starttime_va=read_mux2_py([files[2]],weights)
         msg='ha and va have different gauge meta data'
-        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation),msg
+        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation) and allclose(starttime,starttime_va), msg
 
         self.delete_mux(files)
@@ -6158 +6158 @@
         assert allclose(yvelocity,va)
 
-    def test_read_mux2_py2(self):
+    def test_read_mux2_pyII(self):
         """Spatially varing stage
         """
@@ -6188 +6188 @@
         weights=ones(1,Float)
         #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity=read_mux2_py([files[1]],weights)
+        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
+        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
         msg='ha and ua have different gauge meta data'
-        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation),msg
-        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity=read_mux2_py([files[2]],weights)
+        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
+        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
         msg='ha and va have different gauge meta data'
-        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation),msg
+        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation) and allclose(starttime,starttime_va), msg
 
 
@@ -6216 +6216 @@
         assert allclose(yvelocity,va)
 
-    def test_read_mux2_py3(self):
+    def test_read_mux2_pyIII(self):
         """Varying start and finsh times
         """
@@ -6249 +6249 @@
         weights=ones(1,Float)
         #ensure that files are indeed mux2 files
-        times, latitudes, longitudes, elevation, stage=read_mux2_py([files[0]],weights)
-        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity=read_mux2_py([files[1]],weights)
+        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
+        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
         msg='ha and ua have different gauge meta data'
-        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation),msg
-        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity=read_mux2_py([files[2]],weights)
+        assert allclose(times,ua_times) and allclose(latitudes,ua_latitudes) and allclose(longitudes,ua_longitudes) and allclose(elevation,ua_elevation) and allclose(starttime,starttime_ua), msg
+        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
         msg='ha and va have different gauge meta data'
-        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation),msg
+        assert allclose(times,va_times) and allclose(latitudes,va_latitudes) and allclose(longitudes,va_longitudes) and allclose(elevation,va_elevation) and allclose(starttime,starttime_va), msg
 
         self.delete_mux(files)
@@ -6288 +6288 @@
         assert allclose(yvelocity,va)
 
-    def test_urs2sts(self):
+    def test_urs2stsI(self):
+        """
+        Test single source
+        """
         tide=0
         time_step_count = 3
@@ -6308 +6311 @@
         ua=5*ones((n,time_step_count),Float)
         va=-10*ones((n,time_step_count),Float)
-        #-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,
@@ -6395 +6398 @@
         fid.close()
         self.delete_mux(files)
+        os.remove(sts_file)
+
+    def test_urs2stsII(self):
+        """
+        Test multiple sources
+        """
+        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)
+    
+        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)
+
+        urs2sts([base_nameI,base_nameII],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 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_points[0][0], lat_long_points[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'][:]
+
+        # 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
+        assert allclose(2.0*transpose(ha),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
+
+        # 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*depth),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*depth),ymomentum)
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert allclose(-elevation, gauge_depth)  #Meters
+
+        fid.close()
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
         os.remove(sts_file)
 
@@ -8414 +8549 @@
     #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
 #    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'test_urs2stsII')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
     #suite = unittest.makeSuite(Test_Data_Manager,'covered_')

anuga_core/source/anuga/shallow_water/urs_ext.c

@@ -29 +29 @@
    
     Python call:
-    read_mux2(filenames,weights,numSrc)
+    read_mux2(numSrc,filenames,weights,file_params,write)
 
     NOTE:
@@ -81 +81 @@
   }
 
-  if(PyList_Size(filenames) != pyweights->nd){
-    PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename");
-    return NULL;
+  if(PyList_Size(filenames) != pyweights->dimensions[0]){
+      PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename");
+      return NULL;
   }
 
@@ -115 +115 @@
   }
 
-  cdata=_read_mux2((int)pyweights->nd,muxFileNameArray,weights,(double*)file_params->data,(int)write);
+  cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)write);
 
 
@@ -200 +200 @@
    /* note starting time */
    time(&start_time);
-   
-   
+
    /* allocate space for the names and the weights and pointers to the data*/    
    wt = (float *) malloc(numSrc*sizeof(float));
@@ -274 +273 @@
       mytgs = (struct tgsrwg *)malloc( nsta0*sizeof(struct tgsrwg) );
    } else {
-     /* FIXME (Ole): What should happen in case the are no source files?*/
+     /* FIXME (JJ): What should happen in case the are no source files?*/
      /* If we exit here, tests will break */
      // fprintf(stderr, "No source file specified\n");
@@ -439 +438 @@
    //free(mytgs0);
    
-   if(numSrc>1)
-   {
+   //if(numSrc>1)
+   //{
      //free(data);
      //free(mytgs);
-   } 
+   //} 
    //can't free arrays because I only fill array by making pointer not copy
    //for(isrc=0; isrc<numSrc;isrc++)