Changeset 5473

Timestamp:

Jul 7, 2008, 3:34:26 PM (15 years ago)

Author:

ole

Message:

Allow urs_ext to store only a subset of gauges in mux file according to a list of ordered points

Location:

anuga_core/source/anuga/shallow_water

Files:

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

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4838 +4838 @@
 NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'    
 
-def read_mux2_py(filenames,weights,verbose=False):
-
+def read_mux2_py(filenames,weights,permutation=None,verbose=False):
+    """
+    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
+    files.
+    
+    Input:
+    filenames:   List of filenames specifiying the file containing the
+                 timeseries data (mux2 format) for each source
+    weights:     Weights associated with each source
+    permutation: Specifies the gauge numbers that for which data is to be
+                 extracted
+    """
     from Numeric import ones,Float,compress,zeros,arange
     from urs_ext import read_mux2
@@ -4853 +4865 @@
         verbose=0
         
-    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'
@@ -5034 +5046 @@
                 raise IOError, msg
 
-    #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 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'
+        reference_header_split = reference_header.split(',')
+        for i in range(3):
+            if not file_header[i] == reference_header_split[i]:
+                msg = 'File must contain header: '+reference_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]
+    else:
+        permutation = None  
+
     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 extenstion associated with that quantity
-        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose)
+        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,permutation=permutation,verbose=verbose)
         if quantity!=quantities[0]:
             msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
@@ -5053 +5092 @@
         elevation_old=elevation
         starttime_old=starttime
-
         
     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'
-        reference_header_split = reference_header.split(',')
-        for i in range(3):
-            if not file_header[i] == reference_header_split[i]:
-                msg = 'File must contain header: '+reference_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)
+        latitudes = latitudes_urs
+        longitudes = longitudes_urs
         
         # Self check - can be removed to improve speed
@@ -5116 +5128 @@
         longitudes=longitudes_urs
         permutation = range(latitudes.shape[0])                
-        
 
     msg='File is empty and or clipped region not in file region'
@@ -5178 +5189 @@
     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

@@ -6105 +6105 @@
         return base_name, files
 
-    def test_read_mux2_pyI(self):
+    def test_read_mux2_py0(self):
         """Constant stage,momentum at each gauge
         """
@@ -6543 +6543 @@
         """Test multiple sources with ordering file
         """
-        
         tide=0
         time_step_count = 3
@@ -6945 +6944 @@
             fid.write(line)
         fid.close()
-
+        
         sts_file=base_name
         urs2sts(base_name, basename_out=sts_file,
@@ -8889 +8888 @@
 if __name__ == "__main__":
 
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux2_py0')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')

anuga_core/source/anuga/shallow_water/urs_ext.c

@@ -23 +23 @@
 char isdata(float);
 void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int);
-PyArrayObject *_read_mux2(int, char **, float *, double *, int);
+PyArrayObject *_read_mux2(int, char **, float *, double *, PyObject *, int);
 
 PyObject *read_mux2(PyObject *self, PyObject *args){
@@ -29 +29 @@
    
     Python call:
-    read_mux2(numSrc,filenames,weights,file_params,verbose)
+    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
 
     NOTE:
@@ -38 +38 @@
   PyArrayObject *pyweights,*file_params;
   PyArrayObject *sts_data;
+  PyObject *permutation;
   PyObject *fname;
 
@@ -51 +52 @@
   
   // 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, 
@@ -82 +83 @@
   muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *));
   if (muxFileNameArray == NULL) {
-     printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
-     exit(-1);
+     PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
+     return NULL;
   }
   for (i=0;i<PyList_Size(filenames);i++){
     muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char));
     if (muxFileNameArray[i] == NULL) {
-      printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
-      exit(-1);
+      PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
+      return NULL;
     }
     fname=PyList_GetItem(filenames, i);
     if (!PyString_Check(fname)) {
       PyErr_SetString(PyExc_ValueError, "filename not a string");
+      return NULL;
     }
     muxFileNameArray[i]=PyString_AsString(fname);
@@ -111 +113 @@
 
   // Read in mux2 data from file
-  sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose);
-
+  sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,permutation,(int)verbose);
 
   // Allocate space for return vector
@@ -125 +126 @@
 }
 
-PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
+PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, PyObject *permutation, int verbose)
 {
    FILE *fp;
-   int nsta, nsta0, i, t, isrc, ista;
+   int nsta, nsta0, i, j, t, isrc, ista, N;
    struct tgsrwg *mytgs, *mytgs0;
    char *muxFileName;                                                                  
@@ -144 +145 @@
    float *temp_sts_data;
    PyArrayObject *sts_data;
+   
+   int permuation_dimensions[1];
+   PyArrayObject *permutation_array;
    
    long int offset;
@@ -174 +178 @@
    /* Open the first muxfile */
    if((fp=fopen(muxFileNameArray[0],"r"))==NULL){
-      fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]);
-      exit(-1);  
+     fprintf(stderr,"Cannot open file %s\n", muxFileNameArray[0]);
+     PyErr_SetString(PyExc_RuntimeError,"");
+     return NULL;  
    }
  
@@ -181 +186 @@
    /* First read the number of stations*/   
    fread(&nsta0,sizeof(int),1,fp);
