Changeset 5516 for anuga_core/source/anuga/shallow_water/urs_ext.c
- Timestamp:
- Jul 17, 2008, 1:46:09 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/urs_ext.c
r5485 r5516 11 11 #include "Python.h" 12 12 #include "Numeric/arrayobject.h" 13 #include "structure.h" 13 14 #include "math.h" 14 15 #include <stdio.h> 15 16 #include <float.h> 16 17 #include <time.h> 17 #include "structure.h"18 18 19 19 #define MAX_FILE_NAME_LENGTH 128 … … 26 26 27 27 void fillDataArray(int, int, int, int, int *, int *, float *, int *, int *, float *); 28 long getNumData( int*, int*, int);28 long getNumData(const int *fros, const int *lros, const int nsta); 29 29 char isdata(float); 30 void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int); 31 float** _read_mux2(int, char **, float *, double *, int); 30 float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose); 32 31 33 32 PyObject *read_mux2(PyObject *self, PyObject *args){ 34 /*Read in mux 2 file35 33 /*Read in mux 2 file 34 36 35 Python call: 37 36 read_mux2(numSrc,filenames,weights,file_params,verbose) … … 40 39 A Python int is equivalent to a C long 41 40 A Python double corresponds to a C double 42 */ 43 PyObject *filenames; 44 PyArrayObject *pyweights,*file_params; 45 PyArrayObject *pydata; 46 PyObject *fname; 47 48 char **muxFileNameArray; 49 float *weights; 50 long numSrc; 51 long verbose; 52 53 float **cdata; 54 int dimensions[2]; 55 int nsta0; 56 int nt; 57 double dt; 58 int i,j; 59 int start_tstep; 60 int finish_tstep; 61 int it,time; 62 int num_ts; 63 64 // Convert Python arguments to C 65 if (!PyArg_ParseTuple(args, "iOOOi", 66 &numSrc,&filenames,&pyweights,&file_params,&verbose)) { 67 68 PyErr_SetString(PyExc_RuntimeError, 69 "Input arguments to read_mux2 failed"); 70 return NULL; 71 } 72 73 if(!PyList_Check(filenames)) { 74 PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list"); 75 return NULL; 76 } 77 78 if(PyList_Size(filenames) == 0){ 79 PyErr_SetString(PyExc_ValueError, "empty lists not allowed"); 80 return NULL; 81 } 82 83 if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) { 84 PyErr_SetString(PyExc_ValueError, 85 "pyweights must be one-dimensional and of type double"); 86 return NULL; 87 } 88 89 if(PyList_Size(filenames) != pyweights->dimensions[0]){ 90 PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename"); 91 return NULL; 92 } 93 94 muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *)); 95 if (muxFileNameArray == NULL) { 96 printf("ERROR: Memory for muxFileNameArray could not be allocated.\n"); 97 exit(-1); 98 } 99 for (i=0;i<PyList_Size(filenames);i++){ 100 muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char)); 101 if (muxFileNameArray[i] == NULL) { 102 printf("ERROR: Memory for muxFileNameArray could not be allocated.\n"); 103 exit(-1); 104 } 105 fname=PyList_GetItem(filenames, i); 106 if (!PyString_Check(fname)) { 107 PyErr_SetString(PyExc_ValueError, "filename not a string"); 108 } 109 muxFileNameArray[i]=PyString_AsString(fname); 110 } 111 112 if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) { 113 PyErr_SetString(PyExc_ValueError, 114 "file_params must be one-dimensional and of type double"); 115 return NULL; 116 } 117 118 // Create array for weights which are passed to read_mux2 119 weights = (float *) malloc((int)numSrc*sizeof(float)); 120 for (i=0;i<(int)numSrc;i++){ 121 weights[i]=(float)(*(double *)(pyweights->data+i*pyweights->strides[0])); 122 } 123 124 // Read in mux2 data from file 125 cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose); 126 127 128 // Allocate space for return vector 129 nsta0=(int)*(double *) (file_params -> data+0*file_params->strides[0]); 130 dt=*(double *) (file_params -> data+1*file_params->strides[0]); 131 nt=(int)*(double *) (file_params -> data+2*file_params->strides[0]); 132 133 134 // Find min and max start times of all gauges 135 start_tstep=nt+1; 136 finish_tstep=-1; 137 for (i=0;i<nsta0;i++){ 138 if ((int)cdata[i][nt+3] < start_tstep){ 139 start_tstep=(int)cdata[i][nt+3]; 140 } 141 if ((int)cdata[i][nt+4] > finish_tstep){ 142 finish_tstep=(int)cdata[i][nt+4]; 143 } 144 } 145 146 if ((start_tstep>nt) | (finish_tstep < 0)){ 147 printf("ERROR: Gauge data has incorrect start and finsh times\n"); 148 return NULL; 149 } 150 151 if (start_tstep>=finish_tstep){ 152 printf("ERROR: Gauge data has non-postive_length\n"); 153 return NULL; 154 } 155 156 num_ts=finish_tstep-start_tstep+1; 157 dimensions[0]=nsta0; 158 dimensions[1]=num_ts+POFFSET; 159 pydata = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE); 160 if(pydata == NULL){ 161 printf("ERROR: Memory for pydata array could not be allocated.\n"); 162 return NULL; 163 } 164 165 // Each gauge begins and ends recording at different times. When a gauge is 166 // not recording but at least one other gauge is. Pad the non-recording gauge 167 // array with zeros. 168 for (i=0;i<nsta0;i++){ 169 time=0; 170 for (it=0;it<finish_tstep;it++){ 171 if((it+1>=start_tstep)&&(it+1<=finish_tstep)){ 172 if (it+1>(int)cdata[i][nt+4]){ 173 //This gauge has stopped recording but others are still recording 174 *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=0.0; 175 }else{ 176 *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=cdata[i][it]; 177 } 178 time++; 179 } 180 } 181 // Pass back lat,lon,elevation 182 for (j=0; j<POFFSET; j++){ 183 *(double *) (pydata -> data+i*pydata->strides[0]+(num_ts+j)*pydata->strides[1])=cdata[i][nt+j]; 184 } 185 } 186 187 free(weights); 188 free(muxFileNameArray); 189 free(cdata); 190 return PyArray_Return(pydata); 191 41 */ 42 PyObject *filenames; 43 PyArrayObject *pyweights; 44 PyArrayObject *file_params; 45 PyArrayObject *pydata; 46 PyObject *fname; 47 48 char **muxFileNameArray; 49 float **cdata; 50 float *weights; 51 int dimensions[2]; 52 long numSrc; 53 long verbose; 54 int nsta0; 55 int nt; 56 double dt; 57 int i; 58 int j; 59 int start_tstep; 60 int finish_tstep; 61 int it; 62 int time; 63 int num_ts; 64 65 // Convert Python arguments to C 66 if (!PyArg_ParseTuple(args, "iOOOi", 67 &numSrc, &filenames, &pyweights, &file_params, &verbose)) 68 { 69 PyErr_SetString(PyExc_RuntimeError, 70 "Input arguments to read_mux2 failed"); 71 return NULL; 72 } 73 74 if(!PyList_Check(filenames)) 75 { 76 PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list"); 77 return NULL; 78 } 79 80 if(PyList_Size(filenames) == 0) 81 { 82 PyErr_SetString(PyExc_ValueError, "empty lists not allowed"); 83 return NULL; 84 } 85 86 if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) 87 { 88 PyErr_SetString(PyExc_ValueError, 89 "pyweights must be one-dimensional and of type double"); 90 return NULL; 91 } 92 93 if(PyList_Size(filenames) != pyweights->dimensions[0]) 94 { 95 PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename"); 96 return NULL; 97 } 98 99 muxFileNameArray = (char**)malloc((int)numSrc*sizeof(char*)); 100 if (muxFileNameArray == NULL) 101 { 102 printf("ERROR: Memory for muxFileNameArray could not be allocated.\n"); 103 exit(-1); 104 } 105 106 for (i = 0; i < numSrc; i++) 107 { 108 muxFileNameArray[i] = (char*)malloc((MAX_FILE_NAME_LENGTH + 1)*sizeof(char)); 109 110 if (muxFileNameArray[i] == NULL) 111 { 112 printf("ERROR: Memory for muxFileNameArray could not be allocated.\n"); 113 exit(-1); 114 } 115 116 fname = PyList_GetItem(filenames, i); 117 if (!PyString_Check(fname)) 118 { 119 PyErr_SetString(PyExc_ValueError, "filename not a string"); 120 } 121 122 muxFileNameArray[i] = PyString_AsString(fname); 123 } 124 125 if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) 126 { 127 PyErr_SetString(PyExc_ValueError, 128 "file_params must be one-dimensional and of type double"); 129 return NULL; 130 } 131 132 // Create array for weights which are passed to read_mux2 133 weights = (float*) malloc((int)numSrc*sizeof(float)); 134 for (i = 0; i < (int)numSrc; i++) 135 { 136 weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0])); 137 } 138 139 // Read in mux2 data from file 140 cdata = _read_mux2((int)numSrc, muxFileNameArray, weights, (double*)file_params->data, (int)verbose); 141 142 // Allocate space for return vector 143 nsta0 = (int)*(double*)(file_params->data + 0*file_params->strides[0]); 144 dt = *(double*)(file_params->data + 1*file_params->strides[0]); 145 nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]); 146 147 // Find min and max start times of all gauges 148 start_tstep = nt + 1; 149 finish_tstep = -1; 150 for (i = 0; i < nsta0; i++) 151 { 152 if ((int)cdata[i][nt + 3] < start_tstep) 153 { 154 start_tstep = (int)cdata[i][nt + 3]; 155 } 156 if ((int)cdata[i][nt + 4] > finish_tstep) 157 { 158 finish_tstep = (int)cdata[i][nt + 4]; 159 } 160 } 161 162 if ((start_tstep > nt) | (finish_tstep < 0)) 163 { 164 printf("ERROR: Gauge data has incorrect start and finsh times\n"); 165 return NULL; 166 } 167 168 if (start_tstep >= finish_tstep) 169 { 170 printf("ERROR: Gauge data has non-postive_length\n"); 171 return NULL; 172 } 173 174 num_ts = finish_tstep - start_tstep + 1; 175 dimensions[0] = nsta0; 176 dimensions[1] = num_ts + POFFSET; 177 pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE); 178 if(pydata == NULL) 179 { 180 printf("ERROR: Memory for pydata array could not be allocated.\n"); 181 return NULL; 182 } 183 184 // Each gauge begins and ends recording at different times. When a gauge is 185 // not recording but at least one other gauge is. Pad the non-recording gauge 186 // array with zeros. 187 for (i = 0; i < nsta0; i++) 188 { 189 time = 0; 190 for (it = 0; it < finish_tstep; it++) 191 { 192 if ((it + 1 >= start_tstep) && (it + 1 <= finish_tstep)) 193 { 194 if (it + 1 > (int)cdata[i][nt + 4]) 195 { 196 //This gauge has stopped recording but others are still recording 197 *(double*)(pydata->data + i*pydata->strides[0] + time*pydata->strides[1]) = 0.0; 198 } 199 else 200 { 201 *(double*)(pydata->data + i*pydata->strides[0] + time*pydata->strides[1]) = cdata[i][it]; 202 } 203 time++; 204 } 205 } 206 // Pass back lat,lon,elevation 207 for (j = 0; j < POFFSET; j++) 208 { 209 *(double*)(pydata->data + i*pydata->strides[0] + (num_ts + j)*pydata->strides[1]) = cdata[i][nt + j]; 210 } 211 } 212 213 free(weights); 214 215 for (i = 0; i < numSrc; ++i) 216 { 217 free(muxFileNameArray[i]); 218 } 219 free(muxFileNameArray); 220 221 for (i = 0; i < nsta0; ++i) 222 { 223 free(cdata[i]); 224 } 225 free(cdata); 226 227 return PyArray_Return(pydata); 192 228 } 193 229 194 230 float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose) 195 231 { 196 FILE *fp; 197 int nsta, nsta0, i, isrc, ista; 198 struct tgsrwg *mytgs, *mytgs0; 199 char *muxFileName; 200 int istart, istop; 201 int *fros, *lros; 202 char susMuxFileName; 203 float *muxData; 204 long numData; 205 206 207 int len_sts_data; 208 float **sts_data; 209 float *temp_sts_data; 210 211 long int offset; 212 213 /* Allocate space for the names and the weights and pointers to the data*/ 214 215 /* Check that the input files have mux2 extension*/ 216 susMuxFileName=0; 217 for(isrc=0;isrc<numSrc;isrc++) 218 { 219 muxFileName=muxFileNameArray[isrc]; 220 if(!susMuxFileName && strcmp(muxFileName+strlen(muxFileName)-4,"mux2")!=0){ 221 susMuxFileName=1; 222 } 223 } 224 225 if(susMuxFileName) 226 { 227 printf("\n**************************************************************************\n"); 228 printf(" WARNING: This program operates only on multiplexed files in mux2 format\n"); 229 printf(" At least one input file name does not end with mux2\n"); 230 printf(" Check your results carefully!\n"); 231 printf("**************************************************************************\n\n"); 232 } 233 234 if (verbose){ 235 printf("Reading mux header information\n"); 236 } 237 238 /* Open the first muxfile */ 239 if((fp=fopen(muxFileNameArray[0],"r"))==NULL){ 240 fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]); 241 exit(-1); 242 } 243 244 /* Read in the header */ 245 /* First read the number of stations*/ 246 fread(&nsta0,sizeof(int),1,fp); 247 248 /*now allocate space for, and read in, the structures for each station*/ 249 mytgs0 = (struct tgsrwg *) malloc(nsta0*sizeof(struct tgsrwg)); 250 fread(&mytgs0[0], nsta0*sizeof(struct tgsrwg), 1, fp); 251 252 /*make an array to hold the start and stop steps for each station for each 253 source*/ 254 //FIXME: Should only store start and stop times for one source 255 fros = (int *) malloc(nsta0*numSrc*sizeof(int)); 256 lros = (int *) malloc(nsta0*numSrc*sizeof(int)); 257 258 /* read the start and stop steps for source 0 into the array */ 259 fread(fros,nsta0*sizeof(int),1,fp); 260 fread(lros,nsta0*sizeof(int),1,fp); 261 262 /* compute the size of the data block for source 0 */ 263 numData = getNumData(fros, lros, nsta0); 264 265 /* Burbidge: Added a sanity check here */ 266 if (numData < 0) { 267 fprintf(stderr,"Size of data block appears to be negative!\n"); 268 exit(-1); 269 } 270 fclose(fp); 271 272 // Allocate header array for each remaining source file. 273 // FIXME: only need to store for one source and which is compared to all subsequent 274 // source headers 275 if(numSrc > 1){ 276 /* allocate space for tgsrwg for the other sources */ 277 mytgs = (struct tgsrwg *)malloc( nsta0*sizeof(struct tgsrwg) ); 278 } else { 279 /* FIXME (JJ): What should happen in case the are no source files?*/ 280 /* If we exit here, tests will break */ 281 // fprintf(stderr, "No source file specified\n"); 282 // exit(-1); 283 } 284 285 /* loop over remaining sources, check compatibility, and read them into 286 *muxData */ 287 for(isrc=1; isrc<numSrc; isrc++){ 288 muxFileName = muxFileNameArray[isrc]; 289 290 /* open the mux file */ 291 if((fp=fopen(muxFileName,"r"))==NULL) 292 { 293 fprintf(stderr, "cannot open file %s\n", muxFileName); 294 exit(-1); 295 } 296 297 /* check that the mux files are compatible */ 298 fread(&nsta,sizeof(int),1,fp); 299 if(nsta != nsta0){ 300 fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]); 301 fclose(fp); 302 exit(-1); 303 } 304 fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp); 305 for(ista=0; ista < nsta; ista++){ 306 if(mytgs[ista].dt != mytgs0[ista].dt){ 307 fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]); 308 fclose(fp); 309 exit(-1); 310 } 311 if(mytgs[ista].nt != mytgs0[ista].nt){ 312 fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]); 313 fclose(fp); 314 exit(-1); 315 } 316 } 317 318 /* read the start and stop times for this source */ 319 fread(fros+isrc*nsta0,nsta0*sizeof(int),1,fp); 320 fread(lros+isrc*nsta0,nsta0*sizeof(int),1,fp); 232 FILE *fp; 233 int nsta, nsta0, i, isrc, ista; 234 struct tgsrwg *mytgs=0, *mytgs0=0; 235 char *muxFileName; 236 int istart, istop; 237 int *fros=0, *lros=0; 238 char susMuxFileName; 239 float *muxData; 240 long numData; 241 242 int len_sts_data; 243 float **sts_data; 244 float *temp_sts_data; 245 246 long int offset; 247 248 /* Allocate space for the names and the weights and pointers to the data*/ 249 250 /* Check that the input files have mux2 extension*/ 251 susMuxFileName = 0; 252 for(isrc = 0; isrc < numSrc; isrc++) 253 { 254 muxFileName = muxFileNameArray[isrc]; 255 if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, "mux2") != 0) 256 { 257 susMuxFileName = 1; 258 break; 259 } 260 } 261 262 if(susMuxFileName) 263 { 264 printf("\n**************************************************************************\n"); 265 printf(" WARNING: This program operates only on multiplexed files in mux2 format\n"); 266 printf(" At least one input file name does not end with mux2\n"); 267 printf(" Check your results carefully!\n"); 268 printf("**************************************************************************\n\n"); 269 } 270 271 if (verbose) 272 { 273 printf("Reading mux header information\n"); 274 } 275 276 /* loop over remaining sources, check compatibility, and read them into 277 *muxData */ 278 for (isrc = 0; isrc < numSrc; isrc++) 279 { 280 muxFileName = muxFileNameArray[isrc]; 281 282 /* open the mux file */ 283 if((fp = fopen(muxFileName, "r")) == NULL) 284 { 285 fprintf(stderr, "cannot open file %s\n", muxFileName); 286 exit(-1); 287 } 288 289 if (!isrc) 290 { 291 fread(&nsta0, sizeof(int), 1, fp); 292 293 fros = (int*)malloc(nsta0*numSrc*sizeof(int)); 294 lros = (int*)malloc(nsta0*numSrc*sizeof(int)); 321 295 322 /* compute the size of the data block for this source */ 323 numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0); 324 325 /* Burbidge: Added a sanity check here */ 326 if (numData < 0){ 327 fprintf(stderr,"Size of data block appears to be negative!\n"); 328 exit(-1); 329 } 330 fclose(fp); 331 } 332 params[0]=(double)nsta0; 333 params[1]=(double)mytgs0[0].dt; 334 params[2]=(double)mytgs0[0].nt; 335 336 /* make array(s) to hold the demuxed data */ 337 //FIXME: This can be reduced to only contain stations in order file 338 sts_data = (float **)malloc (nsta0 * sizeof(float *)); 339 340 if (sts_data == NULL){ 341 printf("ERROR: Memory for sts_data could not be allocated.\n"); 342 exit(-1); 343 } 344 345 len_sts_data=mytgs0[0].nt + POFFSET; 346 for (ista=0; ista<nsta0; ista++){ 347 // Initialise sts_data to zero 348 sts_data[ista] = (float *)calloc(len_sts_data, sizeof(float) ); 349 if (sts_data[ista] == NULL){ 350 printf("ERROR: Memory for sts_data could not be allocated.\n"); 351 exit(-1); 352 } 353 sts_data[ista][mytgs0[0].nt]=(float)mytgs0[ista].geolat; 354 sts_data[ista][mytgs0[0].nt+1]=(float)mytgs0[ista].geolon; 355 sts_data[ista][mytgs0[0].nt+2]=(float)mytgs0[ista].z; 356 sts_data[ista][mytgs0[0].nt+3]=(float)fros[ista]; 357 sts_data[ista][mytgs0[0].nt+4]=(float)lros[ista]; 358 } 359 360 temp_sts_data = (float *)calloc(len_sts_data, sizeof(float) ); 361 362 /* Loop over all sources */ 363 //FIXME: remove istart and istop they are not used. 364 istart = -1; 365 istop = -1; 366 for (isrc=0;isrc<numSrc;isrc++){ 367 /* Read in data block from mux2 file */ 368 muxFileName = muxFileNameArray[isrc]; 369 if((fp=fopen(muxFileName,"r"))==NULL) 370 { 371 //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName); 372 fprintf(stderr, "cannot open file %s\n", muxFileName); 373 exit(-1); 374 } 375 376 if (verbose){ 377 printf("Reading mux file %s\n",muxFileName); 378 } 379 380 offset=sizeof(int)+nsta0*(sizeof(struct tgsrwg)+2*sizeof(int)); 381 fseek(fp,offset,0); 382 383 numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0); 384 muxData = (float*) malloc(numData*sizeof(float)); 385 fread(muxData, numData*sizeof(float),1,fp); 386 fclose(fp); 387 388 /* loop over stations */ 389 for(ista=0; ista<nsta0; ista++){ 390 391 /* fill the data0 array from the mux file, and weight it */ 392 fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData); 393 394 /* weight appropriately and add */ 395 for(i=0; i<mytgs0[ista].nt; i++){ 396 if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i])){ 397 sts_data[ista][i] += temp_sts_data[i] * weights[isrc]; 398 //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]); 399 }else{ 400 sts_data[ista][i] = NODATA; 401 } 402 } 403 } 404 } 405 free(muxData); 406 free(temp_sts_data); 407 return sts_data; 408 } 296 mytgs0 = (struct tgsrwg*)malloc(nsta0*sizeof(struct tgsrwg)); 297 mytgs = (struct tgsrwg*)malloc(nsta0*sizeof(struct tgsrwg)); 298 299 fread(mytgs0, nsta0*sizeof(struct tgsrwg), 1, fp); 300 } 301 else 302 { 303 /* check that the mux files are compatible */ 304 fread(&nsta, sizeof(int), 1, fp); 305 if(nsta != nsta0) 306 { 307 fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]); 308 fclose(fp); 309 exit(-1); 310 } 311 312 fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp); 313 314 for (ista = 0; ista < nsta; ista++) 315 { 316 if (mytgs[ista].dt != mytgs0[ista].dt) 317 { 318 fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]); 319 fclose(fp); 320 exit(-1); 321 } 322 if (mytgs[ista].nt != mytgs0[ista].nt) 323 { 324 fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]); 325 fclose(fp); 326 exit(-1); 327 } 328 } 329 } 330 331 /* read the start and stop times for this source */ 332 fread(fros + isrc*nsta0, nsta0*sizeof(int), 1, fp); 333 fread(lros + isrc*nsta0, nsta0*sizeof(int), 1, fp); 334 335 /* compute the size of the data block for this source */ 336 numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0); 337 338 /* Burbidge: Added a sanity check here */ 339 if (numData < 0) 340 { 341 fprintf(stderr,"Size of data block appears to be negative!\n"); 342 exit(-1); 343 } 344 345 fclose(fp); 346 } 347 348 params[0] = (double)nsta0; 349 params[1] = (double)mytgs0[0].dt; 350 params[2] = (double)mytgs0[0].nt; 351 352 /* make array(s) to hold the demuxed data */ 353 //FIXME: This can be reduced to only contain stations in order file 354 sts_data = (float**)malloc(nsta0*sizeof(float *)); 355 if (sts_data == NULL) 356 { 357 printf("ERROR: Memory for sts_data could not be allocated.\n"); 358 exit(-1); 359 } 360 361 len_sts_data = mytgs0[0].nt + POFFSET; 362 for (ista = 0; ista < nsta0; ista++) 363 { 364 // Initialise sts_data to zero 365 sts_data[ista] = (float*)calloc(len_sts_data, sizeof(float)); 366 if (sts_data[ista] == NULL) 367 { 368 printf("ERROR: Memory for sts_data could not be allocated.\n"); 369 exit(-1); 370 } 371 372 sts_data[ista][mytgs0[0].nt] = (float)mytgs0[ista].geolat; 373 sts_data[ista][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon; 374 sts_data[ista][mytgs0[0].nt + 2] = (float)mytgs0[ista].z; 375 sts_data[ista][mytgs0[0].nt + 3] = (float)fros[ista]; 376 sts_data[ista][mytgs0[0].nt + 4] = (float)lros[ista]; 377 } 378 379 temp_sts_data = (float*)calloc(len_sts_data, sizeof(float)); 380 381 /* Loop over all sources */ 382 //FIXME: remove istart and istop they are not used. 383 istart = -1; 384 istop = -1; 385 for (isrc = 0; isrc < numSrc; isrc++) 386 { 387 /* Read in data block from mux2 file */ 388 muxFileName = muxFileNameArray[isrc]; 389 if((fp = fopen(muxFileName, "r")) == NULL) 390 { 391 //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName); 392 fprintf(stderr, "cannot open file %s\n", muxFileName); 393 exit(-1); 394 } 395 396 if (verbose){ 397 printf("Reading mux file %s\n",muxFileName); 398 } 399 400 offset = sizeof(int) + nsta0*(sizeof(struct tgsrwg) + 2*sizeof(int)); 401 fseek(fp, offset, 0); 402 403 numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0); 404 muxData = (float*)malloc(numData*sizeof(float)); 405 fread(muxData, numData*sizeof(float), 1, fp); 406 fclose(fp); 407 408 /* loop over stations */ 409 for(ista = 0; ista < nsta0; ista++) 410 { 411 /* fill the data0 array from the mux file, and weight it */ 412 fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros + isrc*nsta0, lros + isrc*nsta0, temp_sts_data, &istart, &istop, muxData); 413 414 /* weight appropriately and add */ 415 for(i = 0; i < mytgs0[ista].nt; i++) 416 { 417 if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i])) 418 { 419 sts_data[ista][i] += temp_sts_data[i] * weights[isrc]; 420 //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]); 421 } 422 else 423 { 424 sts_data[ista][i] = NODATA; 425 } 426 } 427 } 428 } 429 430 free(muxData); 431 free(temp_sts_data); 432 free(fros); 433 free(lros); 434 free(mytgs0); 435 free(mytgs); 436 437 return sts_data; 438 } 409 439 410 440 … … 412 442 void fillDataArray(int ista, int nsta, int nt, int ig, int *nst, int *nft, float *data, int *istart_p, int *istop_p, float *muxData) 413 443 { 414 int it, last_it, jsta; 415 long int offset=0; 416 417 418 last_it=-1; 419 /* make arrays of starting and finishing time steps for the tide gauges */ 420 /* and fill them from the file */ 421 422 /* update start and stop timesteps for this gauge */ 423 if(nst[ista]!=-1){ 424 if(*istart_p==-1){ 425 *istart_p=nst[ista]; 426 }else{ 427 *istart_p=((nst[ista]<*istart_p)?nst[ista]:*istart_p); 428 } 429 } 430 if(nft[ista]!=-1){ 431 if(*istop_p==-1){ 432 *istop_p=nft[ista]; 433 }else{ 434 *istop_p=((nft[ista]<*istop_p)?nft[ista]:*istop_p); 435 } 436 } 437 if(ig==-1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */ 438 { 439 /* gauge never started recording, or was outside of all grids, fill array with 0 */ 440 for(it=0; it<nt; it++) 441 data[it] = 0.0; 442 } 443 else 444 { 445 for(it=0; it<nt; it++) 446 { 447 last_it = it; 448 /* skip t record of data block */ 449 offset++; 450 /* skip records from earlier tide gauges */ 451 for(jsta=0; jsta<ista; jsta++) 452 if(it+1>=nst[jsta]&&it+1<=nft[jsta]) 453 offset++; 454 455 /* deal with the tide gauge at hand */ 456 if(it+1>=nst[ista]&&it+1<=nft[ista]) 457 /* gauge is recording at this time */ 458 { 459 memcpy(data+it,muxData+offset,sizeof(float)); 444 int it, last_it, jsta; 445 long int offset=0; 446 447 448 last_it = -1; 449 /* make arrays of starting and finishing time steps for the tide gauges */ 450 /* and fill them from the file */ 451 452 /* update start and stop timesteps for this gauge */ 453 if (nst[ista]!= -1) 454 { 455 if(*istart_p == -1) 456 { 457 *istart_p = nst[ista]; 458 } 459 else 460 { 461 *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p); 462 } 463 } 464 if (nft[ista] != -1) 465 { 466 if (*istop_p == -1) 467 { 468 *istop_p = nft[ista]; 469 } 470 else 471 { 472 *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p); 473 } 474 } 475 if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */ 476 { 477 /* gauge never started recording, or was outside of all grids, fill array with 0 */ 478 for(it = 0; it < nt; it++) 479 { 480 data[it] = 0.0; 481 } 482 } 483 else 484 { 485 for(it = 0; it < nt; it++) 486 { 487 last_it = it; 488 /* skip t record of data block */ 460 489 offset++; 461 } 462 else if (it+1<nst[ista]) 463 { 464 /* gauge has not yet started recording */ 465 data[it] = 0.0; 466 } 467 else 468 /* gauge has finished recording */ 469 { 470 data[it] = NODATA; 471 break; 472 } 473 474 /* skip records from later tide gauges */ 475 for(jsta=ista+1; jsta<nsta; jsta++) 476 if(it+1>=nst[jsta]&&it+1<=nft[jsta]) 477 offset++; 478 } 479 480 if(last_it < nt - 1) 481 /* the loop was exited early because the gauge had finished recording */ 482 for(it=last_it+1; it < nt; it++) 483 data[it] = NODATA; 484 } 490 /* skip records from earlier tide gauges */ 491 for(jsta = 0; jsta < ista; jsta++) 492 if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta]) 493 offset++; 494 495 /* deal with the tide gauge at hand */ 496 if(it + 1 >= nst[ista] && it + 1 <= nft[ista]) 497 /* gauge is recording at this time */ 498 { 499 memcpy(data + it, muxData + offset, sizeof(float)); 500 offset++; 501 } 502 else if (it + 1 < nst[ista]) 503 { 504 /* gauge has not yet started recording */ 505 data[it] = 0.0; 506 } 507 else 508 /* gauge has finished recording */ 509 { 510 data[it] = NODATA; 511 break; 512 } 513 514 /* skip records from later tide gauges */ 515 for(jsta = ista + 1; jsta < nsta; jsta++) 516 if(it + 1 >= nst[jsta] && it+1 <= nft[jsta]) 517 offset++; 518 } 519 520 if(last_it < nt - 1) 521 /* the loop was exited early because the gauge had finished recording */ 522 for(it = last_it+1; it < nt; it++) 523 data[it] = NODATA; 524 } 485 525 } 486 526 487 /* Burbidge: No "offset" is sent. Replace with max. Added grid_id */488 void wrttxt( fname, dt, nt, x, beg, lat, lon, max, depth, grid_id, istart, istop )489 char *fname;490 float dt, *x, beg, max, lat, lon, depth;491 int grid_id;492 int nt;493 int istart, istop;494 {495 int it;496 float t;497 FILE *fp;498 499 fp = fopen(fname,"w");500 fprintf(fp,"> lat %g lon %g depth %g max %g start_time %g stop_time %g grid_id %d\n", lat, lon, depth, max,501 istart*dt, istop*dt, grid_id );502 for(it=0;it<nt;it++)503 {504 t=beg+it*dt;505 fprintf(fp,"%9.3f %g\n", t, x[it]);506 }507 fclose(fp);508 }509 510 527 char isdata(float x) 511 528 { 512 //char value; 513 if(x < NODATA + EPSILON && NODATA < x + EPSILON) 514 return 0; 515 else 516 return 1; 529 //char value; 530 if(x < NODATA + EPSILON && NODATA < x + EPSILON) 531 { 532 return 0; 533 } 534 else 535 { 536 return 1; 537 } 517 538 } 518 539 519 long getNumData(int *fros, int *lros, int nsta) 540 541 long getNumData(const int *fros, const int *lros, const int nsta) 520 542 /* calculates the number of data in the data block of a mux file */ 521 543 /* based on the first and last recorded output steps for each gauge */ 522 544 { 523 int ista, last_output_step;524 long numData = 0;525 526 last_output_step = 0;527 for(ista=0; ista < nsta; ista++)528 if(*(fros + ista) != -1)529 {530 numData += *(lros + ista) - *(fros + ista) + 1;531 last_output_step = (last_output_step < *(lros+ista) ? *(lros+ista):last_output_step);532 }533 numData += last_output_step*nsta; /* these are the t records */534 return numData;545 int ista, last_output_step; 546 long numData = 0; 547 548 last_output_step = 0; 549 for(ista = 0; ista < nsta; ista++) 550 if(*(fros + ista) != -1) 551 { 552 numData += *(lros + ista) - *(fros + ista) + 1; 553 last_output_step = (last_output_step < *(lros+ista) ? *(lros+ista):last_output_step); 554 } 555 numData += last_output_step*nsta; /* these are the t records */ 556 return numData; 535 557 } 536 558 … … 541 563 //------------------------------- 542 564 static struct PyMethodDef MethodTable[] = { 543 {"read_mux2", read_mux2, METH_VARARGS, "Print out"},544 {NULL, NULL}565 {"read_mux2", read_mux2, METH_VARARGS, "Print out"}, 566 {NULL, NULL} 545 567 }; 546 568 547 569 // Module initialisation 548 570 void initurs_ext(void){ 549 Py_InitModule("urs_ext", MethodTable);550 551 import_array(); // Necessary for handling of NumPY structures571 Py_InitModule("urs_ext", MethodTable); 572 573 import_array(); // Necessary for handling of NumPY structures 552 574 }
Note: See TracChangeset
for help on using the changeset viewer.