Changeset 5473


Ignore:
Timestamp:
Jul 7, 2008, 3:34:26 PM (11 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:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5470 r5473  
    48384838NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'   
    48394839
    4840 def read_mux2_py(filenames,weights,verbose=False):
    4841 
     4840def read_mux2_py(filenames,weights,permutation=None,verbose=False):
     4841    """
     4842    Access the mux files specified in the filenames list. Combine the
     4843    data found therin as a weighted linear sum as specifed by the weights.
     4844    If permutation is None extract timeseries data for all gauges within the
     4845    files.
     4846   
     4847    Input:
     4848    filenames:   List of filenames specifiying the file containing the
     4849                 timeseries data (mux2 format) for each source
     4850    weights:     Weights associated with each source
     4851    permutation: Specifies the gauge numbers that for which data is to be
     4852                 extracted
     4853    """
    48424854    from Numeric import ones,Float,compress,zeros,arange
    48434855    from urs_ext import read_mux2
     
    48534865        verbose=0
    48544866       
    4855     data=read_mux2(numSrc,filenames,weights,file_params,verbose)
     4867    data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    48564868
    48574869    msg='File parameter values were not read in correctly from c file'
     
    50345046                raise IOError, msg
    50355047
    5036     #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
    5037     #for now set x_momentum and y_moentum quantities to zero
     5048    if ordering_filename is not None:
     5049        # Read ordering file
     5050       
     5051        try:
     5052            fid=open(ordering_filename, 'r')
     5053            file_header=fid.readline().split(',')
     5054            ordering_lines=fid.readlines()
     5055            fid.close()
     5056        except:
     5057            msg = 'Cannot open %s'%ordering_filename
     5058            raise Exception, msg
     5059
     5060        reference_header = 'index, longitude, latitude\n'
     5061        reference_header_split = reference_header.split(',')
     5062        for i in range(3):
     5063            if not file_header[i] == reference_header_split[i]:
     5064                msg = 'File must contain header: '+reference_header+'\n'
     5065                raise Exception, msg
     5066
     5067        if len(ordering_lines)<2:
     5068            msg = 'File must contain at least two points'
     5069            raise Exception, msg
     5070
     5071   
     5072       
     5073        permutation = [int(line.split(',')[0]) for line in ordering_lines]
     5074    else:
     5075        permutation = None 
     5076
    50385077    if (verbose): print 'reading mux2 file'
    50395078    mux={}
    50405079    for i,quantity in enumerate(quantities):
    50415080        # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
    5042         times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose)
     5081        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,permutation=permutation,verbose=verbose)
    50435082        if quantity!=quantities[0]:
    50445083            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
     
    50535092        elevation_old=elevation
    50545093        starttime_old=starttime
    5055 
    50565094       
    50575095    if ordering_filename is not None:
    5058         # Read ordering file
    5059        
    5060         try:
    5061             fid=open(ordering_filename, 'r')
    5062             file_header=fid.readline().split(',')
    5063             ordering_lines=fid.readlines()
    5064             fid.close()
    5065         except:
    5066             msg = 'Cannot open %s'%ordering_filename
    5067             raise Exception, msg
    5068 
    5069         reference_header = 'index, longitude, latitude\n'
    5070         reference_header_split = reference_header.split(',')
    5071         for i in range(3):
    5072             if not file_header[i] == reference_header_split[i]:
    5073                 msg = 'File must contain header: '+reference_header+'\n'
    5074                 raise Exception, msg
    5075 
    5076         if len(ordering_lines)<2:
    5077             msg = 'File must contain at least two points'
    5078             raise Exception, msg
    5079 
    5080    
    5081        
    5082         permutation = [int(line.split(',')[0]) for line in ordering_lines]
    5083    
    5084         latitudes = take(latitudes_urs, permutation)
    5085         longitudes = take(longitudes_urs, permutation)
     5096        latitudes = latitudes_urs
     5097        longitudes = longitudes_urs
    50865098       
    50875099        # Self check - can be removed to improve speed
     
    51165128        longitudes=longitudes_urs
    51175129        permutation = range(latitudes.shape[0])               
    5118        
    51195130
    51205131    msg='File is empty and or clipped region not in file region'
     
    51785189    if verbose: print 'Converting quantities'
    51795190    for j in range(len(times)):
    5180         for i, index in enumerate(permutation):
    5181             w = zscale*mux['HA'][index,j] + mean_stage
    5182             h=w-elevation[index]
     5191        for i in range(number_of_points):
     5192            w = zscale*mux['HA'][i,j] + mean_stage
     5193            h=w-elevation[i]
    51835194            stage[j,i] = w
    51845195
    5185             xmomentum[j,i] = mux['UA'][index,j]*h
    5186             ymomentum[j,i] = mux['VA'][index,j]*h
     5196            xmomentum[j,i] = mux['UA'][i,j]*h
     5197            ymomentum[j,i] = mux['VA'][i,j]*h
    51875198
    51885199    outfile.close()
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5466 r5473  
    61056105        return base_name, files
    61066106
    6107     def test_read_mux2_pyI(self):
     6107    def test_read_mux2_py0(self):
    61086108        """Constant stage,momentum at each gauge
    61096109        """
     
    65436543        """Test multiple sources with ordering file
    65446544        """
    6545        
    65466545        tide=0
    65476546        time_step_count = 3
     
    69456944            fid.write(line)
    69466945        fid.close()
    6947 
     6946       
    69486947        sts_file=base_name
    69496948        urs2sts(base_name, basename_out=sts_file,
     
    88898888if __name__ == "__main__":
    88908889
     8890    #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux2_py0')
    88918891    suite = unittest.makeSuite(Test_Data_Manager,'test')
    88928892    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5471 r5473  
    2323char isdata(float);
    2424void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int);
    25 PyArrayObject *_read_mux2(int, char **, float *, double *, int);
     25PyArrayObject *_read_mux2(int, char **, float *, double *, PyObject *, int);
    2626
    2727PyObject *read_mux2(PyObject *self, PyObject *args){
     
    2929   
    3030    Python call:
    31     read_mux2(numSrc,filenames,weights,file_params,verbose)
     31    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    3232
    3333    NOTE:
     
    3838  PyArrayObject *pyweights,*file_params;
    3939  PyArrayObject *sts_data;
     40  PyObject *permutation;
    4041  PyObject *fname;
    4142
     
    5152 
    5253  // Convert Python arguments to C
    53   if (!PyArg_ParseTuple(args, "iOOOi",
    54                         &numSrc,&filenames,&pyweights,&file_params,&verbose)) {                 
     54  if (!PyArg_ParseTuple(args, "iOOOOi",
     55                        &numSrc,&filenames,&pyweights,&file_params,&permutation,&verbose)) {                   
    5556                       
    5657    PyErr_SetString(PyExc_RuntimeError,
     
    8283  muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *));
    8384  if (muxFileNameArray == NULL) {
    84      printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
    85      exit(-1);
     85     PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
     86     return NULL;
    8687  }
    8788  for (i=0;i<PyList_Size(filenames);i++){
    8889    muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char));
    8990    if (muxFileNameArray[i] == NULL) {
    90       printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
    91       exit(-1);
     91      PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
     92      return NULL;
    9293    }
    9394    fname=PyList_GetItem(filenames, i);
    9495    if (!PyString_Check(fname)) {
    9596      PyErr_SetString(PyExc_ValueError, "filename not a string");
     97      return NULL;
    9698    }
    9799    muxFileNameArray[i]=PyString_AsString(fname);
     
    111113
    112114  // Read in mux2 data from file
    113   sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose);
    114 
     115  sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,permutation,(int)verbose);
    115116
    116117  // Allocate space for return vector
     
    125126}
    126127
    127 PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
     128PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, PyObject *permutation, int verbose)
    128129{
    129130   FILE *fp;
    130    int nsta, nsta0, i, t, isrc, ista;
     131   int nsta, nsta0, i, j, t, isrc, ista, N;
    131132   struct tgsrwg *mytgs, *mytgs0;
    132133   char *muxFileName;                                                                 
     
    144145   float *temp_sts_data;
    145146   PyArrayObject *sts_data;
     147   
     148   int permuation_dimensions[1];
     149   PyArrayObject *permutation_array;
    146150   
    147151   long int offset;
     
    174178   /* Open the first muxfile */
    175179   if((fp=fopen(muxFileNameArray[0],"r"))==NULL){
    176       fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]);
    177       exit(-1); 
     180     fprintf(stderr,"Cannot open file %s\n", muxFileNameArray[0]);
     181     PyErr_SetString(PyExc_RuntimeError,"");
     182     return NULL; 
    178183   }
    179184 
     
    181186   /* First read the number of stations*/   
    182187   fread(&nsta0,sizeof(int),1,fp);
     188   
     189   if (permutation == Py_None){
     190     // if no permutation is specified return the times series of all the gauges in the mux file
     191     permuation_dimensions[0]=nsta0;
     192     permutation_array=(PyArrayObject *) PyArray_FromDims(1, permuation_dimensions, PyArray_INT);
     193     for (ista=0; ista<nsta0; ista++){
     194       *(long *) (permutation_array -> data+ista*permutation_array->strides[0])=ista;
     195     }
     196   }else{
     197     // Specifies the gauge numbers that for which data is to be extracted
     198     permutation_array=(PyArrayObject *) PyArray_ContiguousFromObject(permutation,PyArray_INT,1,1);
     199   }
     200   if(permutation_array == NULL){
     201     PyErr_SetString(PyExc_RuntimeError,"ERROR: Memory for permutation_array array could not be allocated");
     202     return NULL;
     203   }
    183204   
    184205   /*now allocate space for, and read in, the structures for each station*/
     
    203224   if (numData < 0) {
    204225     fprintf(stderr,"Size of data block appears to be negative!\n");
    205      exit(-1);
     226     PyErr_SetString(PyExc_RuntimeError,"");
     227     return NULL; 
    206228   }
    207229   fclose(fp);
     
    229251       {
    230252         fprintf(stderr, "cannot open file %s\n", muxFileName);
    231          exit(-1); 
     253         PyErr_SetString(PyExc_RuntimeError,"");
     254         return NULL;
    232255       }
    233256     
     
    236259     if(nsta != nsta0){
    237260       fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
     261       PyErr_SetString(PyExc_RuntimeError,"");
    238262       fclose(fp);
    239        exit(-1); 
     263       return NULL;
    240264     }
    241265     fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp);
     
    244268         fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
    245269         fclose(fp);
    246          exit(-1);                
     270         return NULL;               
    247271       }   
    248272       if(mytgs[ista].nt != num_ts){
    249273         fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
     274         PyErr_SetString(PyExc_RuntimeError,"");
    250275         fclose(fp);
    251          exit(-1);                
     276         return NULL;             
    252277       }
    253278     }
     
    263288      if (numData < 0){
    264289          fprintf(stderr,"Size of data block appears to be negative!\n");
    265           exit(-1);
     290          PyErr_SetString(PyExc_RuntimeError,"");
     291          return NULL;
    266292      }
    267293      fclose(fp);             
    268294   }
    269    params[0]=(double)nsta0;
     295   N = permutation_array->dimensions[0];
     296   
     297   params[0]=(double)N;
    270298   params[1]=(double)mytgs0[0].dt;
    271    params[2]=(double)mytgs0[0].nt;
     299   params[2]=(double)num_ts;
    272300 
    273301   
    274302   // Find min and max start times of all gauges
     303   // FIXME: ONLY CHECK gauges in permutation for start and finsh times
    275304   start_tstep=num_ts+1;
    276305   finish_tstep=-1;
     
    285314   
    286315   if ((start_tstep>num_ts) | (finish_tstep < 0)){
    287      printf("ERROR: Gauge data has incorrect start and finsh times\n");
     316     PyErr_SetString(PyExc_RuntimeError,"Gauge data has incorrect start and finsh times");
    288317     return NULL;
    289318   }
    290319   
    291320   if (start_tstep>=finish_tstep){
    292      printf("ERROR: Gauge data has non-postive_length\n");
     321     PyErr_SetString(PyExc_RuntimeError,"Gauge data has non-postive_length");
    293322     return NULL;
    294323   }
     
    296325   /* Make array(s) to hold the demuxed data */
    297326   len_sts_data=num_ts+POFFSET;
    298    dimensions[0]=nsta0;
     327   dimensions[0]=N;
    299328   dimensions[1]=len_sts_data;
    300329   sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    301330   if(sts_data == NULL){
    302     printf("ERROR: Memory for pydata array could not be allocated.\n");
     331    PyErr_SetString(PyExc_RuntimeError,"Memory for pydata array could not be allocated");
    303332    return NULL;
    304333   }
    305334   
    306335   /* Initialise sts data to zero */
    307    for (ista=0;ista<nsta0;ista++){
     336   for(i=0; i<N; i++){
     337     ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]);
    308338     for (t=0;t<num_ts;t++){
    309        *(double *) (sts_data -> data+ista*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
    310      }
    311      *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
    312      *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
    313      *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
    314      *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
    315      *(double *) (sts_data -> data+ista*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
     339       *(double *) (sts_data -> data+i*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
     340     }
     341     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
     342     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
     343     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
     344     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
     345     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
    316346   }
    317347
     
    327357     if((fp=fopen(muxFileName,"r"))==NULL)
    328358       {
    329          //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
    330          fprintf(stderr, "cannot open file %s\n", muxFileName);
    331          exit(-1); 
     359         fprintf(stderr, "Cannot open file %s\n", muxFileName);
     360         PyErr_SetString(PyExc_RuntimeError,"");
     361         return NULL;
    332362       }
    333363       
     
    345375
    346376     /* loop over stations */
    347      for(ista=0; ista<nsta0; ista++){               
     377     for(j=0; j<N; j++){
     378       ista = *(long *) (permutation_array -> data+j*permutation_array->strides[0]);
     379       //printf("%d\n",(int)*(long *) (permutation_array -> data+j*permutation_array->strides[0]));
    348380       
    349381       /* fill the data0 array from the mux file, and weight it */
     
    351383       
    352384       /* weight appropriately and add */
    353        for(i=0; i<num_ts; i++){
    354          if((isdata(*(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]))) && isdata(temp_sts_data[i])){
    355            *(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]) += temp_sts_data[i] * weights[isrc];
     385       for(t=0; t<num_ts; t++){
     386         if((isdata(*(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]))) && isdata(temp_sts_data[t])){
     387           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) += temp_sts_data[t] * weights[isrc];
    356388         }else{
    357            *(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]) = 0.0;
     389           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) = 0.0;
    358390         }
    359391       }
    360392     }
    361    }
     393   }   
    362394   free(muxData);
    363395   free(temp_sts_data);
     396   free(fros);
     397   free(lros);
    364398   return sts_data;
    365399}   
Note: See TracChangeset for help on using the changeset viewer.