Changeset 5612
- Timestamp:
- Aug 5, 2008, 3:46:33 PM (15 years ago)
- Location:
- anuga_core/source/anuga/shallow_water
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/data_manager.py
r5611 r5612 4888 4888 numSrc=len(filenames) 4889 4889 4890 file_params=-1*ones(3,Float) #[nsta,dt,nt]4890 file_params=-1*ones(3,Float) #[nsta,dt,nt] 4891 4891 4892 4892 # Convert verbose to int C flag … … 4896 4896 verbose=0 4897 4897 4898 4899 #msg = 'I got the permutation vector:' + str(permutation)4900 #assert permutation is not None, msg4901 4898 if permutation is None: 4902 4899 permutation = ensure_numeric([], Float) 4903 4900 4904 4901 # Call underlying C implementation urs2sts_ext.c 4905 data = read_mux2(numSrc, filenames,weights,file_params,permutation,verbose)4902 data = read_mux2(numSrc, filenames, weights, file_params, permutation, verbose) 4906 4903 4907 4904 msg='File parameter values were not read in correctly from c file' … … 4925 4922 msg='Must have at least one gauge value' 4926 4923 assert nt>0,msg 4927 4924 4928 4925 OFFSET=5 # Number of site parameters p passed back with data 4929 4926 # p=[geolat,geolon,depth,start_tstep,finish_tstep] 4930 4927 4928 # FIXME (Ole): What is the relationship with params and data.shape ? 4929 # It looks as if the following asserts should pass but they don't always 4930 # 4931 4932 #print 4933 #print nsta, nt, data.shape 4934 4935 #msg = 'nt = %d, data.shape[1] == %d' %(nt, data.shape[1]) 4936 #assert nt == data.shape[1] - OFFSET, msg 4937 4938 #msg = 'nsta = %d, data.shape[0] == %d' %(nsta, data.shape[0]) 4939 #assert nsta == data.shape[0], msg 4940 4941 4942 # Number of stations in ordering file 4943 number_of_selected_stations = data.shape[0] 4944 4945 # Index where data ends and parameters begin 4931 4946 parameters_index = data.shape[1]-OFFSET 4932 4947 4933 4948 times=dt*arange(parameters_index) 4934 latitudes=zeros( data.shape[0], Float)4935 longitudes=zeros( data.shape[0], Float)4936 elevation=zeros( data.shape[0], Float)4937 quantity=zeros(( data.shape[0], parameters_index), Float)4949 latitudes=zeros(number_of_selected_stations, Float) 4950 longitudes=zeros(number_of_selected_stations, Float) 4951 elevation=zeros(number_of_selected_stations, Float) 4952 quantity=zeros((number_of_selected_stations, parameters_index), Float) 4938 4953 4939 4954 4940 4955 starttime=1e16 4941 for i in range(0, data.shape[0]): 4956 for i in range(number_of_selected_stations): 4957 quantity[i][:]=data[i][:parameters_index] 4958 4942 4959 latitudes[i]=data[i][parameters_index] 4943 4960 longitudes[i]=data[i][parameters_index+1] 4944 4961 elevation[i]=-data[i][parameters_index+2] 4945 quantity[i][:]=data[i][:parameters_index] # Was data[i][:-OFFSET]4946 4947 4962 first_time_step = data[i][parameters_index+3] 4948 #print 'datamanager:', i, first_time_step, dt*first_time_step4949 starttime=min(dt* data[i][parameters_index+3],starttime)4963 4964 starttime=min(dt*first_time_step, starttime) 4950 4965 4951 4966 return times, latitudes, longitudes, elevation, quantity, starttime -
anuga_core/source/anuga/shallow_water/urs_ext.c
r5611 r5612 203 203 } 204 204 205 / * Loop over all sources, read headers and check compatibility */205 // Loop over all sources, read headers and check compatibility 206 206 for (i = 0; i < numSrc; i++) 207 207 { 208 208 muxFileName = muxFileNameArray[i]; 209 209 210 / * open the mux file */210 // Open the mux file 211 211 if((fp = fopen(muxFileName, "r")) == NULL) 212 212 { … … 219 219 fread(total_number_of_stations, sizeof(int), 1, fp); 220 220 221 fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));222 lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));221 fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 222 lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 223 223 224 mytgs0 = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg));225 mytgs = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg));224 mytgs0 = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg)); 225 mytgs = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg)); 226 226 227 227 fread(mytgs0, *total_number_of_stations*sizeof(struct tgsrwg), 1, fp); … … 229 229 else 230 230 { 231 /* check that the mux files are compatible */ 231 // Check that the mux files are compatible 232 232 fread(&numsta, sizeof(int), 1, fp); 233 233 if(numsta != *total_number_of_stations) … … 246 246 if (mytgs[j].dt != mytgs0[j].dt) 247 247 { 248 fprintf(stderr, "%s has different sampling rate to %s\n",248 fprintf(stderr, "%s has different sampling rate to %s\n", 249 249 muxFileName, 250 250 muxFileNameArray[0]); … … 254 254 if (mytgs[j].nt != mytgs0[j].nt) 255 255 { 256 fprintf(stderr, "%s has different series length to %s\n",256 fprintf(stderr, "%s has different series length to %s\n", 257 257 muxFileName, 258 258 muxFileNameArray[0]); … … 342 342 verbose); 343 343 344 //FIXME: The params can be removed after the python interface is modified.345 params[0] = (double)total_number_of_stations;346 params[1] = (double)delta_t;347 params[2] = (double)number_of_time_steps;348 349 344 // Apply rule that an empty permutation file means 'take all stations' 350 345 // We could change this later by passing in None instead of the empty … … 365 360 } 366 361 } 362 363 // The params array is used only for passing data back to Python. 364 params[0] = (double) number_of_selected_stations; 365 params[1] = (double) delta_t; 366 params[2] = (double) number_of_time_steps; 367 368 367 369 368 370 // Make array(s) to hold demuxed data for stations given in the 369 371 // permutation file 370 sts_data = (float**) malloc(number_of_selected_stations*sizeof(float*));372 sts_data = (float**) malloc(number_of_selected_stations*sizeof(float*)); 371 373 if (sts_data == NULL) 372 374 { … … 380 382 { 381 383 // Initialise sts_data to zero 382 sts_data[i] = (float*) calloc(len_sts_data, sizeof(float));384 sts_data[i] = (float*) calloc(len_sts_data, sizeof(float)); 383 385 if (sts_data[i] == NULL) 384 386 { … … 388 390 } 389 391 390 temp_sts_data = (float*) calloc(len_sts_data, sizeof(float));391 392 muxData = (float*) malloc(numDataMax*sizeof(float));392 temp_sts_data = (float*) calloc(len_sts_data, sizeof(float)); 393 394 muxData = (float*) malloc(numDataMax*sizeof(float)); 393 395 394 396 /* Loop over all sources */ … … 401 403 402 404 403 404 /* Read in data block from mux2 file */ 405 // Read in data block from mux2 file 405 406 muxFileName = muxFileNameArray[isrc]; 406 407 if((fp = fopen(muxFileName, "r")) == NULL) … … 427 428 // use ista with mux data 428 429 // use i with the processed data to be returned 429 430 430 for(i = 0; i < number_of_selected_stations; i++) 431 431 { … … 461 461 N = number_of_time_steps; 462 462 463 464 /*465 printf("station = %d, source = %d, fros = %f, lros = %f\n",466 ista,467 isrc,468 (float) fros_per_source[ista],469 (float) lros_per_source[ista]);470 */471 472 473 474 463 if (isrc == 0) { 475 464 // Assign values for first source … … 485 474 printf("Adjusting start time for station %d and source %d", 486 475 ista, isrc); 487 printf(" from %f to %f\n", sts_data[i][N+3], (float)fros_per_source[ista]); 476 printf(" from %f to %f\n", 477 sts_data[i][N+3], 478 (float) fros_per_source[ista]); 488 479 } 489 sts_data[i][N+3] = (float) fros_per_source[ista];480 sts_data[i][N+3] = (float) fros_per_source[ista]; 490 481 } 491 482 492 493 if (sts_data[i][N+4] < (float)lros_per_source[ista]) { 483 if (sts_data[i][N+4] < (float) lros_per_source[ista]) { 494 484 if (verbose) { 495 485 printf("Adjusting end time for station %d and source %d", 496 486 ista, isrc); 497 printf(" from %f to %f\n", sts_data[i][N+4], (float)lros_per_source[ista]); 487 printf(" from %f to %f\n", 488 sts_data[i][N+4], 489 (float) lros_per_source[ista]); 498 490 } 499 sts_data[i][N+4] = (float) lros_per_source[ista];491 sts_data[i][N+4] = (float) lros_per_source[ista]; 500 492 } 501 493 } … … 593 585 { 594 586 PyErr_SetString(PyExc_ValueError, 595 587 "ERROR: Memory for muxFileNameArray could not be allocated."); 596 588 return NULL; 597 589 } … … 697 689 dimensions[1] = num_ts + POFFSET; 698 690 699 pydata = (PyArrayObject*) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);691 pydata = (PyArrayObject*) PyArray_FromDims(2, dimensions, PyArray_DOUBLE); 700 692 if(pydata == NULL) 701 693 { … … 719 711 { 720 712 // This gauge has stopped recording but others are 721 // still recording713 // still recording 722 714 *(double*)(pydata->data + i*pydata->strides[0] 723 724 715 + time*pydata->strides[1]) = 716 0.0; 725 717 } 726 718 else
Note: See TracChangeset
for help on using the changeset viewer.