Changeset 5469


Ignore:
Timestamp:
Jul 7, 2008, 11:52:19 AM (16 years ago)
Author:
ole
Message:

John and Ole updated urs_ext.c to allow multiple source files to be read in sequence using only the memory to read one file

File:
1 edited

Legend:

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

    r5412 r5469  
    11/*
    2 gcc -fPIC -c urs_ext.c -I/usr/include/python2.4 -o urs_ext.o -Wall -O
     2gcc -fPIC -c urs_ext.c -I/usr/include/python2.5 -o urs_ext.o -Wall -O
    33gcc -shared urs_ext.o  -o urs_ext.so
    44*/
     
    110110  }
    111111
     112  // Create array for weights which are passed to read_mux2
    112113  weights = (float *) malloc((int)numSrc*sizeof(float));
    113114  for (i=0;i<(int)numSrc;i++){
     
    115116  }
    116117
     118  // Read in mux2 data from file
    117119  cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)write);
    118120
     
    123125  nt=(int)*(double *) (file_params -> data+2*file_params->strides[0]);
    124126
     127 
     128  // Find min and max start times of all gauges
    125129  start_tstep=nt+1;
    126130  finish_tstep=-1;
     
    153157  }
    154158
     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.
    155162  for (i=0;i<nsta0;i++){
    156163    time=0;
     
    166173      }
    167174    }
     175    // Pass back lat,lon,elevation
    168176    for (j=0; j<POFFSET; j++){
    169177      *(double *) (pydata -> data+i*pydata->strides[0]+(num_ts+j)*pydata->strides[1])=cdata[i][nt+j];
     
    178186
    179187float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int write)
    180 //void _read_mux2(int numSrc, char **muxFileNameArray, float *weights, float **mydata)
    181188{
    182189   FILE *fp;
    183    int nsta, nsta0, i, isrc, ista;//numSrc;
     190   int nsta, nsta0, i, isrc, ista;
    184191   struct tgsrwg *mytgs, *mytgs0;
    185    //char *muxFileNameArray;
    186192   char *muxFileName;                                                                 
    187    char outFileName[MAX_FILE_NAME_LENGTH+1];
    188    float *wt;
    189    float max, amax;
    190    float *data, *data0;
    191193   int istart, istop;
    192194   int *fros, *lros;
    193195   char susMuxFileName;
    194    float **muxData;
     196   float *muxData;
    195197   long numData;
    196    time_t   start_time, stop_time;
    197 
    198    float **mydata;
    199    
    200    /* note starting time */
    201    time(&start_time);
    202 
    203    /* allocate space for the names and the weights and pointers to the data*/   
    204    wt = (float *) malloc(numSrc*sizeof(float));
    205    muxData = (float**) malloc(numSrc*sizeof(float*));
    206    
    207    /*read the mux file names and the associated weight from stdin*/     
     198
     199
     200   int len_sts_data;
     201   float **sts_data;
     202   float *temp_sts_data;
     203   
     204   long int offset;
     205
     206   /* Allocate space for the names and the weights and pointers to the data*/   
     207   
     208   /* Check that the input files have mux2 extension*/
    208209   susMuxFileName=0;
    209210   for(isrc=0;isrc<numSrc;isrc++)
    210211     {
    211212       muxFileName=muxFileNameArray[isrc];
    212        wt=weights;
    213213       if(!susMuxFileName && strcmp(muxFileName+strlen(muxFileName)-4,"mux2")!=0){
    214214         susMuxFileName=1;
     
    225225   }   
    226226                     
    227    //printf("Demuxing mux files with weights:\n\n");
    228    for(isrc=0;isrc<numSrc;isrc++){
    229      muxFileName = muxFileNameArray[isrc];
    230      //printf("%s %f\n", muxFileName, *(wt+isrc));
    231    }
    232    //printf("\n");
    233227   
    234228   /* open the first muxfile */
     
    242236   /*first read the number of stations*/   
    243237   fread(&nsta0,sizeof(int),1,fp);
     238   
    244239   /*now allocate space for, and read in, the structures for each station*/
    245240   mytgs0 = (struct tgsrwg *) malloc(nsta0*sizeof(struct tgsrwg));
    246241   fread(&mytgs0[0], nsta0*sizeof(struct tgsrwg), 1, fp);
    247242
    248    /*make an array to hold the start and stop steps for each station for each source*/   
     243   /*make an array to hold the start and stop steps for each station for each
     244     source*/   
     245   //FIXME: Should only store start and stop times for one source
    249246   fros = (int *) malloc(nsta0*numSrc*sizeof(int));
    250247   lros = (int *) malloc(nsta0*numSrc*sizeof(int));
     
    260257   if (numData < 0) {
    261258     fprintf(stderr,"Size of data block appears to be negative!\n");
    262      //fprintf(stderr,"numData=%d fros=%d lros=%d nsta0=%d\n",numData,fros,lros,nsta0);
    263259     exit(-1);
    264260   }
    265    
    266    /* allocate space for these data, read them and close the file */   
    267    *muxData = (float*) malloc(numData*sizeof(float));
    268    fread(*muxData, numData*sizeof(float),1,fp);
    269261   fclose(fp);
    270262
     263   // Allocate header array for each remaining source file.
     264   // FIXME: only need to store for one source and which is compared to all subsequent
     265   // source headers
    271266   if(numSrc > 1){
    272267      /* allocate space for tgsrwg for the other sources */
     
    279274   }
    280275   
    281    /* loop over sources, check compatibility, and read them into *muxData */
     276   /* loop over remaining sources, check compatibility, and read them into
     277    *muxData */
    282278   for(isrc=1; isrc<numSrc; isrc++){
    283279     muxFileName = muxFileNameArray[isrc];
     
    286282     if((fp=fopen(muxFileName,"r"))==NULL)
    287283       {
    288          //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
    289284         fprintf(stderr, "cannot open file %s\n", muxFileName);
    290285         exit(-1); 
     
    322317      if (numData < 0){
    323318          fprintf(stderr,"Size of data block appears to be negative!\n");
    324           //fprintf(stderr,"numData=%d fros=%d lros=%d nsta0=%d\n");
    325319          exit(-1);
    326320      }
    327 
    328       /* allocate space, read the data and close the file */
    329       *(muxData+isrc) = (float*) malloc(numData*sizeof(float));
    330       fread(*(muxData+isrc), numData*sizeof(float),1,fp);
    331321      fclose(fp);             
    332322   }
     
    336326   
    337327   /* make array(s) to hold the demuxed data */
    338    mydata = (float **)malloc (nsta0 * sizeof(float *));
    339    if (mydata == NULL){
    340      printf("ERROR: Memory for mydata could not be allocated.\n");
     328   //FIXME: This can be reduced to only contain stations in order file
     329   sts_data = (float **)malloc (nsta0 * sizeof(float *));
     330   
     331   if (sts_data == NULL){
     332     printf("ERROR: Memory for sts_data could not be allocated.\n");
    341333     exit(-1);
    342334   }
     335   
     336   len_sts_data=mytgs0[0].nt + POFFSET;
    343337   for (ista=0; ista<nsta0; ista++){
    344      mydata[ista] = (float *)malloc( (mytgs0[0].nt + POFFSET)* sizeof(float) );
    345      if (mydata[ista] == NULL){
    346        printf("ERROR: Memory for mydata could not be allocated.\n");
     338     // Initialise sts_data to zero
     339     sts_data[ista] = (float *)calloc(len_sts_data, sizeof(float) );
     340     if (sts_data[ista] == NULL){
     341       printf("ERROR: Memory for sts_data could not be allocated.\n");
    347342       exit(-1);
    348343     }
    349      mydata[ista][mytgs0[0].nt]=(float)mytgs0[ista].geolat;
    350      mydata[ista][mytgs0[0].nt+1]=(float)mytgs0[ista].geolon;
    351      mydata[ista][mytgs0[0].nt+2]=(float)mytgs0[ista].z;
    352      mydata[ista][mytgs0[0].nt+3]=(float)fros[ista];
    353      mydata[ista][mytgs0[0].nt+4]=(float)lros[ista];
    354    }
    355 
    356    //data0 = (float *)malloc( mytgs0[0].nt * sizeof(float) );
    357    if(numSrc > 1)
    358       data = (float *)malloc( mytgs0[0].nt * sizeof(float) );
    359          
    360    /* loop over stations */
    361    for(ista=0; ista<nsta0; ista++){             
    362      data0= mydata[ista];
    363      data=data0;
    364      istart = -1;
    365      istop = -1;
    366      /* fill the data0 array from the first mux file, and weight it */
    367      isrc=0;
    368      //muxFileName = muxFileNameArray + isrc*(MAX_FILE_NAME_LENGTH+1);
     344     sts_data[ista][mytgs0[0].nt]=(float)mytgs0[ista].geolat;
     345     sts_data[ista][mytgs0[0].nt+1]=(float)mytgs0[ista].geolon;
     346     sts_data[ista][mytgs0[0].nt+2]=(float)mytgs0[ista].z;
     347     sts_data[ista][mytgs0[0].nt+3]=(float)fros[ista];
     348     sts_data[ista][mytgs0[0].nt+4]=(float)lros[ista];
     349   }
     350
     351   temp_sts_data = (float *)calloc(len_sts_data, sizeof(float) );
     352     
     353   /* Loop over all sources */
     354   //FIXME: remove istart and istop they are not used.
     355   istart = -1;
     356   istop = -1;
     357   for (isrc=0;isrc<numSrc;isrc++){
     358     /* Read in data block from mux2 file */
    369359     muxFileName = muxFileNameArray[isrc];
    370      //printf("Demuxing station %d\n",ista);
    371      //printf("   ... source %s\n",muxFileName);
     360     if((fp=fopen(muxFileName,"r"))==NULL)
     361       {
     362         //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
     363         fprintf(stderr, "cannot open file %s\n", muxFileName);
     364         exit(-1); 
     365       }
     366
     367     offset=sizeof(int)+nsta0*(sizeof(struct tgsrwg)+2*sizeof(int));
     368     fseek(fp,offset,0);
    372369     
    373      fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, data0, &istart, &istop, *(muxData+isrc));
    374      
    375      /* apply the weights */
    376      for(i=0; i<mytgs0[ista].nt; i++)
    377        if(isdata(*(data0+i)))
    378          *(data0+i) *= *wt;
    379      
    380      /* loop over the rest of the sources */
    381      for(isrc=1;isrc<numSrc;isrc++)
    382        {
    383          /* fill the data array */
    384          //muxFileName = muxFileNameArray + isrc*(MAX_FILE_NAME_LENGTH+1);
    385          muxFileName = muxFileNameArray[isrc];
    386          //printf("   ... source %s\n",muxFileName);
    387          fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, data, &istart, &istop, *(muxData+isrc));
    388          
    389          /* weight appropriately and add */
    390          for(i=0; i<mytgs0[ista].nt; i++)
    391            {
    392              if(isdata(*(data0+i))&&isdata(*(data+i)))
    393                *(data0+i) += *(data+i)* *(wt+isrc);                                                                   
    394              else
    395                *(data0+i) = NODATA;
    396            }     
    397          
    398        }  /* end of loop over sources */
    399      
    400      /* now compute the maxima for this station */     
    401      max = 0.0;
    402      amax = 0.0;
    403      for(i=0; i<mytgs0[ista].nt; i++){
    404        if(isdata(*(data0+i)))
    405          {
    406            max = ((*(data0+i) > max) ? *(data0+i):max);
    407            amax = ((fabs(*(data0+i)) > amax) ? fabs(*(data0+i)):amax); 
    408          }   
    409      }
    410      /* write out sac file for the current station  */
    411      /*thomas - instead of passing beg(=0), should pass dt*/
    412      /*thomas - uncomment the following if you want sac files*/
    413      /*sprintf(outFileName,"S%03d.sac", ista );
    414        printf("   ... writing file=%s\n", outFileName);
    415        wrtsac2(outFileName, mytgs0[ista].dt, mytgs0[ista].nt, &data0[0], mytgs0[ista].dt,
    416        mytgs0[ista].geolat, mytgs0[ista].geolon,
    417        mytgs0[ista].center_lat, mytgs0[ista].center_lon,
    418        mytgs0[ista].offset, mytgs0[ista].z );*/
    419      
    420      /* write out text file for the current station  */
    421      /* DB added check for non-zero max */
    422      if (max >0)
    423        {
    424          /* Burbidge: added tide gauge grid id output. no .id field in tgsrwg */
    425          /* thomas: instead of passing beg(=0), pass dt */   
    426          /*      wrttxt(outFileName, mytgs0[ista].dt, mytgs0[ista].nt, &data0[0], beg, mytgs0[ista].geolat,
    427                  mytgs0[ista].geolon, max, mytgs0[ista].z, mytgs0[ista].ig ); */
    428          if (write) {
    429            sprintf(outFileName,"S%5.5d.txt", ista );
    430            //printf("   ... writing file=%s\n", outFileName);
    431            wrttxt(outFileName, mytgs0[ista].dt, mytgs0[ista].nt, &data0[0], mytgs0[ista].dt, mytgs0[ista].geolat,
    432                   mytgs0[ista].geolon, max, mytgs0[ista].z, mytgs0[ista].ig, istart, istop );
     370     numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0);
     371     muxData = (float*) malloc(numData*sizeof(float));
     372     fread(muxData, numData*sizeof(float),1,fp);
     373     fclose(fp);         
     374
     375     /* loop over stations */
     376     for(ista=0; ista<nsta0; ista++){               
     377       
     378       /* fill the data0 array from the mux file, and weight it */
     379       fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
     380       
     381       /* weight appropriately and add */
     382       for(i=0; i<mytgs0[ista].nt; i++){
     383         if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i])){
     384           sts_data[ista][i] += temp_sts_data[i] * weights[isrc];
     385           //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]);
     386         }else{
     387           sts_data[ista][i] = NODATA;
    433388         }
    434389       }
    435      
    436    }  /* end of loop over stations */
    437    //free(data0);
    438    //free(mytgs0);
    439    
    440    //if(numSrc>1)
    441    //{
    442      //free(data);
    443      //free(mytgs);
    444    //}
    445    //can't free arrays because I only fill array by making pointer not copy
    446    //for(isrc=0; isrc<numSrc;isrc++)
    447    //   free(*(muxData+isrc));
    448    //   
    449    //free(muxData);
    450    // free(muxFileNameArray);
    451    
    452    if(susMuxFileName)
    453    {
    454       printf("\n**************************************************************************\n");
    455       printf("   WARNING: This program operates only on multiplexed files in mux2 format\n");
    456       printf("   At least one input filename does not end with mux2\n");
    457       printf("   Check your results carefully!\n");
    458       printf("**************************************************************************\n");
    459    }   
    460    time(&stop_time);
    461    //fprintf(stdout, "\nElapsed time: %10.0f seconds\n", difftime(stop_time, start_time));
    462    return mydata;
     390     }
     391   }
     392   free(muxData);
     393   free(temp_sts_data);
     394   return sts_data;
    463395}   
    464396
Note: See TracChangeset for help on using the changeset viewer.