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

Last change on this file since 6873 was 6752, checked in by duncan, 16 years ago

fix for ticket #304. fopen(muxFileName, "r") changed to fopen(muxFileName, "rb"). Hard to find. Easy to implement.

File size: 25.4 KB
Line 
1/*
2gcc -fPIC -c urs_ext.c -I/usr/include/python2.5 -o urs_ext.o -Wall -O
3gcc -shared urs_ext.o  -o urs_ext.so
4*/
5
6#include "Python.h"
7#include "Numeric/arrayobject.h"
8#include "structure.h"
9#include "math.h"
10#include <stdio.h>
11#include <string.h>
12#include <errno.h>
13#include <float.h>
14#include <time.h>
15
16#define MAX_FILE_NAME_LENGTH 128
17#define NODATA 99.0
18#define EPSILON  0.00001
19
20#define DEBUG 0
21
22#define POFFSET 5 //Number of site_params
23
24static int *fros=NULL;  // First recorded output step
25static int *lros=NULL;  // Last recorded output step
26static struct tgsrwg* mytgs0=NULL;
27
28static long numDataMax=0;
29
30
31/*The MUX file format
32
33
34 */
35
36
37
38/////////////////////////////////////////////////////////////////////////
39//Auxiliary functions
40void fillDataArray(int ista, int total_number_of_stations, int nt, int ig, int *nst, 
41                   int *nft, float *data, int *istart_p, 
42                   int *istop_p, float *muxData)
43{
44    int it, last_it, jsta;
45    long int offset=0;
46
47
48    last_it = -1;
49    /* Make arrays of starting and finishing time steps for the tide gauges */
50    /* and fill them from the file */
51
52    /* Update start and stop timesteps for this gauge */
53    if (nst[ista]!= -1)
54    {
55        if(*istart_p == -1)
56        {
57            *istart_p = nst[ista];
58        }
59        else
60        {
61            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
62        }
63    }
64   
65    if (nft[ista] != -1)
66    {
67        if (*istop_p == -1)
68        {
69            *istop_p = nft[ista];
70        }
71        else
72        {
73            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
74        }
75    }     
76   
77    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
78    {
79        /* gauge never started recording, or was outside of all grids,
80        fill array with 0 */
81        for(it = 0; it < nt; it++)
82        {
83            data[it] = 0.0;
84        }
85    }   
86    else
87    {
88        for(it = 0; it < nt; it++)
89        {
90            last_it = it;
91            /* skip t record of data block */
92            offset++;
93            /* skip records from earlier tide gauges */
94            for(jsta = 0; jsta < ista; jsta++)
95                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
96                    offset++;
97
98            /* deal with the tide gauge at hand */
99            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
100            {
101                /* gauge is recording at this time */
102                memcpy(data + it, muxData + offset, sizeof(float));
103               
104                //printf("%d: muxdata=%f\n", it, muxData[offset]);
105                //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);
106                offset++;
107            }
108            else if (it + 1 < nst[ista])
109            {
110                /* gauge has not yet started recording */
111                data[it] = 0.0;
112            }   
113            else
114                /* gauge has finished recording */
115            {
116                data[it] = NODATA;
117                break;
118            }
119
120            /* skip records from later tide gauges */
121            for(jsta = ista + 1; jsta < total_number_of_stations; jsta++)
122                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
123                    offset++;
124        }
125
126        if(last_it < nt - 1)
127            /* the loop was exited early because the gauge had
128            finished recording */
129            for(it = last_it+1; it < nt; it++)
130                data[it] = NODATA;
131    }
132} 
133
134
135char isdata(float x)
136{
137    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
138    {
139      return 0;
140    }
141    else
142    {
143      return 1; 
144    }
145}
146
147
148long getNumData(const int *fros, const int *lros, const int total_number_of_stations)
149/* calculates the number of data in the data block of a mux file */
150/* based on the first and last recorded output steps for each gauge */ 
151{
152    int ista, last_output_step;
153    long numData = 0;
154
155    last_output_step = 0;   
156    for(ista = 0; ista < total_number_of_stations; ista++)
157        if(*(fros + ista) != -1)
158        {
159            numData += *(lros + ista) - *(fros + ista) + 1;
160            last_output_step = (last_output_step < *(lros+ista) ? 
161                            *(lros+ista):last_output_step);
162        }   
163        numData += last_output_step*total_number_of_stations; /* these are the t records */
164        return numData;
165}
166
167/////////////////////////////////////////////////////////////////////////
168//Internal Functions
169int _read_mux2_headers(int numSrc, 
170                       char **muxFileNameArray, 
171                       int* total_number_of_stations,
172                       int* number_of_time_steps,
173                       double* delta_t,
174                       //long* numDataMax,
175                       int verbose)
176{
177    FILE *fp;
178    int numsta, i, j;
179    struct tgsrwg *mytgs=0;
180    char *muxFileName;                                                                 
181    char susMuxFileName;
182    long numData;
183    size_t elements_read; // fread return value
184    int block_size;
185
186    /* Allocate space for the names and the weights and pointers to the data*/
187
188    /* Check that the input files have mux2 extension*/
189    susMuxFileName = 0;
190    for(i = 0; i < numSrc; i++)
191    { 
192        muxFileName = muxFileNameArray[i];
193        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, 
194                     "mux2") != 0)
195        {
196            susMuxFileName = 1;
197            break;
198        }
199    }
200
201    if(susMuxFileName)
202    {
203        printf("\n**************************************************************************\n");
204        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
205        printf("   At least one input file name does not end with mux2\n");
206        printf("   Check your results carefully!\n");
207        printf("**************************************************************************\n\n");
208    }   
209
210    if (verbose)
211    {
212        printf("Reading mux header information\n");
213    }
214
215    // Loop over all sources, read headers and check compatibility
216    for (i = 0; i < numSrc; i++)
217    {
218        muxFileName = muxFileNameArray[i];
219
220        // Open the mux file
221        if((fp = fopen(muxFileName, "rb")) == NULL)
222        {
223            char *err_msg = strerror(errno);
224
225            fprintf(stderr, "cannot open file '%s': %s\n", muxFileName, err_msg);
226            return -1; 
227        }
228       
229        if (!i)
230        {
231            elements_read = fread(total_number_of_stations, sizeof(int), 1, fp);
232            if ((int) elements_read == 0 && ferror(fp)){
233              fprintf(stderr, "Error reading total number of stations\n");
234              return -2;
235            }
236       
237            fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
238            lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
239     
240            mytgs0 = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg));
241            mytgs = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg));
242
243            block_size = *total_number_of_stations*sizeof(struct tgsrwg);
244            elements_read = fread(mytgs0, block_size , 1, fp);
245            if ((int) elements_read == 0 && ferror(fp)){
246              fprintf(stderr, "Error reading mytgs0\n");
247              return -2;
248            }
249        }
250        else
251        {
252            // Check that the mux files are compatible
253            elements_read = fread(&numsta, sizeof(int), 1, fp);
254            if ((int) elements_read == 0 && ferror(fp)){
255              fprintf(stderr, "Error reading numsta\n");
256              return -2;
257            }
258           
259            if(numsta != *total_number_of_stations)
260            {
261                fprintf(stderr,"%s has different number of stations to %s\n", 
262                muxFileName, 
263                muxFileNameArray[0]);
264                fclose(fp);
265                return -1;   
266            }
267
268            block_size = numsta*sizeof(struct tgsrwg);
269            elements_read = fread(mytgs, block_size, 1, fp); 
270            if ((int) elements_read == 0 && ferror(fp)){
271              fprintf(stderr, "Error reading mgtgs\n");
272              return -2;
273            }       
274           
275           
276            for (j = 0; j < numsta; j++)
277            {
278                if (mytgs[j].dt != mytgs0[j].dt)
279                {
280                    fprintf(stderr, "%s has different sampling rate to %s\n", 
281                    muxFileName, 
282                    muxFileNameArray[0]);
283                    fclose(fp);
284                    return -1;           
285                }   
286                if (mytgs[j].nt != mytgs0[j].nt)
287                {
288                    fprintf(stderr, "%s has different series length to %s\n", 
289                    muxFileName, 
290                    muxFileNameArray[0]);
291                    fclose(fp);
292                    return -1;           
293                }
294
295                if (mytgs[j].nt != mytgs0[0].nt)
296                {
297                    printf("Station 0 has different series length to Station %d\n", j); 
298                }
299            }
300        }
301
302        /* Read the start and stop times for this source */
303        elements_read = fread(fros + i*(*total_number_of_stations), 
304                              *total_number_of_stations*sizeof(int), 1, fp);
305        if ((int) elements_read == 0 && ferror(fp)){
306          fprintf(stderr, "Error reading start times\n");
307          return -3;
308        }           
309                           
310                           
311        elements_read = fread(lros + i*(*total_number_of_stations), 
312                              *total_number_of_stations*sizeof(int), 1, fp);
313        if ((int) elements_read == 0 && ferror(fp)){
314          fprintf(stderr, "Error reading stop times\n");       
315          return -3;
316        }                     
317
318        /* Compute the size of the data block for this source */
319        numData = getNumData(fros + i*(*total_number_of_stations), 
320                             lros + i*(*total_number_of_stations), 
321                             (*total_number_of_stations));
322
323        /* Sanity check */
324        if (numData < 0)
325        {
326            fprintf(stderr,"Size of data block appears to be negative!\n");
327            return -1;       
328        }
329
330        if (numDataMax < numData)
331        {
332            numDataMax = numData;
333        }
334
335        fclose(fp);         
336    }
337
338   
339    // Store time resolution and number of timesteps   
340    // These are the same for all stations as tested above, so
341    // we take the first one.
342    *delta_t = (double)mytgs0[0].dt;
343    *number_of_time_steps = mytgs0[0].nt;
344
345    free(mytgs);
346
347    return 0; // Succesful execution
348}
349
350
351float** _read_mux2(int numSrc, 
352                   char **muxFileNameArray, 
353                   float *weights, 
354                   double *params, 
355                   int *number_of_stations,
356                   long *permutation,
357                   int verbose)
358{
359    FILE *fp;
360    int total_number_of_stations, i, isrc, ista, k;
361    char *muxFileName;
362    int istart=-1, istop=-1;
363    int number_of_selected_stations;
364    float *muxData=NULL; // Suppress warning
365    long numData;
366
367    int len_sts_data, error_code;
368    float **sts_data;
369    float *temp_sts_data;
370
371    long int offset;
372
373    int number_of_time_steps, N;
374    double delta_t;
375   
376    size_t elements_read;
377   
378    // Shorthands pointing to memory blocks for each source
379    int *fros_per_source=NULL;     
380    int *lros_per_source=NULL;         
381
382   
383    error_code = _read_mux2_headers(numSrc, 
384                                    muxFileNameArray, 
385                                    &total_number_of_stations,
386                                    &number_of_time_steps,
387                                    &delta_t,
388                                    verbose);
389    if (error_code != 0) {
390      printf("urs_ext.c: Internal function _read_mux2_headers failed: Error code = %d\n", 
391             error_code);
392     
393      return NULL;
394    }
395   
396
397    // Apply rule that an empty permutation file means 'take all stations'
398    // We could change this later by passing in None instead of the empty
399    // permutation.
400    number_of_selected_stations = *number_of_stations; 
401    if (number_of_selected_stations == 0)
402    {
403        number_of_selected_stations = total_number_of_stations; 
404   
405        // Return possibly updated number of stations
406        *number_of_stations = total_number_of_stations;     
407     
408        // Create the Identity permutation vector
409        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
410        for (i = 0; i < number_of_selected_stations; i++)
411        {
412            permutation[i] = (long) i; 
413        }
414    }
415   
416    // The params array is used only for passing data back to Python.
417    params[0] = (double) number_of_selected_stations;
418    params[1] = (double) delta_t;
419    params[2] = (double) number_of_time_steps;
420   
421   
422   
423    // Make array(s) to hold demuxed data for stations given in the
424    // permutation file
425    sts_data = (float**) malloc(number_of_selected_stations*sizeof(float*));
426    if (sts_data == NULL)
427    {
428        printf("ERROR: Memory for sts_data could not be allocated.\n");
429        return NULL;
430    }
431
432    // For each selected station, allocate space for its data
433    len_sts_data = number_of_time_steps + POFFSET; // Max length of each timeseries?
434    for (i = 0; i < number_of_selected_stations; i++)
435    {
436        // Initialise sts_data to zero
437        sts_data[i] = (float*) calloc(len_sts_data, sizeof(float));
438        if (sts_data[i] == NULL)
439        {
440            printf("ERROR: Memory for sts_data could not be allocated.\n");
441            return NULL;
442        }
443    }
444
445    temp_sts_data = (float*) calloc(len_sts_data, sizeof(float));
446
447    muxData = (float*) calloc(numDataMax, sizeof(float));
448   
449    // Loop over all sources
450    for (isrc = 0; isrc < numSrc; isrc++)
451    {
452   
453        // Shorthands to local memory
454        fros_per_source = (int*) fros + isrc*total_number_of_stations; 
455        lros_per_source = (int*) lros + isrc*total_number_of_stations;     
456           
457   
458        // Read in data block from mux2 file
459        muxFileName = muxFileNameArray[isrc];
460        if((fp = fopen(muxFileName, "rb")) == NULL)
461        {
462            fprintf(stderr, "cannot open file %s\n", muxFileName);
463            return NULL;                   
464        }
465
466        if (verbose){
467            printf("Reading mux file %s\n", muxFileName);
468        }
469
470        offset = (long int)sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
471        //printf("\n offset %i ", (long int)offset);
472        fseek(fp, offset, 0);
473
474        numData = getNumData(fros_per_source, 
475                             lros_per_source, 
476                             total_number_of_stations);
477        // Note numData is larger than what it has to be.                   
478        //elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp);
479        elements_read = fread(muxData, (size_t) sizeof(float), (size_t) numData, fp); 
480        //printf("\n elements_read  %d, ", (int)elements_read);
481        //printf("\n ferror(fp)  %d, ", (int)ferror(fp));
482        if ((int) elements_read == 0 && ferror(fp)) {
483          fprintf(stderr, "Error reading mux data\n");
484          return NULL;
485        }       
486       
487        fclose(fp); 
488       
489        // loop over stations present in the permutation array
490        //     use ista with mux data
491        //     use i with the processed data to be returned         
492        for(i = 0; i < number_of_selected_stations; i++)
493        {               
494   
495            ista = (int) permutation[i]; // Get global index into mux data 
496       
497            // fill the data0 array from the mux file, and weight it
498            fillDataArray(ista, 
499                          total_number_of_stations, 
500                          number_of_time_steps,
501                          mytgs0[ista].ig, // Grid number (if -1 fill with zeros)
502                          fros_per_source, 
503                          lros_per_source, 
504                          temp_sts_data, 
505                          &istart, 
506                          &istop, 
507                          muxData);
508
509            // Weight appropriately and add
510            for(k = 0; k < mytgs0[ista].nt; k++)
511            {
512                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
513                {
514                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
515                }
516                else
517                {
518                    sts_data[i][k] = NODATA;
519                }
520                //printf("%d: temp_sts_data[%d]=%f\n", i, k, temp_sts_data[k]);     
521           
522            }
523           
524           
525            // Update metadata (e.g. start time and end time)
526            N = number_of_time_steps;
527           
528            if (isrc == 0) {
529                // Assign values for first source
530                sts_data[i][N] = (float)mytgs0[ista].geolat;
531                sts_data[i][N+1] = (float)mytgs0[ista].geolon;
532                sts_data[i][N+2] = (float)mytgs0[ista].z;
533                sts_data[i][N+3] = (float)fros_per_source[ista];
534                sts_data[i][N+4] = (float)lros_per_source[ista];
535            } else {
536                // Update first and last timesteps for subsequent sources
537                if (sts_data[i][N+3] > (float)fros_per_source[ista]) {         
538                    if (verbose) {
539                        printf("Adjusting start time for station %d and source %d",
540                               ista, isrc);
541                        printf(" from %f to %f\n", 
542                               sts_data[i][N+3], 
543                               (float) fros_per_source[ista]); 
544                    }
545                    sts_data[i][N+3] = (float) fros_per_source[ista];
546                }
547               
548                if (sts_data[i][N+4] < (float) lros_per_source[ista]) {         
549                    if (verbose) {
550                        printf("Adjusting end time for station %d and source %d",
551                               ista, isrc);
552                        printf(" from %f to %f\n", 
553                               sts_data[i][N+4], 
554                               (float) lros_per_source[ista]); 
555                    }
556                    sts_data[i][N+4] = (float) lros_per_source[ista];
557                }               
558            }
559        }
560    }
561
562    //printf("sts_data[1,8]=%f\n", sts_data[1][8]);
563   
564    free(muxData);
565    free(temp_sts_data);
566    free(fros);
567    free(lros);
568    free(mytgs0);
569
570    return sts_data;
571}
572
573/////////////////////////////////////////////////////////////////////////
574//Python gateways
575PyObject *read_mux2(PyObject *self, PyObject *args)
576{
577    /*Read in mux 2 file
578
579    Python call:
580    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
581
582    NOTE:
583    A Python int is equivalent to a C long
584    (this becomes really important on 64 bit architectures)
585   
586    A Python double corresponds to a C double
587    */
588   
589    PyObject *filenames;
590    PyArrayObject *pyweights;
591    PyArrayObject *file_params;
592    PyArrayObject *permutation;  // Ordering of selected stations   
593    PyArrayObject *pydata;
594    PyObject *fname;
595
596    char **muxFileNameArray;
597    float **cdata;
598    float *weights;
599    int dimensions[2];
600    int numSrc;
601    int verbose;
602    int total_number_of_stations;
603    int number_of_selected_stations;   
604    int nt;
605    double dt;
606    int i;
607    int j;
608    int start_tstep;
609    int finish_tstep;
610    int it;
611    int time;
612    int num_ts;
613   
614    // Convert Python arguments to C
615    if (!PyArg_ParseTuple(args, "iOOOOi",
616              &numSrc, &filenames, &pyweights, &file_params, 
617              &permutation, &verbose)) 
618    {
619            PyErr_SetString(PyExc_RuntimeError, 
620                "Input arguments to read_mux2 failed");
621            return NULL;
622    }
623
624    if(!PyList_Check(filenames)) 
625    {
626        PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list");
627        return NULL;
628    }
629
630    if(PyList_Size(filenames) == 0)
631    {
632        PyErr_SetString(PyExc_ValueError, "empty lists not allowed");
633        return NULL;
634    }
635
636    if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) 
637    {
638        PyErr_SetString(PyExc_ValueError,
639            "pyweights must be one-dimensional and of type double");
640        return NULL; 
641    }
642
643    if(PyList_Size(filenames) != pyweights->dimensions[0])
644    {
645        PyErr_SetString(PyExc_ValueError, 
646            "Must specify one weight for each filename");
647        return NULL;
648    }
649
650    muxFileNameArray = (char**)malloc(numSrc*sizeof(char*));
651    if (muxFileNameArray == NULL) 
652    {
653        PyErr_SetString(PyExc_ValueError, 
654                        "ERROR: Memory for muxFileNameArray could not be allocated.");
655        return NULL;
656    }
657
658    for (i = 0; i < numSrc; i++)
659    {
660
661        fname = PyList_GetItem(filenames, i);
662        if (!fname) 
663        {
664            PyErr_SetString(PyExc_ValueError, "filename not a string");
665            return NULL;
666        }       
667
668        muxFileNameArray[i] = PyString_AsString(fname);
669        if (muxFileNameArray[i] == NULL) 
670        {
671            PyErr_SetString(PyExc_ValueError, 
672          "ERROR: Memory for muxFileNameArray could not be allocated.\n");
673            return NULL;
674        }
675    }
676
677    if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) 
678    {
679        PyErr_SetString(PyExc_ValueError,
680          "file_params must be one-dimensional and of type double");
681        return NULL; 
682    }
683
684
685    // Create array for weights which are passed to read_mux2
686    weights = (float*) malloc(numSrc*sizeof(float));
687    for (i = 0; i < numSrc; i++)
688    {
689        weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0]));
690    }
691   
692    // Desired number of stations
693    number_of_selected_stations = (int) permutation->dimensions[0];
694
695    // Read in mux2 data from file
696    cdata = _read_mux2(numSrc, 
697                       muxFileNameArray, 
698                       weights, 
699                       (double*)file_params->data,
700                       &number_of_selected_stations, 
701                       (long*) permutation->data,
702                       verbose);
703
704    if (!cdata) 
705    {
706        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
707        return NULL;
708    }       
709               
710               
711    // Allocate space for return vector
712    total_number_of_stations = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
713    dt = *(double*)(file_params->data + 1*file_params->strides[0]);
714    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
715
716   
717    // Find min and max start times of all gauges
718    start_tstep = nt + 1;
719    finish_tstep = -1;
720    for (i = 0; i < number_of_selected_stations; i++)
721    {
722        //printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
723        // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);   
724   
725        if ((int)cdata[i][nt + 3] < start_tstep)
726        {
727            start_tstep = (int)cdata[i][nt + 3];
728        }
729        if ((int)cdata[i][nt + 4] > finish_tstep)
730        {
731            finish_tstep = (int)cdata[i][nt + 4]; 
732        }
733    }
734
735    if ((start_tstep > nt) | (finish_tstep < 0))
736    {
737        printf("ERROR: Gauge data has incorrect start and finish times:\n");
738        printf("   start_tstep = %d, max_number_of_steps = %d\n", 
739               start_tstep, nt);
740        printf("   finish_tstep = %d, min_number_of_steps = %d\n", 
741               finish_tstep, 0);   
742           
743        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
744        return NULL;
745    }
746
747    if (start_tstep >= finish_tstep)
748    {
749        PyErr_SetString(PyExc_ValueError,
750                    "ERROR: Gauge data has non-postive_length");
751        return NULL;
752    }
753
754    num_ts = finish_tstep - start_tstep + 1;
755    dimensions[0] = number_of_selected_stations;
756    dimensions[1] = num_ts + POFFSET;
757   
758    pydata = (PyArrayObject*) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
759    if(pydata == NULL)
760    {
761        PyErr_SetString(PyExc_ValueError, 
762          "ERROR: Memory for pydata array could not be allocated.");
763        return NULL;
764    }
765
766   
767    // Each gauge begins and ends recording at different times. When a gauge is
768    // not recording but at least one other gauge is.
769    // Pad the non-recording gauge array with zeros.
770    //printf("\nData put into the cdata array from C code\n");
771    for (i = 0; i < number_of_selected_stations; i++)
772    {
773        time = 0;
774        for (it = 0; it < finish_tstep; it++)
775        {
776            if ((it + 1 >= start_tstep) && (it + 1 <= finish_tstep))
777            {
778                if (it + 1 > (int)cdata[i][nt + 4])
779                {
780                    // This gauge has stopped recording but others are
781                    // still recording
782                    *(double*)(pydata->data + i*pydata->strides[0] 
783                               + time*pydata->strides[1]) = 
784                      0.0;
785                }
786                else
787                {
788                  //printf("cdata[%d][%d] = %f\n", i, it, cdata[i][it]);
789               
790                    *(double*)(pydata->data + i*pydata->strides[0] 
791                                    + time*pydata->strides[1]) = 
792                        cdata[i][it];
793                }
794                time++;
795            }
796        }
797        // Pass back lat,lon,elevation
798        for (j = 0; j < POFFSET; j++)
799        {
800          *(double*)(pydata->data + i*pydata->strides[0] 
801                     + (num_ts + j)*pydata->strides[1]) = 
802            cdata[i][nt + j];
803        }
804    }
805
806    free(weights);
807   
808    // Free filename array, but not independent Python strings
809    // FIXME(Ole): Do we need to update a reference counter in this case?
810    free(muxFileNameArray);
811   
812    for (i = 0; i < number_of_selected_stations; ++i)
813    {
814        free(cdata[i]);
815    }
816    free(cdata);
817
818    return PyArray_Return(pydata);
819}
820
821//-------------------------------
822// Method table for python module
823//-------------------------------
824static struct PyMethodDef MethodTable[] = {
825    {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
826    {NULL, NULL}
827};
828
829// Module initialisation
830void initurs_ext(void){
831    Py_InitModule("urs_ext", MethodTable);
832
833    import_array(); // Necessary for handling of NumPY structures
834}
Note: See TracBrowser for help on using the repository browser.