Changeset 7548 for anuga_core/source/anuga/shallow_water/urs_ext.c
- Timestamp:
- Oct 16, 2009, 11:24:26 AM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/urs_ext.c
r7276 r7548 34 34 35 35 36 36 */ 37 37 38 38 … … 42 42 void fillDataArray(int ista, int total_number_of_stations, int nt, int ig, int *nst, 43 43 int *nft, float *data, int *istart_p, 44 44 int *istop_p, float *muxData) 45 45 { 46 46 int it, last_it, jsta; … … 64 64 } 65 65 } 66 66 67 67 if (nft[ista] != -1) 68 68 { … … 76 76 } 77 77 } 78 78 79 79 if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */ 80 80 { … … 103 103 /* gauge is recording at this time */ 104 104 memcpy(data + it, muxData + offset, sizeof(float)); 105 106 107 105 106 //printf("%d: muxdata=%f\n", it, muxData[offset]); 107 //printf("data[%d]=%f, offset=%d\n", it, data[it], offset); 108 108 offset++; 109 109 } … … 139 139 if(x < NODATA + EPSILON && NODATA < x + EPSILON) 140 140 { 141 return 0;141 return 0; 142 142 } 143 143 else 144 144 { 145 return 1;145 return 1; 146 146 } 147 147 } … … 161 161 numData += *(lros + ista) - *(fros + ista) + 1; 162 162 last_output_step = (last_output_step < *(lros+ista) ? 163 163 *(lros+ista):last_output_step); 164 164 } 165 165 numData += last_output_step*total_number_of_stations; /* these are the t records */ … … 170 170 //Internal Functions 171 171 int _read_mux2_headers(int numSrc, 172 173 174 175 176 177 172 char **muxFileNameArray, 173 int* total_number_of_stations, 174 int* number_of_time_steps, 175 double* delta_t, 176 //long* numDataMax, 177 int verbose) 178 178 { 179 179 FILE *fp; … … 194 194 muxFileName = muxFileNameArray[i]; 195 195 if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, 196 196 "mux2") != 0) 197 197 { 198 198 susMuxFileName = 1; … … 228 228 return -1; 229 229 } 230 230 231 231 if (!i) 232 232 { 233 233 elements_read = fread(total_number_of_stations, sizeof(int), 1, fp); 234 if ((int) elements_read == 0 && ferror(fp)){ 235 fprintf(stderr, "Error reading total number of stations\n"); 236 return -2; 237 } 238 234 if ((int) elements_read == 0 && ferror(fp)){ 235 fprintf(stderr, "Error reading total number of stations\n"); 236 fclose(fp); 237 return -2; 238 } 239 239 240 fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 240 241 lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 241 242 242 243 mytgs0 = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg)); 243 244 mytgs = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg)); 244 245 245 246 block_size = *total_number_of_stations*sizeof(struct tgsrwg); 246 247 elements_read = fread(mytgs0, block_size , 1, fp); 247 if ((int) elements_read == 0 && ferror(fp)){ 248 fprintf(stderr, "Error reading mytgs0\n"); 249 return -2; 250 } 248 if ((int) elements_read == 0 && ferror(fp)){ 249 fprintf(stderr, "Error reading mytgs0\n"); 250 fclose(fp); 251 return -2; 252 } 251 253 } 252 254 else 253 255 { 254 256 // Check that the mux files are compatible 255 257 elements_read = fread(&numsta, sizeof(int), 1, fp); 256 if ((int) elements_read == 0 && ferror(fp)){ 257 fprintf(stderr, "Error reading numsta\n"); 258 return -2; 259 } 260 258 if ((int) elements_read == 0 && ferror(fp)){ 259 fprintf(stderr, "Error reading numsta\n"); 260 fclose(fp); 261 return -2; 262 } 263 261 264 if(numsta != *total_number_of_stations) 262 265 { 263 266 fprintf(stderr,"%s has different number of stations to %s\n", 264 muxFileName,265 muxFileNameArray[0]);267 muxFileName, 268 muxFileNameArray[0]); 266 269 fclose(fp); 267 270 return -1; 268 271 } 269 272 270 273 block_size = numsta*sizeof(struct tgsrwg); 271 274 elements_read = fread(mytgs, block_size, 1, fp); 272 if ((int) elements_read == 0 && ferror(fp)){ 273 fprintf(stderr, "Error reading mgtgs\n"); 274 return -2; 275 } 276 277 275 if ((int) elements_read == 0 && ferror(fp)){ 276 fprintf(stderr, "Error reading mgtgs\n"); 277 fclose(fp); 278 return -2; 279 } 280 281 278 282 for (j = 0; j < numsta; j++) 279 283 { … … 281 285 { 282 286 fprintf(stderr, "%s has different sampling rate to %s\n", 283 muxFileName,284 muxFileNameArray[0]);287 muxFileName, 288 muxFileNameArray[0]); 285 289 fclose(fp); 286 290 return -1; … … 289 293 { 290 294 fprintf(stderr, "%s has different series length to %s\n", 291 muxFileName,292 muxFileNameArray[0]);295 muxFileName, 296 muxFileNameArray[0]); 293 297 fclose(fp); 294 298 return -1; … … 304 308 /* Read the start and stop times for this source */ 305 309 elements_read = fread(fros + i*(*total_number_of_stations), 306 *total_number_of_stations*sizeof(int), 1, fp); 307 if ((int) elements_read == 0 && ferror(fp)){ 308 fprintf(stderr, "Error reading start times\n"); 309 return -3; 310 } 311 312 310 *total_number_of_stations*sizeof(int), 1, fp); 311 if ((int) elements_read == 0 && ferror(fp)){ 312 fprintf(stderr, "Error reading start times\n"); 313 fclose(fp); 314 return -3; 315 } 316 317 313 318 elements_read = fread(lros + i*(*total_number_of_stations), 314 *total_number_of_stations*sizeof(int), 1, fp); 315 if ((int) elements_read == 0 && ferror(fp)){ 316 fprintf(stderr, "Error reading stop times\n"); 317 return -3; 318 } 319 *total_number_of_stations*sizeof(int), 1, fp); 320 if ((int) elements_read == 0 && ferror(fp)){ 321 fprintf(stderr, "Error reading stop times\n"); 322 fclose(fp); 323 return -3; 324 } 319 325 320 326 /* Compute the size of the data block for this source */ 321 327 numData = getNumData(fros + i*(*total_number_of_stations), 322 323 328 lros + i*(*total_number_of_stations), 329 (*total_number_of_stations)); 324 330 325 331 /* Sanity check */ … … 338 344 } 339 345 340 346 341 347 // Store time resolution and number of timesteps 342 348 // These are the same for all stations as tested above, so … … 366 372 float *muxData=NULL; // Suppress warning 367 373 long numData; 374 long *perm = NULL; 375 long *permutation_temp = NULL; 368 376 369 377 int len_sts_data, error_code; … … 375 383 int number_of_time_steps, N; 376 384 double delta_t; 377 385 378 386 size_t elements_read; 379 387 380 388 // Shorthands pointing to memory blocks for each source 381 389 int *fros_per_source=NULL; 382 390 int *lros_per_source=NULL; 383 391 384 392 385 393 error_code = _read_mux2_headers(numSrc, 386 387 388 389 390 394 muxFileNameArray, 395 &total_number_of_stations, 396 &number_of_time_steps, 397 &delta_t, 398 verbose); 391 399 if (error_code != 0) { 392 printf("urs_ext.c: Internal function _read_mux2_headers failed: Error code = %d\n",393 394 395 return NULL;396 } 397 400 printf("urs_ext.c: Internal function _read_mux2_headers failed: Error code = %d\n", 401 error_code); 402 403 return NULL; 404 } 405 398 406 399 407 // Apply rule that an empty permutation file means 'take all stations' … … 404 412 { 405 413 number_of_selected_stations = total_number_of_stations; 406 414 407 415 // Return possibly updated number of stations 408 416 *number_of_stations = total_number_of_stations; 409 417 410 418 // Create the Identity permutation vector 411 permutation = (long *) malloc(number_of_selected_stations*sizeof(long)); 419 permutation_temp = (long *) malloc(number_of_selected_stations*sizeof(long)); 420 if (permutation_temp == NULL) 421 { 422 printf("ERROR: Memory for permutation_temp could not be allocated.\n"); 423 return NULL; 424 } 425 412 426 for (i = 0; i < number_of_selected_stations; i++) 413 427 { 414 permutation[i] = (long) i; 415 } 416 } 417 428 permutation_temp[i] = (long) i; 429 } 430 431 perm = permutation_temp; 432 } 433 else 434 { 435 perm = permutation; 436 } 437 418 438 // The params array is used only for passing data back to Python. 419 439 params[0] = (double) number_of_selected_stations; 420 440 params[1] = (double) delta_t; 421 441 params[2] = (double) number_of_time_steps; 422 423 424 442 425 443 // Make array(s) to hold demuxed data for stations given in the 426 444 // permutation file … … 446 464 447 465 temp_sts_data = (float*) calloc(len_sts_data, sizeof(float)); 466 if (temp_sts_data == NULL) 467 { 468 printf("ERROR: Memory for temp_sts_data could not be allocated.\n"); 469 return NULL; 470 } 448 471 449 472 muxData = (float*) calloc(numDataMax, sizeof(float)); 450 473 if (temp_sts_data == NULL) 474 { 475 printf("ERROR: Memory for muxData could not be allocated.\n"); 476 return NULL; 477 } 478 451 479 // Loop over all sources 452 480 for (isrc = 0; isrc < numSrc; isrc++) 453 481 { 454 482 455 483 // Shorthands to local memory 456 484 fros_per_source = (int*) fros + isrc*total_number_of_stations; 457 485 lros_per_source = (int*) lros + isrc*total_number_of_stations; 458 459 486 487 460 488 // Read in data block from mux2 file 461 489 muxFileName = muxFileNameArray[isrc]; … … 463 491 { 464 492 fprintf(stderr, "cannot open file %s\n", muxFileName); 493 free(muxData); 494 free(temp_sts_data); 495 free(muxData); 496 465 497 return NULL; 466 498 } … … 472 504 offset = (long int)sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int)); 473 505 //printf("\n offset %i ", (long int)offset); 474 506 fseek(fp, offset, 0); 475 507 476 508 numData = getNumData(fros_per_source, 477 478 479 509 lros_per_source, 510 total_number_of_stations); 511 // Note numData is larger than what it has to be. 480 512 //elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp); 481 513 elements_read = fread(muxData, (size_t) sizeof(float), (size_t) numData, fp); 482 //printf("\n elements_read %d, ", (int)elements_read); 483 //printf("\n ferror(fp) %d, ", (int)ferror(fp)); 484 if ((int) elements_read == 0 && ferror(fp)) { 485 fprintf(stderr, "Error reading mux data\n"); 486 return NULL; 487 } 488 514 //printf("\n elements_read %d, ", (int)elements_read); 515 //printf("\n ferror(fp) %d, ", (int)ferror(fp)); 516 if ((int) elements_read == 0 && ferror(fp)) { 517 fprintf(stderr, "Error reading mux data\n"); 518 fclose(fp); 519 free(muxData); 520 free(temp_sts_data); 521 free(muxData); 522 523 return NULL; 524 } 525 489 526 fclose(fp); 490 527 491 528 // loop over stations present in the permutation array 492 529 // use ista with mux data … … 494 531 for(i = 0; i < number_of_selected_stations; i++) 495 532 { 496 497 ista = (int) perm utation[i]; // Get global index into mux data498 533 534 ista = (int) perm[i]; // Get global index into mux data 535 499 536 // fill the data0 array from the mux file, and weight it 500 537 fillDataArray(ista, 501 502 503 504 505 506 507 508 509 538 total_number_of_stations, 539 number_of_time_steps, 540 mytgs0[ista].ig, // Grid number (if -1 fill with zeros) 541 fros_per_source, 542 lros_per_source, 543 temp_sts_data, 544 &istart, 545 &istop, 546 muxData); 510 547 511 548 // Weight appropriately and add … … 520 557 sts_data[i][k] = NODATA; 521 558 } 522 523 559 //printf("%d: temp_sts_data[%d]=%f\n", i, k, temp_sts_data[k]); 560 524 561 } 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 562 563 564 // Update metadata (e.g. start time and end time) 565 N = number_of_time_steps; 566 567 if (isrc == 0) { 568 // Assign values for first source 569 sts_data[i][N] = (float)mytgs0[ista].geolat; 570 sts_data[i][N+1] = (float)mytgs0[ista].geolon; 571 sts_data[i][N+2] = (float)mytgs0[ista].z; 572 sts_data[i][N+3] = (float)fros_per_source[ista]; 573 sts_data[i][N+4] = (float)lros_per_source[ista]; 574 } else { 575 // Update first and last timesteps for subsequent sources 576 if (sts_data[i][N+3] > (float)fros_per_source[ista]) { 577 if (verbose) { 578 printf("Adjusting start time for station %d and source %d", 579 ista, isrc); 580 printf(" from %f to %f\n", 581 sts_data[i][N+3], 582 (float) fros_per_source[ista]); 583 } 584 sts_data[i][N+3] = (float) fros_per_source[ista]; 585 } 586 587 if (sts_data[i][N+4] < (float) lros_per_source[ista]) { 588 if (verbose) { 589 printf("Adjusting end time for station %d and source %d", 590 ista, isrc); 591 printf(" from %f to %f\n", 592 sts_data[i][N+4], 593 (float) lros_per_source[ista]); 594 } 595 sts_data[i][N+4] = (float) lros_per_source[ista]; 596 } 597 } 561 598 } 562 599 } 563 600 564 601 //printf("sts_data[1,8]=%f\n", sts_data[1][8]); 565 602 566 603 free(muxData); 567 604 free(temp_sts_data); … … 570 607 free(mytgs0); 571 608 609 if (permutation_temp) 610 { 611 free(permutation_temp); 612 } 613 572 614 return sts_data; 573 615 } … … 585 627 A Python int is equivalent to a C long 586 628 (this becomes really important on 64 bit architectures) 587 629 588 630 A Python double corresponds to a C double 589 631 */ 590 632 591 633 PyObject *filenames; 592 634 PyArrayObject *pyweights; … … 613 655 int time; 614 656 int num_ts; 615 657 616 658 // Convert Python arguments to C 617 659 if (!PyArg_ParseTuple(args, "iOOOOi", 618 619 620 { 621 622 623 660 &numSrc, &filenames, &pyweights, &file_params, 661 &permutation, &verbose)) 662 { 663 PyErr_SetString(PyExc_RuntimeError, 664 "Input arguments to read_mux2 failed"); 665 return NULL; 624 666 } 625 667 … … 654 696 { 655 697 PyErr_SetString(PyExc_ValueError, 656 "ERROR: Memory for muxFileNameArray could not be allocated."); 698 "ERROR: Memory for muxFileNameArray could not be allocated."); 699 657 700 return NULL; 658 701 } … … 665 708 { 666 709 PyErr_SetString(PyExc_ValueError, "filename not a string"); 710 free(muxFileNameArray); 711 667 712 return NULL; 668 713 } … … 672 717 { 673 718 PyErr_SetString(PyExc_ValueError, 674 "ERROR: Memory for muxFileNameArray could not be allocated.\n"); 719 "ERROR: Memory for muxFileNameArray could not be allocated.\n"); 720 free(muxFileNameArray); 721 675 722 return NULL; 676 723 } … … 680 727 { 681 728 PyErr_SetString(PyExc_ValueError, 682 "file_params must be one-dimensional and of type double"); 729 "file_params must be one-dimensional and of type double"); 730 free(muxFileNameArray); 731 683 732 return NULL; 684 733 } … … 687 736 // Create array for weights which are passed to read_mux2 688 737 weights = (float*) malloc(numSrc*sizeof(float)); 738 if (weights == NULL) 739 { 740 PyErr_SetString(PyExc_ValueError, 741 "ERROR: Memory for weights could not be allocated."); 742 free(muxFileNameArray); 743 return NULL; 744 } 745 746 689 747 for (i = 0; i < numSrc; i++) 690 748 { 691 749 weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0])); 692 750 } 693 751 694 752 // Desired number of stations 695 753 number_of_selected_stations = (int) permutation->dimensions[0]; … … 697 755 // Read in mux2 data from file 698 756 cdata = _read_mux2(numSrc, 699 700 701 702 703 704 757 muxFileNameArray, 758 weights, 759 (double*)file_params->data, 760 &number_of_selected_stations, 761 (long*) permutation->data, 762 verbose); 705 763 706 764 if (!cdata) 707 765 { 708 PyErr_SetString(PyExc_ ValueError, "No STS_DATA returned");766 PyErr_SetString(PyExc_RuntimeError, "No STS_DATA returned"); 709 767 return NULL; 710 768 } 711 712 713 // Allocate space for return vector 714 total_number_of_stations = (int)*(double*)(file_params->data + 0*file_params->strides[0]); 769 770 //Extracting data 771 total_number_of_stations = (int)*((double*)(file_params->data + 0*file_params->strides[0])); 715 772 dt = *(double*)(file_params->data + 1*file_params->strides[0]); 716 773 nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]); 717 774 718 775 719 776 // Find min and max start times of all gauges 720 777 start_tstep = nt + 1; … … 724 781 //printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]); 725 782 // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]); 726 783 727 784 if ((int)cdata[i][nt + 3] < start_tstep) 728 785 { … … 739 796 printf("ERROR: Gauge data has incorrect start and finish times:\n"); 740 797 printf(" start_tstep = %d, max_number_of_steps = %d\n", 741 798 start_tstep, nt); 742 799 printf(" finish_tstep = %d, min_number_of_steps = %d\n", 743 744 800 finish_tstep, 0); 801 745 802 PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 803 804 free(weights); 805 free(muxFileNameArray); 806 for (i = 0; i < number_of_selected_stations; ++i) 807 { 808 free(cdata[i]); 809 } 810 free(cdata); 811 746 812 return NULL; 747 813 } … … 750 816 { 751 817 PyErr_SetString(PyExc_ValueError, 752 "ERROR: Gauge data has non-postive_length"); 818 "ERROR: Gauge data has non-postive_length"); 819 free(weights); 820 free(muxFileNameArray); 821 for (i = 0; i < number_of_selected_stations; ++i) 822 { 823 free(cdata[i]); 824 } 825 free(cdata); 826 753 827 return NULL; 754 828 } … … 757 831 dimensions[0] = number_of_selected_stations; 758 832 dimensions[1] = num_ts + POFFSET; 759 833 760 834 pydata = (PyArrayObject*) anuga_FromDims(2, dimensions, PyArray_DOUBLE); 761 835 if(pydata == NULL) 762 836 { 763 PyErr_SetString(PyExc_ValueError, 764 "ERROR: Memory for pydata array could not be allocated."); 765 return NULL; 766 } 767 768 837 PyErr_SetString(PyExc_RuntimeError, 838 "ERROR: Memory for pydata array could not be allocated."); 839 free(weights); 840 free(muxFileNameArray); 841 for (i = 0; i < number_of_selected_stations; ++i) 842 { 843 free(cdata[i]); 844 } 845 free(cdata); 846 847 return NULL; 848 } 849 850 769 851 // Each gauge begins and ends recording at different times. When a gauge is 770 852 // not recording but at least one other gauge is. … … 783 865 // still recording 784 866 *(double*)(pydata->data + i*pydata->strides[0] 785 786 867 + time*pydata->strides[1]) = 868 0.0; 787 869 } 788 870 else 789 871 { 790 791 872 //printf("cdata[%d][%d] = %f\n", i, it, cdata[i][it]); 873 792 874 *(double*)(pydata->data + i*pydata->strides[0] 793 875 + time*pydata->strides[1]) = 794 876 cdata[i][it]; 795 877 } … … 799 881 // Pass back lat,lon,elevation 800 882 for (j = 0; j < POFFSET; j++) 801 802 803 804 883 { 884 *(double*)(pydata->data + i*pydata->strides[0] 885 + (num_ts + j)*pydata->strides[1]) = 886 cdata[i][nt + j]; 805 887 } 806 888 } 807 889 808 890 free(weights); 809 891 810 892 // Free filename array, but not independent Python strings 811 893 // FIXME(Ole): Do we need to update a reference counter in this case? 812 894 free(muxFileNameArray); 813 895 814 896 for (i = 0; i < number_of_selected_stations; ++i) 815 897 {
Note: See TracChangeset
for help on using the changeset viewer.