Changeset 5537

Timestamp:

Jul 18, 2008, 4:59:46 PM (16 years ago)

Author:

ole

Message:

Work on permutations in urs_ext.c
It does not work yet, but the array is being passed in and the algorithm
has been described.

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (9 diffs)
• test_data_manager.py (modified) (9 diffs)
• urs_ext.c (modified) (23 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4847 +4847 @@
     """Access the mux files specified in the filenames list. Combine the 
        data found therin as a weighted linear sum as specifed by the weights. 
-       If permutation is None extract timeseries data for all gauges within the 
+       If permutation is None or empty extract timeseries data for all gauges within the 
        files. 
          
@@ -4871 +4871 @@
         verbose=0
         
+        
+    #msg = 'I got the permutation vector:' + str(permutation)
+    #assert permutation is not None, msg
+    if permutation is None:
+        permutation = ensure_numeric([], Float)    
+        
     # Call underlying C implementation urs2sts_ext.c    
-    data=read_mux2(numSrc,filenames,weights,file_params,verbose)
+    data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
 
     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
+    
+    msg='The number of stations specifed in the c array and in the file are inconsistent'
+    assert file_params[0] >= len(permutation)
+    
+    msg='The number of stations returned is inconsistent with the requested number'    
+    assert len(permutation) == 0 or len(permutation) == 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'
@@ -4962 +4974 @@
               other lines: index,longitude,latitude
               
-
-              If ordering is None all points are taken in the order they 
+              If ordering is None or ordering file is empty then 
+               all points are taken in the order they 
               appear in the mux2 file.
     
@@ -5021 +5033 @@
                 raise IOError, msg
 
-    # Read MUX2 files
-    if (verbose): print 'reading mux2 file'
-    mux={}
-    for i, quantity in enumerate(quantities):
-    
-        # For each quantity read the associated list of source mux2 file with 
-        # extention associated with that quantity
-        times,\
-        latitudes_urs,\
-        longitudes_urs,\
-        elevation,\
-        mux[quantity],\
-        starttime = read_mux2_py(files_in[i], weights, verbose=verbose)
-    
-        # Check that all quantities have consistent time and space information     
-        if quantity!=quantities[0]:
-            msg='%s, %s and %s have inconsistent gauge data'\
-                %(files_in[0], files_in[1], files_in[2])
-            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), msg
-        times_old=times
-        latitudes_urs_old=latitudes_urs
-        longitudes_urs_old=longitudes_urs
-        elevation_old=elevation
-        starttime_old=starttime
-
-        
+    # Establish permutation array
     if ordering_filename is not None:
+    
         # Read ordering file 
-        
         try:
             fid=open(ordering_filename, 'r')
@@ -5074 +5057 @@
             raise Exception, msg
 
-   
-        
-        permutation = [int(line.split(',')[0]) for line in ordering_lines]
-   
-        latitudes = take(latitudes_urs, permutation)
-        longitudes = take(longitudes_urs, permutation)
+        permutation = ensure_numeric([int(line.split(',')[0]) for line in ordering_lines])
+    else:
+        # Use empty array to signify 'all' points
+        # We could have used 'None' but it got too hard in the C-code ;-)
+        permutation = ensure_numeric([], Float)
+
+                
+    # Read MUX2 files
+    if (verbose): print 'reading mux2 file'
+    mux={}
+    for i, quantity in enumerate(quantities):
+    
+        # For each quantity read the associated list of source mux2 file with 
+        # extention associated with that quantity
+        times,\
+        latitudes,\
+        longitudes,\
+        elevation,\
+        mux[quantity],\
+        starttime = read_mux2_py(files_in[i], weights, permutation, verbose=verbose)
+    
+        # Check that all quantities have consistent time and space information     
+        if quantity!=quantities[0]:
+            msg='%s, %s and %s have inconsistent gauge data'\
+                %(files_in[0], files_in[1], files_in[2])
+            assert allclose(times, times_old), msg
+            assert allclose(latitudes, latitudes_old), msg
+            assert allclose(longitudes, longitudes_old), msg
+            assert allclose(elevation, elevation_old), msg
+            assert allclose(starttime, starttime_old), msg
+        times_old=times
+        latitudes_old=latitudes
+        longitudes_old=longitudes
+        elevation_old=elevation
+        starttime_old=starttime
+
+
+        
+        
         
         # 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
-        
-    else:
-        latitudes=latitudes_urs
-        longitudes=longitudes_urs
-        permutation = range(latitudes.shape[0])                
-        
-        
+
+        #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
+        
+                
         
     # Store timeseries in NetCDF STS file    
@@ -5108 +5120 @@
     assert len(latitudes>0),msg
 
-    number_of_points = latitudes.shape[0]
-    number_of_times = times.shape[0]
-    number_of_latitudes = latitudes.shape[0]
-    number_of_longitudes = longitudes.shape[0]
+    number_of_points = latitudes.shape[0]      # Number of stations retrieved
+    number_of_times = times.shape[0]           # Number of timesteps
+    number_of_latitudes = latitudes.shape[0]   # Number latitudes
+    number_of_longitudes = longitudes.shape[0] # Number longitudes
 
     # NetCDF file definition
@@ -5120 +5132 @@
     
     # Create new file
-    #starttime = times[0]
     sts = Write_sts()
     sts.store_header(outfile, 
@@ -5131 +5142 @@
     # Store
     from anuga.coordinate_transforms.redfearn import redfearn
-    x = zeros(number_of_points, Float)  #Easting
-    y = zeros(number_of_points, Float)  #Northing
+    x = zeros(number_of_points, Float)  # Easting
+    y = zeros(number_of_points, Float)  # Northing
 
     # Check zone boundaries
@@ -5165 +5176 @@
     if verbose: print 'Converting quantities'
     for j in range(len(times)):
-        for i, index in enumerate(permutation):
-            w = zscale*mux['HA'][index,j] + mean_stage
-            h=w-elevation[index]
+        for i in range(number_of_points):
+            w = zscale*mux['HA'][i,j] + mean_stage
+            h=w-elevation[i]
             stage[j,i] = w
 
-            xmomentum[j,i] = mux['UA'][index,j]*h
-            ymomentum[j,i] = mux['VA'][index,j]*h
+            xmomentum[j,i] = mux['UA'][i,j]*h
+            ymomentum[j,i] = mux['VA'][i,j]*h
 
     outfile.close()

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -6107 +6107 @@
         return base_name, files
 
-    def test_read_mux2_pyI(self):
-        """Constant stage,momentum at each gauge
+    def test_urs2sts_read_mux2_pyI(self):
+        """test_urs2sts_read_mux2_pyI(self):
+        Constant stage,momentum at each gauge
         """
         tide = 1
@@ -6159 +6160 @@
         assert allclose(yvelocity,va)
 
-    def test_read_mux2_pyII(self):
+    def test_urs2sts_read_mux2_pyII(self):
         """Spatially varing stage
         """
@@ -6217 +6218 @@
         assert allclose(yvelocity,va)
 
-    def test_read_mux2_pyIII(self):
+    def test_urs2sts_read_mux2_pyIII(self):
         """Varying start and finsh times
         """
@@ -6587 +6588 @@
         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][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_urs2sts_ordering_exception(self):
+        """Test that inconsistent lats and lons in ordering file are caught.
+        """
+        
+        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 and check that exception is being raised
         _, ordering_filename = tempfile.mkstemp('')
         order_fid = open(ordering_filename, 'w')  
@@ -6613 +6790 @@
 
         
-            
-        # 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):
+        
+        
+    def test_urs2sts_file_boundary_stsI(self):
         from anuga.shallow_water import Domain
         from anuga.shallow_water import Reflective_boundary
@@ -6818 +6867 @@
         os.remove(meshname)
 
-    def test_file_boundary_stsII(self):
+    def test_urs2sts_file_boundary_stsII(self):
         """ mux2 file has points not included in boundary
          mux2 gauges are not stored with the same order as they are 
@@ -6901 +6950 @@
         os.remove(meshname)
 
-    def test_file_boundary_stsIII(self):
+    def test_urs2sts_file_boundary_stsIII_ordering(self):
         """Read correct points from ordering file
         """
@@ -8891 +8940 @@
 if __name__ == "__main__":
 
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts0')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_read_mux2_pyI')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
     #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
@@ -8906 +8955 @@
     else:
         pass
-    runner = unittest.TextTestRunner() # verbosity=2)
+    runner = unittest.TextTestRunner(verbosity=2)
     runner.run(suite)
 

anuga_core/source/anuga/shallow_water/urs_ext.c

@@ -28 +28 @@
 long getNumData(const int *fros, const int *lros, const int nsta);
 char isdata(float);
-float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose);
+float** _read_mux2(int numSrc, 
+                   char **muxFileNameArray, 
+                   float *weights, 
+                   double *params, 
+                   int number_of_selected_stations,
+                   int *permutation,
+                   int verbose);
 
 PyObject *read_mux2(PyObject *self, PyObject *args){
@@ -34 +40 @@
 
     Python call:
-    read_mux2(numSrc,filenames,weights,file_params,verbose)
+    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
 
     NOTE:
@@ -43 +49 @@
     PyArrayObject *pyweights;
     PyArrayObject *file_params;
+    PyArrayObject *permutation;    
     PyArrayObject *pydata;
     PyObject *fname;
@@ -53 +60 @@
     int verbose;
     int nsta0;
+    int number_of_selected_stations;    
     int nt;
     double dt;
@@ -66 +74 @@
 
     // Convert Python arguments to C
-    if (!PyArg_ParseTuple(args, "iOOOi",
-        &numSrc, &filenames, &pyweights, &file_params, &verbose)) 
+    if (!PyArg_ParseTuple(args, "iOOOOi",
+                          &numSrc, &filenames, &pyweights, &file_params, &permutation, &verbose)) 
     {
             PyErr_SetString(PyExc_RuntimeError, 
@@ -106 +114 @@
     }
 
-    //printf("Checkpoint 2 for urs2sts_ext.c\n");    
+    printf("Checkpoint 2 for urs2sts_ext.c\n");    
     //printf("numSrc = %d\n", numSrc);
     for (i = 0; i < numSrc; i++)
@@ -128 +136 @@
         if (muxFileNameArray[i] == NULL) 
         {
-            printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
+            PyErr_SetString(PyExc_ValueError, 
+                            "ERROR: Memory for muxFileNameArray could not be allocated.\n");
             return NULL;
         }
@@ -137 +146 @@
     {
         PyErr_SetString(PyExc_ValueError,
-            "file_params must be one-dimensional and of type double");
+                        "file_params must be one-dimensional and of type double");
         return NULL; 
     }
 
-    //printf("Checkpoint 4 for urs2sts_ext.c\n");        
+    printf("Checkpoint 4 for urs2sts_ext.c\n");        
     // Create array for weights which are passed to read_mux2
     weights = (float*) malloc(numSrc*sizeof(float));
@@ -149 +158 @@
     }
     
-    //printf("Checkpoint 5 for urs2sts_ext.c\n");            
+   
+    number_of_selected_stations = (int) permutation->dimensions[0];
+    printf("Checkpoint 5 for urs2sts_ext.c\n");        
     // Read in mux2 data from file
     cdata = _read_mux2(numSrc, 
                        muxFileNameArray, 
                        weights, 
-                       (double*)file_params->data, 
+                       (double*)file_params->data,
+                       number_of_selected_stations, // Desired number of stations
+                       (int*) permutation->data, // Ordering of selected stations
                        verbose);
 
+    if (!cdata) 
+    {
+        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
+        return NULL;
+    }       
+                       
+                       
     //printf("Checkpoint 6 for urs2sts_ext.c\n");                
     // Allocate space for return vector
@@ -163 +183 @@
     nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
 
-    //printf("Checkpoint 7 for urs2sts_ext.c\n");                    
+    printf("Checkpoint 7 for urs2sts_ext.c\n");                    
+    
     // Find min and max start times of all gauges
     start_tstep = nt + 1;
     finish_tstep = -1;
-    for (i = 0; i < nsta0; i++)
-    {
+    for (i = 0; i < number_of_selected_stations; i++)
+    {
+        printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
+        printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);  
+        
         if ((int)cdata[i][nt + 3] < start_tstep)
         {
@@ -179 +203 @@
     }
 
-    //printf("Checkpoint 8 for urs2sts_ext.c\n");                        
+    printf("Checkpoint 8 for urs2sts_ext.c\n");                        
     if ((start_tstep > nt) | (finish_tstep < 0))
     {
-        printf("ERROR: Gauge data has incorrect start and finsh times\n");
+        printf("ERROR: Gauge data has incorrect start and finish times:\n");
+        printf("   start_tstep = %d, max_number_of_steps = %d\n", start_tstep, nt);
+        printf("   finish_tstep = %d, min_number_of_steps = %d\n", finish_tstep, 0);    
+        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times");  
         return NULL;
     }
@@ -188 +215 @@
     if (start_tstep >= finish_tstep)
     {
-        printf("ERROR: Gauge data has non-postive_length\n");
+        PyErr_SetString(PyExc_ValueError, "ERROR: Gauge data has non-postive_length");
         return NULL;
     }
 
     num_ts = finish_tstep - start_tstep + 1;
-    dimensions[0] = nsta0;
+    dimensions[0] = number_of_selected_stations;
     dimensions[1] = num_ts + POFFSET;
     pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
     if(pydata == NULL)
     {
-        printf("ERROR: Memory for pydata array could not be allocated.\n");
-        return NULL;
-    }
-
-    //printf("Checkpoint 9 for urs2sts_ext.c\n");                            
+        PyErr_SetString(PyExc_ValueError, "ERROR: Memory for pydata array could not be allocated.");
+        return NULL;
+    }
+
+    printf("Checkpoint 9 for urs2sts_ext.c\n");                            
     // Each gauge begins and ends recording at different times. When a gauge is
     // not recording but at least one other gauge is. Pad the non-recording gauge
     // array with zeros.
-    for (i = 0; i < nsta0; i++)
+    for (i = 0; i < number_of_selected_stations; i++)
     {
         time = 0;
@@ -232 +259 @@
     }
 
-    //printf("Checkpoint 10 for urs2sts_ext.c\n");                                
+    printf("Checkpoint 10 for urs2sts_ext.c\n");                                
     free(weights);
     
@@ -239 +266 @@
     free(muxFileNameArray);
     
-    for (i = 0; i < nsta0; ++i)
+    for (i = 0; i < number_of_selected_stations; ++i)
     {
         free(cdata[i]);
@@ -245 +272 @@
     free(cdata);
 
-    //printf("Checkpoint 11 for urs2sts_ext.c\n");                                    
+    printf("Checkpoint 11 for urs2sts_ext.c\n");                                    
     return  PyArray_Return(pydata);
 }
@@ -253 +280 @@
                    float *weights, 
                    double *params, 
+                   int number_of_selected_stations,
+                   int *permutation,
                    int verbose)
 {
@@ -299 +328 @@
     }
 
-    /* loop over all sources, read headers and check compatibility */
+    printf("P1\n");
+    /* Loop over all sources, read headers and check compatibility */
     for (isrc = 0; isrc < numSrc; isrc++)
     {
         muxFileName = muxFileNameArray[isrc];
 
-                /* open the mux file */
+        /* open the mux file */
         if((fp = fopen(muxFileName, "r")) == NULL)
         {
             fprintf(stderr, "cannot open file %s\n", muxFileName);
-            exit(-1);  
+            return NULL;  
         }
         
@@ -329 +359 @@
             if(nsta != nsta0)
             {
-                fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
+                fprintf(stderr,"%s has different number of stations to %s\n", 
+                        muxFileName, 
+                        muxFileNameArray[0]);
                 fclose(fp);
-                exit(-1);   
+                return NULL;   
             }
 
@@ -340 +372 @@
                 if (mytgs[ista].dt != mytgs0[ista].dt)
                 {
-                    fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
+                    fprintf(stderr,"%s has different sampling rate to %s\n", 
+                            muxFileName, 
+                            muxFileNameArray[0]);
                     fclose(fp);
-                    exit(-1);                 
+                    return NULL;                    
                 }   
                 if (mytgs[ista].nt != mytgs0[ista].nt)
                 {
-                    fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
+                    fprintf(stderr,"%s has different series length to %s\n", 
+                            muxFileName, 
+                            muxFileNameArray[0]);
                     fclose(fp);
-                    exit(-1);                 
+                    return NULL;                    
                 }
             }
         }
 
-        /* read the start and stop times for this source */
+        /* Read the start and stop times for this source */
         fread(fros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
         fread(lros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
 
-        /* compute the size of the data block for this source */
+        /* Compute the size of the data block for this source */
         numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
 
-        /* Burbidge: Added a sanity check here */
+        /* Sanity check */
         if (numData < 0)
         {
             fprintf(stderr,"Size of data block appears to be negative!\n");
-            exit(-1);
+            return NULL;            
         }
 
@@ -370 +406 @@
     }
 
+    printf("P2\n");    
     params[0] = (double)nsta0;
     params[1] = (double)mytgs0[0].dt;
     params[2] = (double)mytgs0[0].nt;
 
-    /* make array(s) to hold the demuxed data */
-    //FIXME: This can be reduced to only contain stations given in permutation file
-    sts_data = (float**)malloc(nsta0*sizeof(float *));
+    // Apply rule that an empty permutation file means 'take all stations'
+    // We can change this later by passing in None instead of the empty permutation.
+    
+    printf("number_of_selected_stations B4 = %d\n", number_of_selected_stations);    
+    if (number_of_selected_stations == 0)
+    {
+        number_of_selected_stations = nsta0;  
+      
+        printf("Creating identity permutation vector of length = %d\n", nsta0);
+        // Create the Identity permutation vector
+        permutation = (int *) malloc(number_of_selected_stations*sizeof(int));
+        for (i = 0; i < number_of_selected_stations; i++)
+        {
+          permutation[i] = i;  
+        }
+      
+    }
+    
+    printf("P3\n");        
+    printf("Number of selected stations = %d\n", number_of_selected_stations);
+    /* Make array(s) to hold demuxed data for stations given in the permutation file */
+    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
     if (sts_data == NULL)
     {
         printf("ERROR: Memory for sts_data could not be allocated.\n");
-        exit(-1);
-    }
-
-    len_sts_data = mytgs0[0].nt + POFFSET;
-    for (ista = 0; ista < nsta0; ista++)
-    {
+        return NULL;
+    }
+
+    printf("P4\n");            
+    // For each selected station, allocate space for its data
+    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries (I think)
+    for (i = 0; i < number_of_selected_stations; i++)
+    {
+        printf("P4.1, i=%d\n", i);                
+    
         // Initialise sts_data to zero
-        sts_data[ista] = (float*)calloc(len_sts_data, sizeof(float));
-        if (sts_data[ista] == NULL)
+        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
+        if (sts_data[i] == NULL)
         {
             printf("ERROR: Memory for sts_data could not be allocated.\n");
-            exit(-1);
-        }
-
-        sts_data[ista][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
-        sts_data[ista][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
-        sts_data[ista][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
-        sts_data[ista][mytgs0[0].nt + 3] = (float)fros[ista];
-        sts_data[ista][mytgs0[0].nt + 4] = (float)lros[ista];
-    }
-
+            return NULL;
+        }
+
+        
+        printf("P4.2, i=%d\n", i);                      
+        ista = permutation[i]; // Get global index into mux data
+        
+
+        printf("P4.3, i=%d\n", i);                      
+        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
+        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
+        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
+        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
+        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
+    }
+
+    printf("P5\n");                
     temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
 
@@ -415 +482 @@
             //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
             fprintf(stderr, "cannot open file %s\n", muxFileName);
-            exit(-1);  
+            return NULL;                            
         }
 
         if (verbose){
-            printf("Reading mux file %s\n",muxFileName);
+            printf("Reading mux file %s\n", muxFileName);
         }
 
@@ -430 +497 @@
         fclose(fp);
 
-        /* loop over stations */
-        /* This is where we should only take stations with indices in the permutation array 
-           (and in that order) */
-        
-        //for i in range(len(permution)):
-        //     ista = permutation[i]
+        // loop over stations present in the permutation array 
         //     use ista with mux data
         //     use i with the processed data to be returned         
-        i=0;
         
-        for(ista = 0; ista < nsta0; ista++)
+        for(i = 0; i < number_of_selected_stations; i++)
         {               
+        
+            ista = permutation[i]; // Get global index into mux data    
+            
             /* fill the data0 array from the mux file, and weight it */
-            fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros + isrc*nsta0, lros + isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
+            fillDataArray(ista, 
+                          nsta0, 
+                          mytgs0[ista].nt, 
+                          mytgs0[ista].ig, 
+                          fros + isrc*nsta0, 
+                          lros + isrc*nsta0, 
+                          temp_sts_data, 
+                          &istart, 
+                          &istop, 
+                          muxData);
 
             /* weight appropriately and add */
             for(k = 0; k < mytgs0[ista].nt; k++)
             {
-                if((isdata(sts_data[ista][k])) && isdata(temp_sts_data[k]))
+                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
                 {
-                    sts_data[ista][k] += temp_sts_data[k] * weights[isrc];
+                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
                     //printf("%d,%d,%f\n",ista,k,sts_data[ista][k]);
                 }
                 else
                 {
-                    sts_data[ista][k] = NODATA;
+                    sts_data[i][k] = NODATA;
                 }
             }
@@ -461 +534 @@
     }
 
+    printf("P6\n");                    
     free(muxData);
     free(temp_sts_data);