Changeset 5471


Ignore:
Timestamp:
Jul 7, 2008, 1:23:52 PM (11 years ago)
Author:
ole
Message:

Removed duplication of sts data originally used to pass back to python

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5470 r5471  
    2323char isdata(float);
    2424void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int);
    25 float** _read_mux2(int, char **, float *, double *, int);
     25PyArrayObject *_read_mux2(int, char **, float *, double *, int);
    2626
    2727PyObject *read_mux2(PyObject *self, PyObject *args){
     
    3737  PyObject *filenames;
    3838  PyArrayObject *pyweights,*file_params;
    39   PyArrayObject *pydata;
     39  PyArrayObject *sts_data;
    4040  PyObject *fname;
    4141
     
    4545  long verbose;
    4646
    47   float **cdata;
    48   int dimensions[2];
    4947  int nsta0;
    5048  int nt;
    5149  double dt;
    52   int i,j;
    53   int start_tstep;
    54   int finish_tstep;
    55   int it,time;
    56   int num_ts;
     50  int i;
    5751 
    5852  // Convert Python arguments to C
     
    117111
    118112  // Read in mux2 data from file
    119   cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose);
     113  sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose);
    120114
    121115
     
    125119  nt=(int)*(double *) (file_params -> data+2*file_params->strides[0]);
    126120
    127  
    128   // Find min and max start times of all gauges
    129   start_tstep=nt+1;
    130   finish_tstep=-1;
    131   for (i=0;i<nsta0;i++){
    132     if ((int)cdata[i][nt+3] < start_tstep){
    133       start_tstep=(int)cdata[i][nt+3];
    134     }
    135     if ((int)cdata[i][nt+4] > finish_tstep){
    136       finish_tstep=(int)cdata[i][nt+4];
    137     }
    138   }
    139 
    140   if ((start_tstep>nt) | (finish_tstep < 0)){
    141     printf("ERROR: Gauge data has incorrect start and finsh times\n");
    142     return NULL;
    143   }
    144 
    145   if (start_tstep>=finish_tstep){
    146     printf("ERROR: Gauge data has non-postive_length\n");
    147     return NULL;
    148   }
    149 
    150   num_ts=finish_tstep-start_tstep+1;
    151   dimensions[0]=nsta0;
    152   dimensions[1]=num_ts+POFFSET;
    153   pydata = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    154   if(pydata == NULL){
    155     printf("ERROR: Memory for pydata array could not be allocated.\n");
    156     return NULL;
    157   }
    158 
    159   // Each gauge begins and ends recording at different times. When a gauge is
    160   // not recording but at least one other gauge is. Pad the non-recording gauge
    161   // array with zeros.
    162   for (i=0;i<nsta0;i++){
    163     time=0;
    164     for (it=0;it<finish_tstep;it++){
    165       if((it+1>=start_tstep)&&(it+1<=finish_tstep)){
    166         if (it+1>(int)cdata[i][nt+4]){
    167           //This gauge has stopped recording but others are still recording
    168         *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=0.0;
    169         }else{
    170           *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=cdata[i][it];
    171         }
    172         time++;
    173       }
    174     }
    175     // Pass back lat,lon,elevation
    176     for (j=0; j<POFFSET; j++){
    177       *(double *) (pydata -> data+i*pydata->strides[0]+(num_ts+j)*pydata->strides[1])=cdata[i][nt+j];
    178      }
    179   }
    180 
    181121  free(weights);
    182122  free(muxFileNameArray);
    183   free(cdata);
    184   return  PyArray_Return(pydata);
     123  return  PyArray_Return(sts_data);
    185124
    186125}
    187126
    188 float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
     127PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
    189128{
    190129   FILE *fp;
    191    int nsta, nsta0, i, isrc, ista;
     130   int nsta, nsta0, i, t, isrc, ista;
    192131   struct tgsrwg *mytgs, *mytgs0;
    193132   char *muxFileName;                                                                 
     
    198137   long numData;
    199138
    200 
     139   int start_tstep;
     140   int finish_tstep;
     141   int num_ts;
     142   int dimensions[2];
    201143   int len_sts_data;
    202    float **sts_data;
    203144   float *temp_sts_data;
     145   PyArrayObject *sts_data;
    204146   
    205147   long int offset;
     
    256198   /* compute the size of the data block for source 0 */   
    257199   numData = getNumData(fros, lros, nsta0);
    258 
     200   num_ts = mytgs0[0].nt;
     201   
    259202   /* Burbidge: Added a sanity check here */
    260203   if (numData < 0) {
     
    303246         exit(-1);                 
    304247       }   
    305        if(mytgs[ista].nt != mytgs0[ista].nt){
     248       if(mytgs[ista].nt != num_ts){
    306249         fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
    307250         fclose(fp);
     
    327270   params[1]=(double)mytgs0[0].dt;
    328271   params[2]=(double)mytgs0[0].nt;
    329    
    330    /* make array(s) to hold the demuxed data */
    331    //FIXME: This can be reduced to only contain stations in order file
    332    sts_data = (float **)malloc (nsta0 * sizeof(float *));
    333    
    334    if (sts_data == NULL){
    335      printf("ERROR: Memory for sts_data could not be allocated.\n");
    336      exit(-1);
    337    }
    338    
    339    len_sts_data=mytgs0[0].nt + POFFSET;
    340    for (ista=0; ista<nsta0; ista++){
    341      // Initialise sts_data to zero
    342      sts_data[ista] = (float *)calloc(len_sts_data, sizeof(float) );
    343      if (sts_data[ista] == NULL){
    344        printf("ERROR: Memory for sts_data could not be allocated.\n");
    345        exit(-1);
    346      }
    347      sts_data[ista][mytgs0[0].nt]=(float)mytgs0[ista].geolat;
    348      sts_data[ista][mytgs0[0].nt+1]=(float)mytgs0[ista].geolon;
    349      sts_data[ista][mytgs0[0].nt+2]=(float)mytgs0[ista].z;
    350      sts_data[ista][mytgs0[0].nt+3]=(float)fros[ista];
    351      sts_data[ista][mytgs0[0].nt+4]=(float)lros[ista];
    352    }
     272 
     273   
     274   // Find min and max start times of all gauges
     275   start_tstep=num_ts+1;
     276   finish_tstep=-1;
     277   for (ista=0;ista<nsta0;ista++){
     278     if (fros[ista]< start_tstep){
     279       start_tstep=fros[ista];
     280     }
     281     if (lros[0+nsta0*ista] > finish_tstep){
     282       finish_tstep=lros[ista];
     283     }
     284   }
     285   
     286   if ((start_tstep>num_ts) | (finish_tstep < 0)){
     287     printf("ERROR: Gauge data has incorrect start and finsh times\n");
     288     return NULL;
     289   }
     290   
     291   if (start_tstep>=finish_tstep){
     292     printf("ERROR: Gauge data has non-postive_length\n");
     293     return NULL;
     294   }
     295   
     296   /* Make array(s) to hold the demuxed data */
     297   len_sts_data=num_ts+POFFSET;
     298   dimensions[0]=nsta0;
     299   dimensions[1]=len_sts_data;
     300   sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
     301   if(sts_data == NULL){
     302    printf("ERROR: Memory for pydata array could not be allocated.\n");
     303    return NULL;
     304   }
     305   
     306   /* Initialise sts data to zero */
     307   for (ista=0;ista<nsta0;ista++){
     308     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];
     316   }
    353317
    354318   temp_sts_data = (float *)calloc(len_sts_data, sizeof(float) );
     
    384348       
    385349       /* fill the data0 array from the mux file, and weight it */
    386        fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
     350       fillDataArray(ista, nsta0, num_ts, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
    387351       
    388352       /* weight appropriately and add */
    389        for(i=0; i<mytgs0[ista].nt; i++){
    390          if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i])){
    391            sts_data[ista][i] += temp_sts_data[i] * weights[isrc];
    392            //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]);
     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];
    393356         }else{
    394            sts_data[ista][i] = NODATA;
     357           *(double *) (sts_data -> data+ista*sts_data->strides[0]+i*sts_data->strides[1]) = 0.0;
    395358         }
    396359       }
Note: See TracChangeset for help on using the changeset viewer.