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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.