source: anuga_core/source/anuga/shallow_water/urs_ext.c @ 5418

Last change on this file since 5418 was 5412, checked in by jakeman, 16 years ago

allow urs2sts to accept multiple sources

File size: 19.1 KB
Line 
1/*
2gcc -fPIC -c urs_ext.c -I/usr/include/python2.4 -o urs_ext.o -Wall -O
3gcc -shared urs_ext.o  -o urs_ext.so
4*/
5#include "Python.h"
6#include "Numeric/arrayobject.h"
7#include "math.h"
8#include <stdio.h>
9#include <float.h>
10#include <time.h>
11#include "structure.h"
12
13#define MAX_FILE_NAME_LENGTH 128
14#define NODATA 99.0
15#define EPSILON  0.00001
16
17#define DEBUG 0
18
19#define POFFSET 5 //Number of site_params
20
21void fillDataArray(int, int, int, int, int *, int *, float *, int *, int *, float *);
22long getNumData(int*, int*, int);
23char isdata(float);
24void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int);
25float** _read_mux2(int, char **, float *, double *, int);
26
27PyObject *read_mux2(PyObject *self, PyObject *args){
28/*Read in mux 2 file
29   
30    Python call:
31    read_mux2(numSrc,filenames,weights,file_params,write)
32
33    NOTE:
34    A Python int is equivalent to a C long
35    A Python double corresponds to a C double
36*/
37  PyObject *filenames;
38  PyArrayObject *pyweights,*file_params;
39  PyArrayObject *pydata;
40  PyObject *fname;
41
42  char **muxFileNameArray;
43  float *weights;
44  long numSrc;
45  long write;
46
47  float **cdata;
48  int dimensions[2];
49  int nsta0;
50  int nt;
51  double dt;
52  int i,j;
53  int start_tstep;
54  int finish_tstep;
55  int it,time;
56  int num_ts;
57 
58  // Convert Python arguments to C
59  if (!PyArg_ParseTuple(args, "iOOOi",
60                        &numSrc,&filenames,&pyweights,&file_params,&write)) {                   
61                       
62    PyErr_SetString(PyExc_RuntimeError, 
63                    "Input arguments to read_mux2 failed");
64    return NULL;
65  }
66
67  if(!PyList_Check(filenames)) {
68    PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list");
69    return NULL;
70  }
71
72  if(PyList_Size(filenames) == 0){
73    PyErr_SetString(PyExc_ValueError, "empty lists not allowed");
74    return NULL;
75  }
76
77  if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) {
78    PyErr_SetString(PyExc_ValueError,
79                    "pyweights must be one-dimensional and of type double");
80    return NULL; 
81  }
82
83  if(PyList_Size(filenames) != pyweights->dimensions[0]){
84      PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename");
85      return NULL;
86  }
87
88  muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *));
89  if (muxFileNameArray == NULL) {
90     printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
91     exit(-1);
92  }
93  for (i=0;i<PyList_Size(filenames);i++){
94    muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char));
95    if (muxFileNameArray[i] == NULL) {
96      printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
97      exit(-1);
98    }
99    fname=PyList_GetItem(filenames, i);
100    if (!PyString_Check(fname)) {
101      PyErr_SetString(PyExc_ValueError, "filename not a string");
102    }
103    muxFileNameArray[i]=PyString_AsString(fname);
104  }
105
106  if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) {
107    PyErr_SetString(PyExc_ValueError,
108                    "file_params must be one-dimensional and of type double");
109    return NULL; 
110  }
111
112  weights = (float *) malloc((int)numSrc*sizeof(float));
113  for (i=0;i<(int)numSrc;i++){
114    weights[i]=(float)(*(double *)(pyweights->data+i*pyweights->strides[0]));
115  }
116
117  cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)write);
118
119
120  // Allocate space for return vector
121  nsta0=(int)*(double *) (file_params -> data+0*file_params->strides[0]);
122  dt=*(double *) (file_params -> data+1*file_params->strides[0]);
123  nt=(int)*(double *) (file_params -> data+2*file_params->strides[0]);
124
125  start_tstep=nt+1;
126  finish_tstep=-1;
127  for (i=0;i<nsta0;i++){
128    if ((int)cdata[i][nt+3] < start_tstep){
129      start_tstep=(int)cdata[i][nt+3];
130    }
131    if ((int)cdata[i][nt+4] > finish_tstep){
132      finish_tstep=(int)cdata[i][nt+4]; 
133    }
134  }
135
136  if ((start_tstep>nt) | (finish_tstep < 0)){
137    printf("ERROR: Gauge data has incorrect start and finsh times\n");
138    return NULL;
139  }
140
141  if (start_tstep>=finish_tstep){
142    printf("ERROR: Gauge data has non-postive_length\n");
143    return NULL;
144  }
145
146  num_ts=finish_tstep-start_tstep+1;
147  dimensions[0]=nsta0;
148  dimensions[1]=num_ts+POFFSET;
149  pydata = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
150  if(pydata == NULL){
151    printf("ERROR: Memory for pydata array could not be allocated.\n");
152    return NULL;
153  }
154
155  for (i=0;i<nsta0;i++){
156    time=0;
157    for (it=0;it<finish_tstep;it++){
158      if((it+1>=start_tstep)&&(it+1<=finish_tstep)){
159        if (it+1>(int)cdata[i][nt+4]){
160          //This gauge has stopped recording but others are still recording
161        *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=0.0;
162        }else{
163          *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=cdata[i][it];
164        }
165        time++;
166      }
167    }
168    for (j=0; j<POFFSET; j++){
169      *(double *) (pydata -> data+i*pydata->strides[0]+(num_ts+j)*pydata->strides[1])=cdata[i][nt+j];
170     }
171  }
172
173  free(weights);
174  free(muxFileNameArray);
175  return  PyArray_Return(pydata);
176
177}
178
179float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int write)
180//void _read_mux2(int numSrc, char **muxFileNameArray, float *weights, float **mydata)
181{
182   FILE *fp;
183   int nsta, nsta0, i, isrc, ista;//numSrc;
184   struct tgsrwg *mytgs, *mytgs0;
185   //char *muxFileNameArray;
186   char *muxFileName;                                                                 
187   char outFileName[MAX_FILE_NAME_LENGTH+1];
188   float *wt;
189   float max, amax;
190   float *data, *data0;
191   int istart, istop;
192   int *fros, *lros;
193   char susMuxFileName;
194   float **muxData;
195   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*/     
208   susMuxFileName=0;
209   for(isrc=0;isrc<numSrc;isrc++)
210     { 
211       muxFileName=muxFileNameArray[isrc];
212       wt=weights;
213       if(!susMuxFileName && strcmp(muxFileName+strlen(muxFileName)-4,"mux2")!=0){
214         susMuxFileName=1;
215       }
216     }
217   
218   if(susMuxFileName)
219   {
220      printf("\n**************************************************************************\n");
221      printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
222      printf("   At least one input file name does not end with mux2\n");
223      printf("   Check your results carefully!\n");
224      printf("**************************************************************************\n\n");
225   }   
226                     
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");
233   
234   /* open the first muxfile */
235   if((fp=fopen(muxFileNameArray[0],"r"))==NULL)
236   {
237      fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]);
238      exit(-1); 
239   }
240 
241   /* read in the header */
242   /*first read the number of stations*/   
243   fread(&nsta0,sizeof(int),1,fp);
244   /*now allocate space for, and read in, the structures for each station*/
245   mytgs0 = (struct tgsrwg *) malloc(nsta0*sizeof(struct tgsrwg));
246   fread(&mytgs0[0], nsta0*sizeof(struct tgsrwg), 1, fp);
247
248   /*make an array to hold the start and stop steps for each station for each source*/   
249   fros = (int *) malloc(nsta0*numSrc*sizeof(int));
250   lros = (int *) malloc(nsta0*numSrc*sizeof(int));
251   
252   /* read the start and stop steps for source 0 into the array */   
253   fread(fros,nsta0*sizeof(int),1,fp);
254   fread(lros,nsta0*sizeof(int),1,fp);
255
256   /* compute the size of the data block for source 0 */   
257   numData = getNumData(fros, lros, nsta0);
258
259   /* Burbidge: Added a sanity check here */
260   if (numData < 0) {
261     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);
263     exit(-1);
264   }
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);
269   fclose(fp); 
270
271   if(numSrc > 1){
272      /* allocate space for tgsrwg for the other sources */
273      mytgs = (struct tgsrwg *)malloc( nsta0*sizeof(struct tgsrwg) );
274   } else {
275     /* FIXME (JJ): What should happen in case the are no source files?*/
276     /* If we exit here, tests will break */
277     // fprintf(stderr, "No source file specified\n");
278     // exit(-1);       
279   }
280   
281   /* loop over sources, check compatibility, and read them into *muxData */
282   for(isrc=1; isrc<numSrc; isrc++){
283     muxFileName = muxFileNameArray[isrc];
284     
285     /* open the mux file */
286     if((fp=fopen(muxFileName,"r"))==NULL)
287       {
288         //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
289         fprintf(stderr, "cannot open file %s\n", muxFileName);
290         exit(-1); 
291       }
292     
293     /* check that the mux files are compatible */     
294     fread(&nsta,sizeof(int),1,fp);
295     if(nsta != nsta0){
296       fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
297       fclose(fp);
298       exit(-1);   
299     }
300     fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp);
301     for(ista=0; ista < nsta; ista++){
302       if(mytgs[ista].dt != mytgs0[ista].dt){
303         fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
304         fclose(fp);
305         exit(-1);                 
306       }   
307       if(mytgs[ista].nt != mytgs0[ista].nt){
308         fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
309         fclose(fp);
310         exit(-1);                 
311       }
312     }
313
314      /* read the start and stop times for this source */
315      fread(fros+isrc*nsta0,nsta0*sizeof(int),1,fp);
316      fread(lros+isrc*nsta0,nsta0*sizeof(int),1,fp);
317     
318      /* compute the size of the data block for this source */
319      numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0);
320
321      /* Burbidge: Added a sanity check here */
322      if (numData < 0){
323          fprintf(stderr,"Size of data block appears to be negative!\n");
324          //fprintf(stderr,"numData=%d fros=%d lros=%d nsta0=%d\n");
325          exit(-1);
326      }
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);
331      fclose(fp);             
332   }
333   params[0]=(double)nsta0;
334   params[1]=(double)mytgs0[0].dt;
335   params[2]=(double)mytgs0[0].nt;
336   
337   /* 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");
341     exit(-1);
342   }
343   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");
347       exit(-1);
348     }
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);
369     muxFileName = muxFileNameArray[isrc];
370     //printf("Demuxing station %d\n",ista);
371     //printf("   ... source %s\n",muxFileName);
372     
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 );
433         }
434       }
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;
463}   
464
465
466/* thomas */
467void fillDataArray(int ista, int nsta, int nt, int ig, int *nst, int *nft, float *data, int *istart_p, int *istop_p, float *muxData)
468{
469   int it, last_it, jsta;
470   long int offset=0;
471
472
473   last_it=-1;
474   /* make arrays of starting and finishing time steps for the tide gauges */
475   /* and fill them from the file */
476     
477   /* update start and stop timesteps for this gauge */
478   if(nst[ista]!=-1){
479     if(*istart_p==-1){
480         *istart_p=nst[ista];
481     }else{
482         *istart_p=((nst[ista]<*istart_p)?nst[ista]:*istart_p);
483     }
484   }
485   if(nft[ista]!=-1){
486     if(*istop_p==-1){
487         *istop_p=nft[ista];
488     }else{
489         *istop_p=((nft[ista]<*istop_p)?nft[ista]:*istop_p);
490     }
491   }     
492   if(ig==-1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
493   {
494      /* gauge never started recording, or was outside of all grids, fill array with 0 */
495      for(it=0; it<nt; it++)
496         data[it] = 0.0;
497   }   
498   else
499   {
500      for(it=0; it<nt; it++)
501      {
502         last_it = it;
503         /* skip t record of data block */
504         offset++;
505         /* skip records from earlier tide gauges */
506         for(jsta=0; jsta<ista; jsta++)
507            if(it+1>=nst[jsta]&&it+1<=nft[jsta])
508               offset++;
509               
510         /* deal with the tide gauge at hand */
511         if(it+1>=nst[ista]&&it+1<=nft[ista])
512         /* gauge is recording at this time */
513         {
514            memcpy(data+it,muxData+offset,sizeof(float));
515            offset++;
516         }
517         else if (it+1<nst[ista])
518         {
519            /* gauge has not yet started recording */
520            data[it] = 0.0;
521         }   
522         else
523         /* gauge has finished recording */                                           
524         {
525            data[it] = NODATA;
526            break;
527         }
528   
529         /* skip records from later tide gauges */
530         for(jsta=ista+1; jsta<nsta; jsta++)
531            if(it+1>=nst[jsta]&&it+1<=nft[jsta])
532               offset++;
533      }
534   
535      if(last_it < nt - 1)
536         /* the loop was exited early because the gauge had finished recording */
537         for(it=last_it+1; it < nt; it++)
538            data[it] = NODATA;
539   }
540} 
541
542/* Burbidge: No "offset" is sent. Replace with max. Added grid_id */
543void wrttxt( fname, dt, nt, x, beg, lat, lon, max, depth, grid_id, istart, istop )
544char *fname;
545float dt, *x, beg, max, lat, lon, depth;
546int grid_id;
547int nt;
548int istart, istop;
549{
550   int it;
551   float t;
552   FILE *fp;
553
554   fp = fopen(fname,"w");
555   fprintf(fp,"> lat %g lon %g depth %g max %g start_time %g stop_time %g grid_id %d\n", lat, lon, depth, max, 
556            istart*dt, istop*dt, grid_id );
557   for(it=0;it<nt;it++)
558   {
559      t=beg+it*dt;
560      fprintf(fp,"%9.3f %g\n", t, x[it]);
561   }
562   fclose(fp);
563}
564 
565char isdata(float x)
566{
567  //char value;
568   if(x < NODATA + EPSILON && NODATA < x + EPSILON)
569      return 0;
570   else
571      return 1; 
572}
573
574long getNumData(int *fros, int *lros, int nsta)
575/* calculates the number of data in the data block of a mux file */
576/* based on the first and last recorded output steps for each gauge */ 
577{
578   int ista, last_output_step;
579   long numData = 0;
580
581   last_output_step = 0;   
582   for(ista=0; ista < nsta; ista++)
583      if(*(fros + ista) != -1)
584      {
585         numData += *(lros + ista) - *(fros + ista) + 1;
586         last_output_step = (last_output_step < *(lros+ista) ? *(lros+ista):last_output_step);
587      }   
588   numData += last_output_step*nsta; /* these are the t records */
589   return numData;
590}   
591
592
593
594//-------------------------------
595// Method table for python module
596//-------------------------------
597static struct PyMethodDef MethodTable[] = {
598  {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
599  {NULL, NULL}
600};
601
602// Module initialisation
603void initurs_ext(void){
604  Py_InitModule("urs_ext", MethodTable);
605
606  import_array(); // Necessary for handling of NumPY structures
607}
Note: See TracBrowser for help on using the repository browser.