+   
+   if (permutation == Py_None){
+     // if no permutation is specified return the times series of all the gauges in the mux file
+     permuation_dimensions[0]=nsta0;
+     permutation_array=(PyArrayObject *) PyArray_FromDims(1, permuation_dimensions, PyArray_INT);
+     for (ista=0; ista<nsta0; ista++){
+       *(long *) (permutation_array -> data+ista*permutation_array->strides[0])=ista;
+     } 
+   }else{
+     // Specifies the gauge numbers that for which data is to be extracted
+     permutation_array=(PyArrayObject *) PyArray_ContiguousFromObject(permutation,PyArray_INT,1,1);
+   }
+   if(permutation_array == NULL){
+     PyErr_SetString(PyExc_RuntimeError,"ERROR: Memory for permutation_array array could not be allocated");
+     return NULL;
+   }
    
    /*now allocate space for, and read in, the structures for each station*/
@@ -203 +224 @@
    if (numData < 0) {
      fprintf(stderr,"Size of data block appears to be negative!\n");
-     exit(-1);
+     PyErr_SetString(PyExc_RuntimeError,"");
+     return NULL;  
    }
    fclose(fp); 
@@ -229 +251 @@
        {
          fprintf(stderr, "cannot open file %s\n", muxFileName);
-         exit(-1);  
+         PyErr_SetString(PyExc_RuntimeError,"");
+         return NULL; 
        }
      
@@ -236 +259 @@
      if(nsta != nsta0){
        fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
+       PyErr_SetString(PyExc_RuntimeError,"");
        fclose(fp);
-       exit(-1);   
+       return NULL; 
      }
      fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp);
@@ -244 +268 @@
          fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
          fclose(fp);
-         exit(-1);                 
+         return NULL;                
        }   
        if(mytgs[ista].nt != num_ts){
          fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
+         PyErr_SetString(PyExc_RuntimeError,"");
          fclose(fp);
-         exit(-1);                 
+         return NULL;              
        }
      }
@@ -263 +288 @@
       if (numData < 0){
           fprintf(stderr,"Size of data block appears to be negative!\n");
-          exit(-1);
+          PyErr_SetString(PyExc_RuntimeError,"");
+          return NULL;
       }
       fclose(fp);             
    }
-   params[0]=(double)nsta0;
+   N = permutation_array->dimensions[0];
+   
+   params[0]=(double)N;
    params[1]=(double)mytgs0[0].dt;
-   params[2]=(double)mytgs0[0].nt;
+   params[2]=(double)num_ts;
   
     
    // Find min and max start times of all gauges
+   // FIXME: ONLY CHECK gauges in permutation for start and finsh times
    start_tstep=num_ts+1;
    finish_tstep=-1;
@@ -285 +314 @@
    
    if ((start_tstep>num_ts) | (finish_tstep < 0)){
-     printf("ERROR: Gauge data has incorrect start and finsh times\n");
+     PyErr_SetString(PyExc_RuntimeError,"Gauge data has incorrect start and finsh times");
      return NULL;
    }
    
    if (start_tstep>=finish_tstep){
-     printf("ERROR: Gauge data has non-postive_length\n");
+     PyErr_SetString(PyExc_RuntimeError,"Gauge data has non-postive_length");
      return NULL;
    }
@@ -296 +325 @@
    /* Make array(s) to hold the demuxed data */
    len_sts_data=num_ts+POFFSET;
-   dimensions[0]=nsta0;
+   dimensions[0]=N;
    dimensions[1]=len_sts_data;
    sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    if(sts_data == NULL){
-    printf("ERROR: Memory for pydata array could not be allocated.\n");
+    PyErr_SetString(PyExc_RuntimeError,"Memory for pydata array could not be allocated");
     return NULL;
    }
    
    /* Initialise sts data to zero */
-   for (ista=0;ista<nsta0;ista++){
+   for(i=0; i<N; i++){ 
+     ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]);
      for (t=0;t<num_ts;t++){
-       *(double *) (sts_data -> data+ista*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
-     }
-     *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
-     *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
-     *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
-     *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
-     *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
+       *(double *) (sts_data -> data+i*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
+     }
+     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
+     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
+     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
+     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
+     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
    } 
 
@@ -327 +357 @@
      if((fp=fopen(muxFileName,"r"))==NULL)
        {
-         //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
-         fprintf(stderr, "cannot open file %s\n", muxFileName);
-         exit(-1);  
+         fprintf(stderr, "Cannot open file %s\n", muxFileName);
+         PyErr_SetString(PyExc_RuntimeError,"");
+         return NULL;
        }
        
@@ -345 +375 @@
 
      /* loop over stations */
-     for(ista=0; ista<nsta0; ista++){               
+     for(j=0; j<N; j++){ 
+       ista = *(long *) (permutation_array -> data+j*permutation_array->strides[0]);
+       //printf("%d\n",(int)*(long *) (permutation_array -> data+j*permutation_array->strides[0]));
        
        /* fill the data0 array from the mux file, and weight it */
@@ -351 +383 @@
        
        /* weight appropriately and add */
-       for(i=0; i<num_ts; i++){
-         if((isdata(*(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]))) && isdata(temp_sts_data[i])){
-           *(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]) += temp_sts_data[i] * weights[isrc];
+       for(t=0; t<num_ts; t++){
+         if((isdata(*(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]))) && isdata(temp_sts_data[t])){
+           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) += temp_sts_data[t] * weights[isrc];
          }else{
-           *(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]) = 0.0;
+           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) = 0.0;
          }
        }
      }
-   }
+   }   
    free(muxData);
    free(temp_sts_data);
+   free(fros);
+   free(lros);
    return sts_data;
 }