Changeset 3720
- Timestamp:
- Oct 10, 2006, 11:02:18 AM (18 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
r3694 r3720 60 60 import csv 61 61 import os 62 from struct import unpack 63 import array as p_array 62 64 63 65 from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \ 64 searchsorted, zeros, allclose, around 66 searchsorted, zeros, allclose, around, reshape 67 from Scientific.IO.NetCDF import NetCDFFile 65 68 66 69 … … 68 71 from anuga.geospatial_data.geospatial_data import Geospatial_data 69 72 from anuga.config import minimum_storable_height as default_minimum_storable_height 73 from anuga.utilities.numerical_tools import ensure_numeric 74 from anuga.caching.caching import myhash 75 from anuga.utilities.anuga_exceptions import ANUGAError 70 76 71 77 def make_filename(s): … … 4043 4049 return 4044 4050 4051 #### URS 2 SWW ### 4052 4053 lon_name = 'LON' 4054 lat_name = 'LAT' 4055 time_name = 'TIME' 4056 precision = Float # So if we want to change the precision its done here 4057 class Write_nc: 4058 """ 4059 Write an nc file. 4060 4061 Note, this should be checked to meet cdc netcdf conventions for gridded 4062 data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml 4063 4064 """ 4065 def __init__(self, 4066 quantity_name, 4067 file_name, 4068 time_step_count, 4069 time_step, 4070 lon, 4071 lat): 4072 """ 4073 time_step_count is the number of time steps. 4074 time_step is the time step size 4075 4076 pre-condition: quantity_name must be 'HA' 'UA'or 'VA'. 4077 """ 4078 self.quantity_name = quantity_name 4079 quantity_units= {'HA':'METERS', 4080 'UA':'METERS/SECOND', 4081 'VA':'METERS/SECOND'} 4082 4083 #self.file_name = file_name 4084 self.time_step_count = time_step_count 4085 self.time_step = time_step 4086 4087 # NetCDF file definition 4088 self.outfile = NetCDFFile(file_name, 'w') 4089 outfile = self.outfile 4090 4091 #Create new file 4092 nc_lon_lat_header(outfile, lon, lat) 4093 4094 # TIME 4095 outfile.createDimension(time_name, None) 4096 outfile.createVariable(time_name, precision, (time_name,)) 4097 4098 #QUANTITY 4099 outfile.createVariable(self.quantity_name, precision, 4100 (time_name, lat_name, lon_name)) 4101 outfile.variables[self.quantity_name].missing_value=-1.e+034 4102 outfile.variables[self.quantity_name].units= \ 4103 quantity_units[self.quantity_name] 4104 outfile.variables[lon_name][:]= ensure_numeric(lon) 4105 outfile.variables[lat_name][:]= ensure_numeric(lat) 4106 4107 #Assume no one will be wanting to read this, while we are writing 4108 #outfile.close() 4109 4110 def store_timestep(self, quantity_slice): 4111 """ 4112 quantity_slice is the data to be stored at this time step 4113 """ 4114 4115 outfile = self.outfile 4116 4117 # Get the variables 4118 time = outfile.variables[time_name] 4119 quantity = outfile.variables[self.quantity_name] 4120 4121 i = len(time) 4122 4123 #Store time 4124 time[i] = i*self.time_step #self.domain.time 4125 quantity[i,:] = quantity_slice*100 # To convert from m to cm 4126 # And m/s to cm/sec 4127 4128 def close(self): 4129 self.outfile.close() 4130 4131 def urs2sww(basename_in='o', basename_out=None, verbose=False, 4132 remove_nc_files=True, 4133 minlat=None, maxlat=None, 4134 minlon= None, maxlon=None, 4135 mint=None, maxt=None, 4136 mean_stage=0, 4137 zscale=1, 4138 fail_on_NaN=True, 4139 NaN_filler=0, 4140 elevation=None): 4141 """ 4142 Convert URS C binary format for wave propagation to 4143 sww format native to abstract_2d_finite_volumes. 4144 4145 Specify only basename_in and read files of the form 4146 basefilename-z-mux, basefilename-e-mux and basefilename-n-mux containing 4147 relative height, x-velocity and y-velocity, respectively. 4148 4149 Also convert latitude and longitude to UTM. All coordinates are 4150 assumed to be given in the GDA94 datum. The latitude and longitude 4151 information is for a grid. 4152 4153 min's and max's: If omitted - full extend is used. 4154 To include a value min may equal it, while max must exceed it. 4155 Lat and lon are assuemd to be in decimal degrees 4156 4157 URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE 4158 which means that latitude is the fastest 4159 varying dimension (row major order, so to speak) 4160 4161 In URS C binary the latitudes and longitudes are in assending order. 4162 """ 4163 4164 if basename_out == None: 4165 basename_out = basename_in 4166 files_out = urs2nc(basename_in, basename_out) 4167 ferret2sww(basename_out, 4168 minlat=minlat, 4169 maxlat=maxlat, 4170 minlon=minlon, 4171 mint=mint, 4172 maxt=maxt, 4173 mean_stage=mean_stage, 4174 zscale=zscale, 4175 fail_on_NaN=fail_on_NaN, 4176 NaN_filler=NaN_filler, 4177 inverted_bathymetry=True, 4178 verbose=verbose) 4179 #print "files_out",files_out 4180 if remove_nc_files: 4181 for file_out in files_out: 4182 os.remove(file_out) 4183 4184 def urs2nc(basename_in = 'o', basename_out = 'urs'): 4185 files_in = [basename_in+'-z-mux', 4186 basename_in+'-e-mux', 4187 basename_in+'-n-mux'] 4188 files_out = [basename_out+'_ha.nc', 4189 basename_out+'_ua.nc', 4190 basename_out+'_va.nc'] 4191 quantities = ['HA','UA','VA'] 4192 4193 hashed_elevation = None 4194 for file_in, file_out, quantity in map(None, files_in, 4195 files_out, 4196 quantities): 4197 lonlatdep, lon, lat = _binary_c2nc(file_in, 4198 file_out, 4199 quantity) 4200 #print "lon",lon 4201 #print "lat",lat 4202 if hashed_elevation == None: 4203 elevation_file = basename_out+'_e.nc' 4204 write_elevation_sww(elevation_file, 4205 lon, 4206 lat, 4207 lonlatdep[:,2]) 4208 hashed_elevation = myhash(lonlatdep) 4209 else: 4210 msg = "The elevation information in the mux files is inconsistent" 4211 assert hashed_elevation == myhash(lonlatdep), msg 4212 files_out.append(elevation_file) 4213 return files_out 4214 4215 def _binary_c2nc(file_in, file_out, quantity): 4216 """ 4217 4218 return the depth info, so it can be written to a file 4219 """ 4220 columns = 3 # long, lat , depth 4221 f = open(file_in, 'rb') 4222 4223 # Number of points/stations 4224 (points_num,)= unpack('i',f.read(4)) 4225 4226 # nt, int - Number of time steps 4227 (time_step_count,)= unpack('i',f.read(4)) 4228 4229 #dt, float - time step, seconds 4230 (time_step,) = unpack('f', f.read(4)) 4231 4232 msg = "Bad data in the mux file." 4233 if points_num < 0: 4234 raise ANUGAError, msg 4235 if time_step_count < 0: 4236 raise ANUGAError, msg 4237 if time_step < 0: 4238 raise ANUGAError, msg 4239 4240 lonlatdep = p_array.array('f') 4241 lonlatdep.read(f, columns * points_num) 4242 lonlatdep = array(lonlatdep, typecode=Float) 4243 lonlatdep = reshape(lonlatdep, ( points_num, columns)) 4244 4245 lon, lat = lon_lat2grid(lonlatdep) 4246 nc_file = Write_nc(quantity, 4247 file_out, 4248 time_step_count, 4249 time_step, 4250 lon, 4251 lat) 4252 4253 for i in range(time_step_count): 4254 #Read in a time slice 4255 hz_p_array = p_array.array('f') 4256 hz_p_array.read(f, points_num) 4257 hz_p = array(hz_p_array, typecode=Float) 4258 hz_p = reshape(hz_p, ( len(lat), len(lon))) 4259 4260 nc_file.store_timestep(hz_p) 4261 4262 nc_file.close() 4263 4264 return lonlatdep, lon, lat 4265 4266 4267 def write_elevation_sww(file_out, lon, lat, depth_vector): 4268 4269 # NetCDF file definition 4270 outfile = NetCDFFile(file_out, 'w') 4271 4272 #Create new file 4273 nc_lon_lat_header(outfile, lon, lat) 4274 4275 # ELEVATION 4276 zname = 'ELEVATION' 4277 outfile.createVariable(zname, precision, (lat_name, lon_name)) 4278 outfile.variables[zname].units='CENTIMETERS' 4279 outfile.variables[zname].missing_value=-1.e+034 4280 4281 outfile.variables[lon_name][:]= ensure_numeric(lon) 4282 outfile.variables[lat_name][:]= ensure_numeric(lat) 4283 4284 depth = reshape(depth_vector, ( len(lat), len(lon))) 4285 #print "depth",depth 4286 outfile.variables[zname][:]= depth 4287 4288 outfile.close() 4289 4290 def nc_lon_lat_header(outfile, lon, lat): 4291 4292 outfile.institution = 'Geoscience Australia' 4293 outfile.description = 'Converted from URS binary C' 4294 4295 # Longitude 4296 outfile.createDimension(lon_name, len(lon)) 4297 outfile.createVariable(lon_name, precision, (lon_name,)) 4298 outfile.variables[lon_name].point_spacing='uneven' 4299 outfile.variables[lon_name].units='degrees_east' 4300 outfile.variables[lon_name].assignValue(lon) 4301 4302 4303 # Latitude 4304 outfile.createDimension(lat_name, len(lat)) 4305 outfile.createVariable(lat_name, precision, (lat_name,)) 4306 outfile.variables[lat_name].point_spacing='uneven' 4307 outfile.variables[lat_name].units='degrees_north' 4308 outfile.variables[lat_name].assignValue(lat) 4309 4310 4311 4312 def lon_lat2grid(long_lat_dep): 4313 """ 4314 given a list of points that are assumed to be an a grid, 4315 return the long's and lat's of the grid. 4316 long_lat_dep is an array where each row is a position. 4317 The first column is longitudes. 4318 The second column is latitudes. 4319 """ 4320 LONG = 0 4321 LAT = 1 4322 points_num = len(long_lat_dep) 4323 lat = [long_lat_dep[0][LAT]] 4324 this_rows_long = long_lat_dep[0][LONG] 4325 i = 1 # Index of long_lat_dep 4326 4327 #Working out the lat's 4328 while long_lat_dep[i][LONG] == this_rows_long and i < points_num: 4329 lat.append(long_lat_dep[i][LAT]) 4330 i += 1 4331 4332 lats_per_long = i 4333 long = [long_lat_dep[0][LONG]] 4334 #Working out the longs 4335 while i < points_num: 4336 msg = 'Data is not gridded. It must be for this operation' 4337 assert long_lat_dep[i][LAT] == lat[i%lats_per_long], msg 4338 if i%lats_per_long == 0: 4339 long.append(long_lat_dep[i][LONG]) 4340 i += 1 4341 4342 msg = 'Our of range latitudes/longitudes' 4343 for l in lat:assert -90 < l < 90 , msg 4344 for l in long:assert -180 < l < 180 , msg 4345 4346 return long, lat 4347 4348 #### END URS 2 SWW ### 4045 4349 #------------------------------------------------------------- 4046 4350 if __name__ == "__main__": -
anuga_core/source/anuga/shallow_water/test_data_manager.py
r3694 r3720 9 9 import os 10 10 from Scientific.IO.NetCDF import NetCDFFile 11 from struct import pack 11 12 12 13 from anuga.shallow_water.data_manager import * … … 16 17 # This is needed to run the tests of local functions 17 18 import data_manager 18 19 from anuga.coordinate_transforms.redfearn import redfearn 19 20 from anuga.coordinate_transforms.geo_reference import Geo_reference 20 21 … … 71 72 self.F = bed 72 73 73 74 75 76 74 #Write A testfile (not realistic. Values aren't realistic) 77 78 75 self.test_MOST_file = 'most_small' 79 76 … … 4984 4981 4985 4982 os.remove(csv_file) 4983 4984 #### TESTS FOR URS 2 SWW ### 4985 4986 def create_mux(self, points_num=None): 4987 # write all the mux stuff. 4988 time_step_count = 3 4989 time_step = 0.5 4990 4991 longitudes = [150.66667, 150.83334, 151., 151.16667] 4992 latitudes = [-34.5, -34.33333, -34.16667, -34] 4993 4994 if points_num == None: 4995 points_num = len(longitudes) * len(latitudes) 4996 4997 lonlatdeps = [] 4998 quantities = ['HA','UA','VA'] 4999 mux_names = ['-z-mux','-e-mux','-n-mux'] 5000 quantities_init = [[],[],[]] 5001 for i,lon in enumerate(longitudes): 5002 for j,lat in enumerate(latitudes): 5003 _ , e, n = redfearn(lat, lon) 5004 lonlatdeps.append([lon, lat, n]) 5005 quantities_init[0].append(e) # HA 5006 quantities_init[1].append(n ) # UA 5007 quantities_init[2].append(e) # VA 5008 #print "lonlatdeps",lonlatdeps 5009 base_name = str(id(self)) 5010 files = [] 5011 for i,q in enumerate(quantities): 5012 quantities_init[i] = ensure_numeric(quantities_init[i]) 5013 #print "HA_init", HA_init 5014 q_time = zeros((time_step_count, points_num), Float) 5015 for time in range(time_step_count): 5016 q_time[time,:] = quantities_init[i] #* time * 4 5017 5018 #Write C files 5019 columns = 3 # long, lat , depth 5020 file = base_name + mux_names[i] 5021 f = open(file, 'wb') 5022 files.append(file) 5023 f.write(pack('i',points_num)) 5024 f.write(pack('i',time_step_count)) 5025 f.write(pack('f',time_step)) 5026 5027 #write lat/long info 5028 for lonlatdep in lonlatdeps: 5029 for float in lonlatdep: 5030 f.write(pack('f',float)) 5031 5032 # Write quantity info 5033 for time in range(time_step_count): 5034 for i in range(points_num): 5035 f.write(pack('f',q_time[time,i])) 5036 f.close() 5037 return base_name, files 5038 5039 5040 def delete_mux(self, files): 5041 for file in files: 5042 os.remove(file) 5043 5044 def test_urs2sww_test_fail(self): 5045 points_num = -100 5046 time_step_count = 45 5047 time_step = -7 5048 base_name = str(id(self)) 5049 files = [] 5050 quantities = ['HA','UA','VA'] 5051 mux_names = ['-z-mux','-e-mux','-n-mux'] 5052 for i,q in enumerate(quantities): 5053 #Write C files 5054 columns = 3 # long, lat , depth 5055 file = base_name + mux_names[i] 5056 f = open(file, 'wb') 5057 files.append(file) 5058 f.write(pack('i',points_num)) 5059 f.write(pack('i',time_step_count)) 5060 f.write(pack('f',time_step)) 5061 5062 f.close() 5063 tide = 1 5064 try: 5065 urs2sww(base_name, remove_nc_files=True, mean_stage=tide) 5066 except ANUGAError: 5067 pass 5068 else: 5069 self.delete_mux(files) 5070 msg = 'Should have raised exception' 5071 raise msg 5072 sww_file = base_name + '.sww' 5073 5074 self.delete_mux(files) 5075 5076 def test_urs2sww(self): 5077 tide = 1 5078 base_name, files = self.create_mux() 5079 urs2sww(base_name, remove_nc_files=True, mean_stage=tide) 5080 sww_file = base_name + '.sww' 5081 5082 #Let's interigate the sww file 5083 fid = NetCDFFile(sww_file) 5084 5085 x = fid.variables['x'][:] 5086 y = fid.variables['y'][:] 5087 geo_reference = Geo_reference(NetCDFObject=fid) 5088 5089 #Check that first coordinate is correctly represented 5090 #Work out the UTM coordinates for first point 5091 zone, e, n = redfearn(-34.5, 150.66667) 5092 5093 assert allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n]) 5094 5095 #Check first value 5096 stage = fid.variables['stage'][:] 5097 xmomentum = fid.variables['xmomentum'][:] 5098 ymomentum = fid.variables['ymomentum'][:] 5099 5100 #print ymomentum 5101 assert allclose(stage[0,0], e +tide) #Meters 5102 5103 #Check the momentums - ua 5104 #momentum = velocity*(stage-elevation) 5105 #momentum = velocity*(stage+elevation) 5106 # -(-elevation) since elevation is inverted in mux files 5107 # = n*(e+tide+n) based on how I'm writing these files 5108 answer = n*(e+tide+n) 5109 actual = xmomentum[0,0] 5110 #print "answer",answer 5111 #print "actual",actual 5112 assert allclose(answer, actual) #Meters 5113 5114 fid.close() 5115 5116 5117 self.delete_mux(files) 5118 os.remove(sww_file) 5119 5120 def bad_test_assuming_theres_a_mux_file(self): 5121 # these mux files aren't in the repository, plus they are bad! 5122 base_name = 'o-z-mux' 5123 base_name = 'o-e-mux' 5124 file_name = base_name + '.nc' 5125 lonlatdep_numeric, lon, lat = _binary_c2nc(base_name, file_name, 'HA') 5126 5127 #os.remove(file_name) 5128 5129 def test_lon_lat2grid(self): 5130 lonlatdep = [ 5131 [ 113.06700134 , -26.06669998 , 0. ] , 5132 [ 113.06700134 , -26.33329964 , 0. ] , 5133 [ 113.19999695 , -26.06669998 , 0. ] , 5134 [ 113.19999695 , -26.33329964 , 0. ] ] 5135 5136 long, lat = lon_lat2grid(lonlatdep) 5137 5138 for i, result in enumerate(lat): 5139 assert lonlatdep [i][1] == result 5140 assert len(lat) == 2 5141 5142 for i, result in enumerate(long): 5143 assert lonlatdep [i*2][0] == result 5144 assert len(long) == 2 5145 5146 def test_lon_lat2grid_bad(self): 5147 lonlatdep = [ 5148 [ -26.06669998, 113.06700134, 1. ], 5149 [ -26.06669998 , 113.19999695 , 2. ], 5150 [ -26.06669998 , 113.33300018, 3. ], 5151 [ -26.06669998 , 113.43299866 , 4. ], 5152 [ -26.20000076 , 113.06700134, 5. ], 5153 [ -26.20000076 , 113.19999695 , 6. ], 5154 [ -26.20000076 , 113.33300018 , 7. ], 5155 [ -26.20000076 , 113.43299866 , 8. ], 5156 [ -26.33329964 , 113.06700134, 9. ], 5157 [ -26.33329964 , 113.19999695 , 10. ], 5158 [ -26.33329964 , 113.33300018 , 11. ], 5159 [ -26.33329964 , 113.43299866 , 12. ], 5160 [ -26.43330002 , 113.06700134 , 13 ], 5161 [ -26.43330002 , 113.19999695 , 14. ], 5162 [ -26.43330002 , 113.33300018, 15. ], 5163 [ -26.43330002 , 113.43299866, 16. ]] 5164 try: 5165 long, lat = lon_lat2grid(lonlatdep) 5166 except AssertionError: 5167 pass 5168 else: 5169 msg = 'Should have raised exception' 5170 raise msg 5171 5172 def test_lon_lat2gridII(self): 5173 lonlatdep = [ 5174 [ 113.06700134 , -26.06669998 , 0. ] , 5175 [ 113.06700134 , -26.33329964 , 0. ] , 5176 [ 113.19999695 , -26.06669998 , 0. ] , 5177 [ 113.19999695 , -26.344329964 , 0. ] ] 5178 try: 5179 long, lat = lon_lat2grid(lonlatdep) 5180 except AssertionError: 5181 pass 5182 else: 5183 msg = 'Should have raised exception' 5184 raise msg 5185 5186 def trial_loading(self): 5187 basename_in = 'karratha' 5188 basename_out = basename_in 5189 urs2sww(basename_in, basename_out) 5190 #### END TESTS FOR URS 2 SWW ### 5191 4986 5192 4987 5193 #-------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.