Changeset 5537 for anuga_core/source/anuga/shallow_water/data_manager.py
- Timestamp:
- Jul 18, 2008, 4:59:46 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/data_manager.py
r5534 r5537 4847 4847 """Access the mux files specified in the filenames list. Combine the 4848 4848 data found therin as a weighted linear sum as specifed by the weights. 4849 If permutation is None extract timeseries data for all gauges within the4849 If permutation is None or empty extract timeseries data for all gauges within the 4850 4850 files. 4851 4851 … … 4871 4871 verbose=0 4872 4872 4873 4874 #msg = 'I got the permutation vector:' + str(permutation) 4875 #assert permutation is not None, msg 4876 if permutation is None: 4877 permutation = ensure_numeric([], Float) 4878 4873 4879 # Call underlying C implementation urs2sts_ext.c 4874 data=read_mux2(numSrc,filenames,weights,file_params, verbose)4880 data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose) 4875 4881 4876 4882 msg='File parameter values were not read in correctly from c file' 4877 4883 assert len(compress(file_params>0,file_params))!=0,msg 4878 msg='The number of stations specifed in the c array and in the file are inconsitent' 4879 assert file_params[0]==data.shape[0],msg 4884 4885 msg='The number of stations specifed in the c array and in the file are inconsistent' 4886 assert file_params[0] >= len(permutation) 4887 4888 msg='The number of stations returned is inconsistent with the requested number' 4889 assert len(permutation) == 0 or len(permutation) == data.shape[0], msg 4880 4890 4881 4891 nsta=int(file_params[0]) 4882 4892 msg='Must have at least one station' 4883 4893 assert nsta>0,msg 4894 4884 4895 dt=file_params[1] 4885 4896 msg='Must have a postive timestep' 4886 4897 assert dt>0,msg 4898 4887 4899 nt=int(file_params[2]) 4888 4900 msg='Must have at least one gauge value' … … 4962 4974 other lines: index,longitude,latitude 4963 4975 4964 4965 If ordering is Noneall points are taken in the order they4976 If ordering is None or ordering file is empty then 4977 all points are taken in the order they 4966 4978 appear in the mux2 file. 4967 4979 … … 5021 5033 raise IOError, msg 5022 5034 5023 # Read MUX2 files 5024 if (verbose): print 'reading mux2 file' 5025 mux={} 5026 for i, quantity in enumerate(quantities): 5027 5028 # For each quantity read the associated list of source mux2 file with 5029 # extention associated with that quantity 5030 times,\ 5031 latitudes_urs,\ 5032 longitudes_urs,\ 5033 elevation,\ 5034 mux[quantity],\ 5035 starttime = read_mux2_py(files_in[i], weights, verbose=verbose) 5036 5037 # Check that all quantities have consistent time and space information 5038 if quantity!=quantities[0]: 5039 msg='%s, %s and %s have inconsistent gauge data'\ 5040 %(files_in[0], files_in[1], files_in[2]) 5041 assert allclose(times, times_old), msg 5042 assert allclose(latitudes_urs, latitudes_urs_old), msg 5043 assert allclose(longitudes_urs, longitudes_urs_old), msg 5044 assert allclose(elevation, elevation_old), msg 5045 assert allclose(starttime, starttime_old), msg 5046 times_old=times 5047 latitudes_urs_old=latitudes_urs 5048 longitudes_urs_old=longitudes_urs 5049 elevation_old=elevation 5050 starttime_old=starttime 5051 5052 5035 # Establish permutation array 5053 5036 if ordering_filename is not None: 5037 5054 5038 # Read ordering file 5055 5056 5039 try: 5057 5040 fid=open(ordering_filename, 'r') … … 5074 5057 raise Exception, msg 5075 5058 5076 5077 5078 permutation = [int(line.split(',')[0]) for line in ordering_lines] 5079 5080 latitudes = take(latitudes_urs, permutation) 5081 longitudes = take(longitudes_urs, permutation) 5059 permutation = ensure_numeric([int(line.split(',')[0]) for line in ordering_lines]) 5060 else: 5061 # Use empty array to signify 'all' points 5062 # We could have used 'None' but it got too hard in the C-code ;-) 5063 permutation = ensure_numeric([], Float) 5064 5065 5066 # Read MUX2 files 5067 if (verbose): print 'reading mux2 file' 5068 mux={} 5069 for i, quantity in enumerate(quantities): 5070 5071 # For each quantity read the associated list of source mux2 file with 5072 # extention associated with that quantity 5073 times,\ 5074 latitudes,\ 5075 longitudes,\ 5076 elevation,\ 5077 mux[quantity],\ 5078 starttime = read_mux2_py(files_in[i], weights, permutation, verbose=verbose) 5079 5080 # Check that all quantities have consistent time and space information 5081 if quantity!=quantities[0]: 5082 msg='%s, %s and %s have inconsistent gauge data'\ 5083 %(files_in[0], files_in[1], files_in[2]) 5084 assert allclose(times, times_old), msg 5085 assert allclose(latitudes, latitudes_old), msg 5086 assert allclose(longitudes, longitudes_old), msg 5087 assert allclose(elevation, elevation_old), msg 5088 assert allclose(starttime, starttime_old), msg 5089 times_old=times 5090 latitudes_old=latitudes 5091 longitudes_old=longitudes 5092 elevation_old=elevation 5093 starttime_old=starttime 5094 5095 5096 5097 5082 5098 5083 5099 # Self check - can be removed to improve speed 5084 ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines] 5085 ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines] 5086 5087 msg = 'Longitudes specified in ordering file do not match those found in mux files' 5088 msg += 'I got %s instead of %s (only beginning shown)'\ 5089 %(str(longitudes[:10]) + '...', 5090 str(ref_longitudes[:10]) + '...') 5091 assert allclose(longitudes, ref_longitudes), msg 5092 5093 msg = 'Latitudes specified in ordering file do not match those found in mux files. ' 5094 msg += 'I got %s instead of %s (only beginning shown)'\ 5095 %(str(latitudes[:10]) + '...', 5096 str(ref_latitudes[:10]) + '...') 5097 assert allclose(latitudes, ref_latitudes), msg 5098 5099 else: 5100 latitudes=latitudes_urs 5101 longitudes=longitudes_urs 5102 permutation = range(latitudes.shape[0]) 5103 5104 5100 5101 #ref_longitudes = [float(line.split(',')[1]) for line in ordering_lines] 5102 #ref_latitudes = [float(line.split(',')[2]) for line in ordering_lines] 5103 # 5104 #msg = 'Longitudes specified in ordering file do not match those found in mux files' 5105 #msg += 'I got %s instead of %s (only beginning shown)'\ 5106 # %(str(longitudes[:10]) + '...', 5107 # str(ref_longitudes[:10]) + '...') 5108 #assert allclose(longitudes, ref_longitudes), msg 5109 # 5110 #msg = 'Latitudes specified in ordering file do not match those found in mux files. ' 5111 #msg += 'I got %s instead of %s (only beginning shown)'\ 5112 # %(str(latitudes[:10]) + '...', 5113 # str(ref_latitudes[:10]) + '...') 5114 #assert allclose(latitudes, ref_latitudes), msg 5115 5116 5105 5117 5106 5118 # Store timeseries in NetCDF STS file … … 5108 5120 assert len(latitudes>0),msg 5109 5121 5110 number_of_points = latitudes.shape[0] 5111 number_of_times = times.shape[0] 5112 number_of_latitudes = latitudes.shape[0] 5113 number_of_longitudes = longitudes.shape[0] 5122 number_of_points = latitudes.shape[0] # Number of stations retrieved 5123 number_of_times = times.shape[0] # Number of timesteps 5124 number_of_latitudes = latitudes.shape[0] # Number latitudes 5125 number_of_longitudes = longitudes.shape[0] # Number longitudes 5114 5126 5115 5127 # NetCDF file definition … … 5120 5132 5121 5133 # Create new file 5122 #starttime = times[0]5123 5134 sts = Write_sts() 5124 5135 sts.store_header(outfile, … … 5131 5142 # Store 5132 5143 from anuga.coordinate_transforms.redfearn import redfearn 5133 x = zeros(number_of_points, Float) # Easting5134 y = zeros(number_of_points, Float) # Northing5144 x = zeros(number_of_points, Float) # Easting 5145 y = zeros(number_of_points, Float) # Northing 5135 5146 5136 5147 # Check zone boundaries … … 5165 5176 if verbose: print 'Converting quantities' 5166 5177 for j in range(len(times)): 5167 for i , index in enumerate(permutation):5168 w = zscale*mux['HA'][i ndex,j] + mean_stage5169 h=w-elevation[i ndex]5178 for i in range(number_of_points): 5179 w = zscale*mux['HA'][i,j] + mean_stage 5180 h=w-elevation[i] 5170 5181 stage[j,i] = w 5171 5182 5172 xmomentum[j,i] = mux['UA'][i ndex,j]*h5173 ymomentum[j,i] = mux['VA'][i ndex,j]*h5183 xmomentum[j,i] = mux['UA'][i,j]*h 5184 ymomentum[j,i] = mux['VA'][i,j]*h 5174 5185 5175 5186 outfile.close()
Note: See TracChangeset
for help on using the changeset viewer.