1 | """ Conversion routines. |
---|
2 | ANUGA needs to deal with many different file formats, and this |
---|
3 | module provides routines for easily converting between them. |
---|
4 | |
---|
5 | These routines are necessarily high level, sitting above the various |
---|
6 | ANUGA modules. They take a file as input, and output a file. |
---|
7 | """ |
---|
8 | |
---|
9 | #non ANUGA imports |
---|
10 | from Scientific.IO.NetCDF import NetCDFFile |
---|
11 | import numpy as num |
---|
12 | import os.path |
---|
13 | import os |
---|
14 | |
---|
15 | #ANUGA imports |
---|
16 | from anuga.coordinate_transforms.geo_reference import Geo_reference, \ |
---|
17 | write_NetCDF_georeference, ensure_geo_reference |
---|
18 | from anuga.abstract_2d_finite_volumes.pmesh2domain import \ |
---|
19 | pmesh_to_domain_instance |
---|
20 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
21 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \ |
---|
22 | netcdf_float |
---|
23 | |
---|
24 | |
---|
25 | #shallow water imports |
---|
26 | from anuga.shallow_water.sww_file import Read_sww, Write_sww |
---|
27 | from anuga.shallow_water.shallow_water_domain import Domain |
---|
28 | from anuga.shallow_water.shallow_water_domain import Domain |
---|
29 | |
---|
30 | |
---|
31 | def sww2obj(basefilename, size): |
---|
32 | """ Convert netcdf based data output to obj |
---|
33 | |
---|
34 | Convert SWW data to OBJ data. |
---|
35 | basefilename Stem of filename, needs size and extension added. |
---|
36 | size The number of lines to write. |
---|
37 | """ |
---|
38 | |
---|
39 | from Scientific.IO.NetCDF import NetCDFFile |
---|
40 | |
---|
41 | # Get NetCDF |
---|
42 | FN = create_filename('.', basefilename, 'sww', size) |
---|
43 | log.critical('Reading from %s' % FN) |
---|
44 | fid = NetCDFFile(FN, netcdf_mode_r) #Open existing file for read |
---|
45 | |
---|
46 | # Get the variables |
---|
47 | x = fid.variables['x'] |
---|
48 | y = fid.variables['y'] |
---|
49 | z = fid.variables['elevation'] |
---|
50 | time = fid.variables['time'] |
---|
51 | stage = fid.variables['stage'] |
---|
52 | |
---|
53 | M = size #Number of lines |
---|
54 | xx = num.zeros((M,3), num.float) |
---|
55 | yy = num.zeros((M,3), num.float) |
---|
56 | zz = num.zeros((M,3), num.float) |
---|
57 | |
---|
58 | for i in range(M): |
---|
59 | for j in range(3): |
---|
60 | xx[i,j] = x[i+j*M] |
---|
61 | yy[i,j] = y[i+j*M] |
---|
62 | zz[i,j] = z[i+j*M] |
---|
63 | |
---|
64 | # Write obj for bathymetry |
---|
65 | FN = create_filename('.', basefilename, 'obj', size) |
---|
66 | write_obj(FN,xx,yy,zz) |
---|
67 | |
---|
68 | # Now read all the data with variable information, combine with |
---|
69 | # x,y info and store as obj |
---|
70 | for k in range(len(time)): |
---|
71 | t = time[k] |
---|
72 | log.critical('Processing timestep %f' % t) |
---|
73 | |
---|
74 | for i in range(M): |
---|
75 | for j in range(3): |
---|
76 | zz[i,j] = stage[k,i+j*M] |
---|
77 | |
---|
78 | #Write obj for variable data |
---|
79 | #FN = create_filename(basefilename, 'obj', size, time=t) |
---|
80 | FN = create_filename('.', basefilename[:5], 'obj', size, time=t) |
---|
81 | write_obj(FN, xx, yy, zz) |
---|
82 | |
---|
83 | |
---|
84 | ## |
---|
85 | # @brief |
---|
86 | # @param basefilename Stem of filename, needs size and extension added. |
---|
87 | def dat2obj(basefilename): |
---|
88 | """Convert line based data output to obj |
---|
89 | FIXME: Obsolete? |
---|
90 | """ |
---|
91 | |
---|
92 | import glob, os |
---|
93 | from anuga.config import data_dir |
---|
94 | |
---|
95 | # Get bathymetry and x,y's |
---|
96 | lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines() |
---|
97 | |
---|
98 | M = len(lines) #Number of lines |
---|
99 | x = num.zeros((M,3), num.float) |
---|
100 | y = num.zeros((M,3), num.float) |
---|
101 | z = num.zeros((M,3), num.float) |
---|
102 | |
---|
103 | for i, line in enumerate(lines): |
---|
104 | tokens = line.split() |
---|
105 | values = map(float, tokens) |
---|
106 | |
---|
107 | for j in range(3): |
---|
108 | x[i,j] = values[j*3] |
---|
109 | y[i,j] = values[j*3+1] |
---|
110 | z[i,j] = values[j*3+2] |
---|
111 | |
---|
112 | # Write obj for bathymetry |
---|
113 | write_obj(data_dir + os.sep + basefilename + '_geometry', x, y, z) |
---|
114 | |
---|
115 | # Now read all the data files with variable information, combine with |
---|
116 | # x,y info and store as obj. |
---|
117 | |
---|
118 | files = glob.glob(data_dir + os.sep + basefilename + '*.dat') |
---|
119 | for filename in files: |
---|
120 | log.critical('Processing %s' % filename) |
---|
121 | |
---|
122 | lines = open(data_dir + os.sep + filename, 'r').readlines() |
---|
123 | assert len(lines) == M |
---|
124 | root, ext = os.path.splitext(filename) |
---|
125 | |
---|
126 | # Get time from filename |
---|
127 | i0 = filename.find('_time=') |
---|
128 | if i0 == -1: |
---|
129 | #Skip bathymetry file |
---|
130 | continue |
---|
131 | |
---|
132 | i0 += 6 #Position where time starts |
---|
133 | i1 = filename.find('.dat') |
---|
134 | |
---|
135 | if i1 > i0: |
---|
136 | t = float(filename[i0:i1]) |
---|
137 | else: |
---|
138 | raise DataTimeError, 'Hmmmm' |
---|
139 | |
---|
140 | for i, line in enumerate(lines): |
---|
141 | tokens = line.split() |
---|
142 | values = map(float,tokens) |
---|
143 | |
---|
144 | for j in range(3): |
---|
145 | z[i,j] = values[j] |
---|
146 | |
---|
147 | # Write obj for variable data |
---|
148 | write_obj(data_dir + os.sep + basefilename + '_time=%.4f' % t, x, y, z) |
---|
149 | |
---|
150 | ## |
---|
151 | # @brief Convert time-series text file to TMS file. |
---|
152 | # @param filename |
---|
153 | # @param quantity_names |
---|
154 | # @param time_as_seconds |
---|
155 | def timefile2netcdf(filename, quantity_names=None, time_as_seconds=False): |
---|
156 | """Template for converting typical text files with time series to |
---|
157 | NetCDF tms file. |
---|
158 | |
---|
159 | The file format is assumed to be either two fields separated by a comma: |
---|
160 | |
---|
161 | time [DD/MM/YY hh:mm:ss], value0 value1 value2 ... |
---|
162 | |
---|
163 | E.g |
---|
164 | |
---|
165 | 31/08/04 00:00:00, 1.328223 0 0 |
---|
166 | 31/08/04 00:15:00, 1.292912 0 0 |
---|
167 | |
---|
168 | or time (seconds), value0 value1 value2 ... |
---|
169 | |
---|
170 | 0.0, 1.328223 0 0 |
---|
171 | 0.1, 1.292912 0 0 |
---|
172 | |
---|
173 | will provide a time dependent function f(t) with three attributes |
---|
174 | |
---|
175 | filename is assumed to be the rootname with extenisons .txt and .sww |
---|
176 | """ |
---|
177 | |
---|
178 | import time, calendar |
---|
179 | from anuga.config import time_format |
---|
180 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
181 | |
---|
182 | file_text = filename + '.txt' |
---|
183 | fid = open(file_text) |
---|
184 | line = fid.readline() |
---|
185 | fid.close() |
---|
186 | |
---|
187 | fields = line.split(',') |
---|
188 | msg = "File %s must have the format 'datetime, value0 value1 value2 ...'" \ |
---|
189 | % file_text |
---|
190 | assert len(fields) == 2, msg |
---|
191 | |
---|
192 | if not time_as_seconds: |
---|
193 | try: |
---|
194 | starttime = calendar.timegm(time.strptime(fields[0], time_format)) |
---|
195 | except ValueError: |
---|
196 | msg = 'First field in file %s must be' % file_text |
---|
197 | msg += ' date-time with format %s.\n' % time_format |
---|
198 | msg += 'I got %s instead.' % fields[0] |
---|
199 | raise DataTimeError, msg |
---|
200 | else: |
---|
201 | try: |
---|
202 | starttime = float(fields[0]) |
---|
203 | except Error: |
---|
204 | msg = "Bad time format" |
---|
205 | raise DataTimeError, msg |
---|
206 | |
---|
207 | # Split values |
---|
208 | values = [] |
---|
209 | for value in fields[1].split(): |
---|
210 | values.append(float(value)) |
---|
211 | |
---|
212 | q = ensure_numeric(values) |
---|
213 | |
---|
214 | msg = 'ERROR: File must contain at least one independent value' |
---|
215 | assert len(q.shape) == 1, msg |
---|
216 | |
---|
217 | # Read times proper |
---|
218 | from anuga.config import time_format |
---|
219 | import time, calendar |
---|
220 | |
---|
221 | fid = open(file_text) |
---|
222 | lines = fid.readlines() |
---|
223 | fid.close() |
---|
224 | |
---|
225 | N = len(lines) |
---|
226 | d = len(q) |
---|
227 | |
---|
228 | T = num.zeros(N, num.float) # Time |
---|
229 | Q = num.zeros((N, d), num.float) # Values |
---|
230 | |
---|
231 | for i, line in enumerate(lines): |
---|
232 | fields = line.split(',') |
---|
233 | if not time_as_seconds: |
---|
234 | realtime = calendar.timegm(time.strptime(fields[0], time_format)) |
---|
235 | else: |
---|
236 | realtime = float(fields[0]) |
---|
237 | T[i] = realtime - starttime |
---|
238 | |
---|
239 | for j, value in enumerate(fields[1].split()): |
---|
240 | Q[i, j] = float(value) |
---|
241 | |
---|
242 | msg = 'File %s must list time as a monotonuosly ' % filename |
---|
243 | msg += 'increasing sequence' |
---|
244 | assert num.alltrue(T[1:] - T[:-1] > 0), msg |
---|
245 | |
---|
246 | #Create NetCDF file |
---|
247 | from Scientific.IO.NetCDF import NetCDFFile |
---|
248 | |
---|
249 | fid = NetCDFFile(filename + '.tms', netcdf_mode_w) |
---|
250 | |
---|
251 | fid.institution = 'Geoscience Australia' |
---|
252 | fid.description = 'Time series' |
---|
253 | |
---|
254 | #Reference point |
---|
255 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
256 | #FIXME: Use Georef |
---|
257 | fid.starttime = starttime |
---|
258 | |
---|
259 | # dimension definitions |
---|
260 | #fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
261 | #fid.createDimension('number_of_vertices', 3) |
---|
262 | |
---|
263 | fid.createDimension('number_of_timesteps', len(T)) |
---|
264 | |
---|
265 | fid.createVariable('time', netcdf_float, ('number_of_timesteps',)) |
---|
266 | |
---|
267 | fid.variables['time'][:] = T |
---|
268 | |
---|
269 | for i in range(Q.shape[1]): |
---|
270 | try: |
---|
271 | name = quantity_names[i] |
---|
272 | except: |
---|
273 | name = 'Attribute%d' % i |
---|
274 | |
---|
275 | fid.createVariable(name, netcdf_float, ('number_of_timesteps',)) |
---|
276 | fid.variables[name][:] = Q[:,i] |
---|
277 | |
---|
278 | fid.close() |
---|
279 | |
---|
280 | |
---|
281 | |
---|
282 | ## |
---|
283 | # @brief |
---|
284 | # @param filename |
---|
285 | # @param verbose |
---|
286 | def tsh2sww(filename, verbose=False): |
---|
287 | """ |
---|
288 | to check if a tsh/msh file 'looks' good. |
---|
289 | """ |
---|
290 | |
---|
291 | if verbose == True: log.critical('Creating domain from %s' % filename) |
---|
292 | |
---|
293 | domain = pmesh_to_domain_instance(filename, Domain) |
---|
294 | |
---|
295 | if verbose == True: log.critical("Number of triangles = %s" % len(domain)) |
---|
296 | |
---|
297 | domain.smooth = True |
---|
298 | domain.format = 'sww' #Native netcdf visualisation format |
---|
299 | file_path, filename = path.split(filename) |
---|
300 | filename, ext = path.splitext(filename) |
---|
301 | domain.set_name(filename) |
---|
302 | domain.reduction = mean |
---|
303 | |
---|
304 | if verbose == True: log.critical("file_path = %s" % file_path) |
---|
305 | |
---|
306 | if file_path == "": |
---|
307 | file_path = "." |
---|
308 | domain.set_datadir(file_path) |
---|
309 | |
---|
310 | if verbose == True: |
---|
311 | log.critical("Output written to %s%s%s.%s" |
---|
312 | % (domain.get_datadir(), sep, domain.get_name(), |
---|
313 | domain.format)) |
---|
314 | |
---|
315 | sww = SWW_file(domain) |
---|
316 | sww.store_connectivity() |
---|
317 | sww.store_timestep() |
---|
318 | |
---|
319 | |
---|
320 | ## |
---|
321 | # @brief Convert CSIRO ESRI file to an SWW boundary file. |
---|
322 | # @param bath_dir |
---|
323 | # @param elevation_dir |
---|
324 | # @param ucur_dir |
---|
325 | # @param vcur_dir |
---|
326 | # @param sww_file |
---|
327 | # @param minlat |
---|
328 | # @param maxlat |
---|
329 | # @param minlon |
---|
330 | # @param maxlon |
---|
331 | # @param zscale |
---|
332 | # @param mean_stage |
---|
333 | # @param fail_on_NaN |
---|
334 | # @param elevation_NaN_filler |
---|
335 | # @param bath_prefix |
---|
336 | # @param elevation_prefix |
---|
337 | # @param verbose |
---|
338 | # @note Also convert latitude and longitude to UTM. All coordinates are |
---|
339 | # assumed to be given in the GDA94 datum. |
---|
340 | def asc_csiro2sww(bath_dir, |
---|
341 | elevation_dir, |
---|
342 | ucur_dir, |
---|
343 | vcur_dir, |
---|
344 | sww_file, |
---|
345 | minlat=None, maxlat=None, |
---|
346 | minlon=None, maxlon=None, |
---|
347 | zscale=1, |
---|
348 | mean_stage=0, |
---|
349 | fail_on_NaN=True, |
---|
350 | elevation_NaN_filler=0, |
---|
351 | bath_prefix='ba', |
---|
352 | elevation_prefix='el', |
---|
353 | verbose=False): |
---|
354 | """ |
---|
355 | Produce an sww boundary file, from esri ascii data from CSIRO. |
---|
356 | |
---|
357 | Also convert latitude and longitude to UTM. All coordinates are |
---|
358 | assumed to be given in the GDA94 datum. |
---|
359 | |
---|
360 | assume: |
---|
361 | All files are in esri ascii format |
---|
362 | |
---|
363 | 4 types of information |
---|
364 | bathymetry |
---|
365 | elevation |
---|
366 | u velocity |
---|
367 | v velocity |
---|
368 | |
---|
369 | Assumptions |
---|
370 | The metadata of all the files is the same |
---|
371 | Each type is in a seperate directory |
---|
372 | One bath file with extention .000 |
---|
373 | The time period is less than 24hrs and uniform. |
---|
374 | """ |
---|
375 | |
---|
376 | from Scientific.IO.NetCDF import NetCDFFile |
---|
377 | |
---|
378 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
379 | |
---|
380 | # So if we want to change the precision it's done here |
---|
381 | precision = netcdf_float |
---|
382 | |
---|
383 | # go in to the bath dir and load the only file, |
---|
384 | bath_files = os.listdir(bath_dir) |
---|
385 | bath_file = bath_files[0] |
---|
386 | bath_dir_file = bath_dir + os.sep + bath_file |
---|
387 | bath_metadata, bath_grid = _read_asc(bath_dir_file) |
---|
388 | |
---|
389 | #Use the date.time of the bath file as a basis for |
---|
390 | #the start time for other files |
---|
391 | base_start = bath_file[-12:] |
---|
392 | |
---|
393 | #go into the elevation dir and load the 000 file |
---|
394 | elevation_dir_file = elevation_dir + os.sep + elevation_prefix \ |
---|
395 | + base_start |
---|
396 | |
---|
397 | elevation_files = os.listdir(elevation_dir) |
---|
398 | ucur_files = os.listdir(ucur_dir) |
---|
399 | vcur_files = os.listdir(vcur_dir) |
---|
400 | elevation_files.sort() |
---|
401 | |
---|
402 | # the first elevation file should be the |
---|
403 | # file with the same base name as the bath data |
---|
404 | assert elevation_files[0] == 'el' + base_start |
---|
405 | |
---|
406 | number_of_latitudes = bath_grid.shape[0] |
---|
407 | number_of_longitudes = bath_grid.shape[1] |
---|
408 | number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2 |
---|
409 | |
---|
410 | longitudes = [bath_metadata['xllcorner'] + x*bath_metadata['cellsize'] \ |
---|
411 | for x in range(number_of_longitudes)] |
---|
412 | latitudes = [bath_metadata['yllcorner'] + y*bath_metadata['cellsize'] \ |
---|
413 | for y in range(number_of_latitudes)] |
---|
414 | |
---|
415 | # reverse order of lat, so the first lat represents the first grid row |
---|
416 | latitudes.reverse() |
---|
417 | |
---|
418 | kmin, kmax, lmin, lmax = get_min_max_indexes(latitudes[:],longitudes[:], |
---|
419 | minlat=minlat, maxlat=maxlat, |
---|
420 | minlon=minlon, maxlon=maxlon) |
---|
421 | |
---|
422 | bath_grid = bath_grid[kmin:kmax,lmin:lmax] |
---|
423 | latitudes = latitudes[kmin:kmax] |
---|
424 | longitudes = longitudes[lmin:lmax] |
---|
425 | number_of_latitudes = len(latitudes) |
---|
426 | number_of_longitudes = len(longitudes) |
---|
427 | number_of_times = len(os.listdir(elevation_dir)) |
---|
428 | number_of_points = number_of_latitudes * number_of_longitudes |
---|
429 | number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2 |
---|
430 | |
---|
431 | #Work out the times |
---|
432 | if len(elevation_files) > 1: |
---|
433 | # Assume: The time period is less than 24hrs. |
---|
434 | time_period = (int(elevation_files[1][-3:]) \ |
---|
435 | - int(elevation_files[0][-3:])) * 60*60 |
---|
436 | times = [x*time_period for x in range(len(elevation_files))] |
---|
437 | else: |
---|
438 | times = [0.0] |
---|
439 | |
---|
440 | if verbose: |
---|
441 | log.critical('------------------------------------------------') |
---|
442 | log.critical('Statistics:') |
---|
443 | log.critical(' Extent (lat/lon):') |
---|
444 | log.critical(' lat in [%f, %f], len(lat) == %d' |
---|
445 | % (min(latitudes), max(latitudes), len(latitudes))) |
---|
446 | log.critical(' lon in [%f, %f], len(lon) == %d' |
---|
447 | % (min(longitudes), max(longitudes), len(longitudes))) |
---|
448 | log.critical(' t in [%f, %f], len(t) == %d' |
---|
449 | % (min(times), max(times), len(times))) |
---|
450 | |
---|
451 | ######### WRITE THE SWW FILE ############# |
---|
452 | |
---|
453 | # NetCDF file definition |
---|
454 | outfile = NetCDFFile(sww_file, netcdf_mode_w) |
---|
455 | |
---|
456 | #Create new file |
---|
457 | outfile.institution = 'Geoscience Australia' |
---|
458 | outfile.description = 'Converted from XXX' |
---|
459 | |
---|
460 | #For sww compatibility |
---|
461 | outfile.smoothing = 'Yes' |
---|
462 | outfile.order = 1 |
---|
463 | |
---|
464 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
465 | outfile.starttime = starttime = times[0] |
---|
466 | |
---|
467 | # dimension definitions |
---|
468 | outfile.createDimension('number_of_volumes', number_of_volumes) |
---|
469 | outfile.createDimension('number_of_vertices', 3) |
---|
470 | outfile.createDimension('number_of_points', number_of_points) |
---|
471 | outfile.createDimension('number_of_timesteps', number_of_times) |
---|
472 | |
---|
473 | # variable definitions |
---|
474 | outfile.createVariable('x', precision, ('number_of_points',)) |
---|
475 | outfile.createVariable('y', precision, ('number_of_points',)) |
---|
476 | outfile.createVariable('elevation', precision, ('number_of_points',)) |
---|
477 | |
---|
478 | #FIXME: Backwards compatibility |
---|
479 | #outfile.createVariable('z', precision, ('number_of_points',)) |
---|
480 | ################################# |
---|
481 | |
---|
482 | outfile.createVariable('volumes', netcdf_int, ('number_of_volumes', |
---|
483 | 'number_of_vertices')) |
---|
484 | |
---|
485 | outfile.createVariable('time', precision, ('number_of_timesteps',)) |
---|
486 | |
---|
487 | outfile.createVariable('stage', precision, ('number_of_timesteps', |
---|
488 | 'number_of_points')) |
---|
489 | |
---|
490 | outfile.createVariable('xmomentum', precision, ('number_of_timesteps', |
---|
491 | 'number_of_points')) |
---|
492 | |
---|
493 | outfile.createVariable('ymomentum', precision, ('number_of_timesteps', |
---|
494 | 'number_of_points')) |
---|
495 | |
---|
496 | #Store |
---|
497 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
498 | |
---|
499 | x = num.zeros(number_of_points, num.float) #Easting |
---|
500 | y = num.zeros(number_of_points, num.float) #Northing |
---|
501 | |
---|
502 | if verbose: log.critical('Making triangular grid') |
---|
503 | |
---|
504 | #Get zone of 1st point. |
---|
505 | refzone, _, _ = redfearn(latitudes[0], longitudes[0]) |
---|
506 | |
---|
507 | vertices = {} |
---|
508 | i = 0 |
---|
509 | for k, lat in enumerate(latitudes): |
---|
510 | for l, lon in enumerate(longitudes): |
---|
511 | vertices[l,k] = i |
---|
512 | |
---|
513 | zone, easting, northing = redfearn(lat, lon) |
---|
514 | |
---|
515 | #msg = 'Zone boundary crossed at longitude =', lon |
---|
516 | #assert zone == refzone, msg |
---|
517 | #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing) |
---|
518 | x[i] = easting |
---|
519 | y[i] = northing |
---|
520 | i += 1 |
---|
521 | |
---|
522 | #Construct 2 triangles per 'rectangular' element |
---|
523 | volumes = [] |
---|
524 | for l in range(number_of_longitudes-1): #X direction |
---|
525 | for k in range(number_of_latitudes-1): #Y direction |
---|
526 | v1 = vertices[l,k+1] |
---|
527 | v2 = vertices[l,k] |
---|
528 | v3 = vertices[l+1,k+1] |
---|
529 | v4 = vertices[l+1,k] |
---|
530 | |
---|
531 | #Note, this is different to the ferrit2sww code |
---|
532 | #since the order of the lats is reversed. |
---|
533 | volumes.append([v1,v3,v2]) #Upper element |
---|
534 | volumes.append([v4,v2,v3]) #Lower element |
---|
535 | |
---|
536 | volumes = num.array(volumes, num.int) #array default# |
---|
537 | |
---|
538 | geo_ref = Geo_reference(refzone, min(x), min(y)) |
---|
539 | geo_ref.write_NetCDF(outfile) |
---|
540 | |
---|
541 | # This will put the geo ref in the middle |
---|
542 | #geo_ref = Geo_reference(refzone, (max(x)+min(x))/2., (max(x)+min(y))/2.) |
---|
543 | |
---|
544 | if verbose: |
---|
545 | log.critical('------------------------------------------------') |
---|
546 | log.critical('More Statistics:') |
---|
547 | log.critical(' Extent (/lon):') |
---|
548 | log.critical(' x in [%f, %f], len(lat) == %d' |
---|
549 | % (min(x), max(x), len(x))) |
---|
550 | log.critical(' y in [%f, %f], len(lon) == %d' |
---|
551 | % (min(y), max(y), len(y))) |
---|
552 | log.critical('geo_ref: ', geo_ref) |
---|
553 | |
---|
554 | z = num.resize(bath_grid,outfile.variables['elevation'][:].shape) |
---|
555 | outfile.variables['x'][:] = x - geo_ref.get_xllcorner() |
---|
556 | outfile.variables['y'][:] = y - geo_ref.get_yllcorner() |
---|
557 | # FIXME (Ole): Remove once viewer has been recompiled and changed |
---|
558 | # to use elevation instead of z |
---|
559 | #outfile.variables['z'][:] = z |
---|
560 | outfile.variables['elevation'][:] = z |
---|
561 | outfile.variables['volumes'][:] = volumes.astype(num.int32) # On Opteron 64 |
---|
562 | |
---|
563 | stage = outfile.variables['stage'] |
---|
564 | xmomentum = outfile.variables['xmomentum'] |
---|
565 | ymomentum = outfile.variables['ymomentum'] |
---|
566 | |
---|
567 | outfile.variables['time'][:] = times #Store time relative |
---|
568 | |
---|
569 | if verbose: log.critical('Converting quantities') |
---|
570 | |
---|
571 | n = number_of_times |
---|
572 | for j in range(number_of_times): |
---|
573 | # load in files |
---|
574 | elevation_meta, elevation_grid = \ |
---|
575 | _read_asc(elevation_dir + os.sep + elevation_files[j]) |
---|
576 | |
---|
577 | _, u_momentum_grid = _read_asc(ucur_dir + os.sep + ucur_files[j]) |
---|
578 | _, v_momentum_grid = _read_asc(vcur_dir + os.sep + vcur_files[j]) |
---|
579 | |
---|
580 | #cut matrix to desired size |
---|
581 | elevation_grid = elevation_grid[kmin:kmax,lmin:lmax] |
---|
582 | u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax] |
---|
583 | v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax] |
---|
584 | |
---|
585 | # handle missing values |
---|
586 | missing = (elevation_grid == elevation_meta['NODATA_value']) |
---|
587 | if num.sometrue (missing): |
---|
588 | if fail_on_NaN: |
---|
589 | msg = 'File %s contains missing values' \ |
---|
590 | % (elevation_files[j]) |
---|
591 | raise DataMissingValuesError, msg |
---|
592 | else: |
---|
593 | elevation_grid = elevation_grid*(missing==0) \ |
---|
594 | + missing*elevation_NaN_filler |
---|
595 | |
---|
596 | if verbose and j % ((n+10)/10) == 0: log.critical(' Doing %d of %d' |
---|
597 | % (j, n)) |
---|
598 | |
---|
599 | i = 0 |
---|
600 | for k in range(number_of_latitudes): #Y direction |
---|
601 | for l in range(number_of_longitudes): #X direction |
---|
602 | w = zscale*elevation_grid[k,l] + mean_stage |
---|
603 | stage[j,i] = w |
---|
604 | h = w - z[i] |
---|
605 | xmomentum[j,i] = u_momentum_grid[k,l]*h |
---|
606 | ymomentum[j,i] = v_momentum_grid[k,l]*h |
---|
607 | i += 1 |
---|
608 | |
---|
609 | outfile.close() |
---|
610 | |
---|
611 | |
---|
612 | |
---|
613 | def _read_asc(filename, verbose=False): |
---|
614 | """Read esri file from the following ASCII format (.asc) |
---|
615 | |
---|
616 | Example: |
---|
617 | |
---|
618 | ncols 3121 |
---|
619 | nrows 1800 |
---|
620 | xllcorner 722000 |
---|
621 | yllcorner 5893000 |
---|
622 | cellsize 25 |
---|
623 | NODATA_value -9999 |
---|
624 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
625 | |
---|
626 | filename Path to the file to read. |
---|
627 | verbose True if this function is to be verbose. |
---|
628 | """ |
---|
629 | |
---|
630 | datafile = open(filename) |
---|
631 | |
---|
632 | if verbose: log.critical('Reading DEM from %s' % filename) |
---|
633 | |
---|
634 | lines = datafile.readlines() |
---|
635 | datafile.close() |
---|
636 | |
---|
637 | if verbose: log.critical('Got %d lines' % len(lines)) |
---|
638 | |
---|
639 | ncols = int(lines.pop(0).split()[1].strip()) |
---|
640 | nrows = int(lines.pop(0).split()[1].strip()) |
---|
641 | xllcorner = float(lines.pop(0).split()[1].strip()) |
---|
642 | yllcorner = float(lines.pop(0).split()[1].strip()) |
---|
643 | cellsize = float(lines.pop(0).split()[1].strip()) |
---|
644 | NODATA_value = float(lines.pop(0).split()[1].strip()) |
---|
645 | |
---|
646 | assert len(lines) == nrows |
---|
647 | |
---|
648 | #Store data |
---|
649 | grid = [] |
---|
650 | |
---|
651 | n = len(lines) |
---|
652 | for i, line in enumerate(lines): |
---|
653 | cells = line.split() |
---|
654 | assert len(cells) == ncols |
---|
655 | grid.append(num.array([float(x) for x in cells])) |
---|
656 | grid = num.array(grid) |
---|
657 | |
---|
658 | return {'xllcorner':xllcorner, |
---|
659 | 'yllcorner':yllcorner, |
---|
660 | 'cellsize':cellsize, |
---|
661 | 'NODATA_value':NODATA_value}, grid |
---|
662 | |
---|
663 | |
---|
664 | ## |
---|
665 | # @brief Convert URS file to SWW file. |
---|
666 | # @param basename_in Stem of the input filename. |
---|
667 | # @param basename_out Stem of the output filename. |
---|
668 | # @param verbose True if this function is to be verbose. |
---|
669 | # @param remove_nc_files |
---|
670 | # @param minlat Sets extent of area to be used. If not supplied, full extent. |
---|
671 | # @param maxlat Sets extent of area to be used. If not supplied, full extent. |
---|
672 | # @param minlon Sets extent of area to be used. If not supplied, full extent. |
---|
673 | # @param maxlon Sets extent of area to be used. If not supplied, full extent. |
---|
674 | # @param mint |
---|
675 | # @param maxt |
---|
676 | # @param mean_stage |
---|
677 | # @param origin A 3-tuple with geo referenced UTM coordinates |
---|
678 | # @param zscale |
---|
679 | # @param fail_on_NaN |
---|
680 | # @param NaN_filler |
---|
681 | # @param elevation |
---|
682 | # @note Also convert latitude and longitude to UTM. All coordinates are |
---|
683 | # assumed to be given in the GDA94 datum. |
---|
684 | def urs2sww(basename_in='o', basename_out=None, verbose=False, |
---|
685 | remove_nc_files=True, |
---|
686 | minlat=None, maxlat=None, |
---|
687 | minlon=None, maxlon=None, |
---|
688 | mint=None, maxt=None, |
---|
689 | mean_stage=0, |
---|
690 | origin=None, |
---|
691 | zscale=1, |
---|
692 | fail_on_NaN=True, |
---|
693 | NaN_filler=0): |
---|
694 | """Convert a URS file to an SWW file. |
---|
695 | Convert URS C binary format for wave propagation to |
---|
696 | sww format native to abstract_2d_finite_volumes. |
---|
697 | |
---|
698 | Specify only basename_in and read files of the form |
---|
699 | basefilename-z-mux2, basefilename-e-mux2 and |
---|
700 | basefilename-n-mux2 containing relative height, |
---|
701 | x-velocity and y-velocity, respectively. |
---|
702 | |
---|
703 | Also convert latitude and longitude to UTM. All coordinates are |
---|
704 | assumed to be given in the GDA94 datum. The latitude and longitude |
---|
705 | information is for a grid. |
---|
706 | |
---|
707 | min's and max's: If omitted - full extend is used. |
---|
708 | To include a value min may equal it, while max must exceed it. |
---|
709 | Lat and lon are assumed to be in decimal degrees. |
---|
710 | NOTE: minlon is the most east boundary. |
---|
711 | |
---|
712 | origin is a 3-tuple with geo referenced |
---|
713 | UTM coordinates (zone, easting, northing) |
---|
714 | It will be the origin of the sww file. This shouldn't be used, |
---|
715 | since all of anuga should be able to handle an arbitary origin. |
---|
716 | |
---|
717 | URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE |
---|
718 | which means that latitude is the fastest |
---|
719 | varying dimension (row major order, so to speak) |
---|
720 | |
---|
721 | In URS C binary the latitudes and longitudes are in assending order. |
---|
722 | """ |
---|
723 | |
---|
724 | if basename_out == None: |
---|
725 | basename_out = basename_in |
---|
726 | |
---|
727 | files_out = urs2nc(basename_in, basename_out) |
---|
728 | |
---|
729 | ferret2sww(basename_out, |
---|
730 | minlat=minlat, |
---|
731 | maxlat=maxlat, |
---|
732 | minlon=minlon, |
---|
733 | maxlon=maxlon, |
---|
734 | mint=mint, |
---|
735 | maxt=maxt, |
---|
736 | mean_stage=mean_stage, |
---|
737 | origin=origin, |
---|
738 | zscale=zscale, |
---|
739 | fail_on_NaN=fail_on_NaN, |
---|
740 | NaN_filler=NaN_filler, |
---|
741 | inverted_bathymetry=True, |
---|
742 | verbose=verbose) |
---|
743 | |
---|
744 | if remove_nc_files: |
---|
745 | for file_out in files_out: |
---|
746 | os.remove(file_out) |
---|
747 | |
---|
748 | |
---|
749 | ## |
---|
750 | # @brief Convert 3 URS files back to 4 NC files. |
---|
751 | # @param basename_in Stem of the input filenames. |
---|
752 | # @param basename_outStem of the output filenames. |
---|
753 | # @note The name of the urs file names must be: |
---|
754 | # [basename_in]-z-mux |
---|
755 | # [basename_in]-e-mux |
---|
756 | # [basename_in]-n-mux |
---|
757 | def urs2nc(basename_in='o', basename_out='urs'): |
---|
758 | """Convert the 3 urs files to 4 nc files. |
---|
759 | |
---|
760 | The name of the urs file names must be; |
---|
761 | [basename_in]-z-mux |
---|
762 | [basename_in]-e-mux |
---|
763 | [basename_in]-n-mux |
---|
764 | """ |
---|
765 | |
---|
766 | files_in = [basename_in + WAVEHEIGHT_MUX_LABEL, |
---|
767 | basename_in + EAST_VELOCITY_LABEL, |
---|
768 | basename_in + NORTH_VELOCITY_LABEL] |
---|
769 | files_out = [basename_out + '_ha.nc', |
---|
770 | basename_out + '_ua.nc', |
---|
771 | basename_out + '_va.nc'] |
---|
772 | quantities = ['HA', 'UA', 'VA'] |
---|
773 | |
---|
774 | #if os.access(files_in[0]+'.mux', os.F_OK) == 0 : |
---|
775 | for i, file_name in enumerate(files_in): |
---|
776 | if os.access(file_name, os.F_OK) == 0: |
---|
777 | if os.access(file_name + '.mux', os.F_OK) == 0 : |
---|
778 | msg = 'File %s does not exist or is not accessible' % file_name |
---|
779 | raise IOError, msg |
---|
780 | else: |
---|
781 | files_in[i] += '.mux' |
---|
782 | log.critical("file_name %s" % file_name) |
---|
783 | |
---|
784 | hashed_elevation = None |
---|
785 | for file_in, file_out, quantity in map(None, files_in, |
---|
786 | files_out, |
---|
787 | quantities): |
---|
788 | lonlatdep, lon, lat, depth = _binary_c2nc(file_in, |
---|
789 | file_out, |
---|
790 | quantity) |
---|
791 | if hashed_elevation == None: |
---|
792 | elevation_file = basename_out + '_e.nc' |
---|
793 | write_elevation_nc(elevation_file, |
---|
794 | lon, |
---|
795 | lat, |
---|
796 | depth) |
---|
797 | hashed_elevation = myhash(lonlatdep) |
---|
798 | else: |
---|
799 | msg = "The elevation information in the mux files is inconsistent" |
---|
800 | assert hashed_elevation == myhash(lonlatdep), msg |
---|
801 | |
---|
802 | files_out.append(elevation_file) |
---|
803 | |
---|
804 | return files_out |
---|
805 | |
---|
806 | |
---|
807 | ## |
---|
808 | # @brief Convert a quantity URS file to a NetCDF file. |
---|
809 | # @param file_in Path to input URS file. |
---|
810 | # @param file_out Path to the output file. |
---|
811 | # @param quantity Name of the quantity to be written to the output file. |
---|
812 | # @return A tuple (lonlatdep, lon, lat, depth). |
---|
813 | def _binary_c2nc(file_in, file_out, quantity): |
---|
814 | """Reads in a quantity urs file and writes a quantity nc file. |
---|
815 | Additionally, returns the depth and lat, long info, |
---|
816 | so it can be written to a file. |
---|
817 | """ |
---|
818 | |
---|
819 | columns = 3 # long, lat , depth |
---|
820 | mux_file = open(file_in, 'rb') |
---|
821 | |
---|
822 | # Number of points/stations |
---|
823 | (points_num,) = unpack('i', mux_file.read(4)) |
---|
824 | |
---|
825 | # nt, int - Number of time steps |
---|
826 | (time_step_count,) = unpack('i', mux_file.read(4)) |
---|
827 | |
---|
828 | #dt, float - time step, seconds |
---|
829 | (time_step,) = unpack('f', mux_file.read(4)) |
---|
830 | |
---|
831 | msg = "Bad data in the mux file." |
---|
832 | if points_num < 0: |
---|
833 | mux_file.close() |
---|
834 | raise ANUGAError, msg |
---|
835 | if time_step_count < 0: |
---|
836 | mux_file.close() |
---|
837 | raise ANUGAError, msg |
---|
838 | if time_step < 0: |
---|
839 | mux_file.close() |
---|
840 | raise ANUGAError, msg |
---|
841 | |
---|
842 | lonlatdep = p_array.array('f') |
---|
843 | lonlatdep.read(mux_file, columns * points_num) |
---|
844 | lonlatdep = num.array(lonlatdep, dtype=num.float) |
---|
845 | lonlatdep = num.reshape(lonlatdep, (points_num, columns)) |
---|
846 | |
---|
847 | lon, lat, depth = lon_lat2grid(lonlatdep) |
---|
848 | lon_sorted = list(lon) |
---|
849 | lon_sorted.sort() |
---|
850 | |
---|
851 | if not num.alltrue(lon == lon_sorted): |
---|
852 | msg = "Longitudes in mux file are not in ascending order" |
---|
853 | raise IOError, msg |
---|
854 | |
---|
855 | lat_sorted = list(lat) |
---|
856 | lat_sorted.sort() |
---|
857 | |
---|
858 | nc_file = Write_nc(quantity, |
---|
859 | file_out, |
---|
860 | time_step_count, |
---|
861 | time_step, |
---|
862 | lon, |
---|
863 | lat) |
---|
864 | |
---|
865 | for i in range(time_step_count): |
---|
866 | #Read in a time slice from mux file |
---|
867 | hz_p_array = p_array.array('f') |
---|
868 | hz_p_array.read(mux_file, points_num) |
---|
869 | hz_p = num.array(hz_p_array, dtype=num.float) |
---|
870 | hz_p = num.reshape(hz_p, (len(lon), len(lat))) |
---|
871 | hz_p = num.transpose(hz_p) # mux has lat varying fastest, nc has long v.f. |
---|
872 | |
---|
873 | #write time slice to nc file |
---|
874 | nc_file.store_timestep(hz_p) |
---|
875 | |
---|
876 | mux_file.close() |
---|
877 | nc_file.close() |
---|
878 | |
---|
879 | return lonlatdep, lon, lat, depth |
---|
880 | |
---|
881 | |
---|
882 | |
---|
883 | ## |
---|
884 | # @brief Return max&min indexes (for slicing) of area specified. |
---|
885 | # @param latitudes_ref ?? |
---|
886 | # @param longitudes_ref ?? |
---|
887 | # @param minlat Minimum latitude of specified area. |
---|
888 | # @param maxlat Maximum latitude of specified area. |
---|
889 | # @param minlon Minimum longitude of specified area. |
---|
890 | # @param maxlon Maximum longitude of specified area. |
---|
891 | # @return Tuple (lat_min_index, lat_max_index, lon_min_index, lon_max_index) |
---|
892 | def get_min_max_indexes(latitudes_ref, longitudes_ref, |
---|
893 | minlat=None, maxlat=None, |
---|
894 | minlon=None, maxlon=None): |
---|
895 | """ |
---|
896 | Return max, min indexes (for slicing) of the lat and long arrays to cover |
---|
897 | the area specified with min/max lat/long. |
---|
898 | |
---|
899 | Think of the latitudes and longitudes describing a 2d surface. |
---|
900 | The area returned is, if possible, just big enough to cover the |
---|
901 | inputed max/min area. (This will not be possible if the max/min area |
---|
902 | has a section outside of the latitudes/longitudes area.) |
---|
903 | |
---|
904 | asset longitudes are sorted, |
---|
905 | long - from low to high (west to east, eg 148 - 151) |
---|
906 | assert latitudes are sorted, ascending or decending |
---|
907 | """ |
---|
908 | |
---|
909 | latitudes = latitudes_ref[:] |
---|
910 | longitudes = longitudes_ref[:] |
---|
911 | |
---|
912 | latitudes = ensure_numeric(latitudes) |
---|
913 | longitudes = ensure_numeric(longitudes) |
---|
914 | |
---|
915 | assert num.allclose(num.sort(longitudes), longitudes) |
---|
916 | |
---|
917 | #print latitudes[0],longitudes[0] |
---|
918 | #print len(latitudes),len(longitudes) |
---|
919 | #print latitudes[len(latitudes)-1],longitudes[len(longitudes)-1] |
---|
920 | |
---|
921 | lat_ascending = True |
---|
922 | if not num.allclose(num.sort(latitudes), latitudes): |
---|
923 | lat_ascending = False |
---|
924 | # reverse order of lat, so it's in ascending order |
---|
925 | latitudes = latitudes[::-1] |
---|
926 | assert num.allclose(num.sort(latitudes), latitudes) |
---|
927 | |
---|
928 | largest_lat_index = len(latitudes)-1 |
---|
929 | |
---|
930 | #Cut out a smaller extent. |
---|
931 | if minlat == None: |
---|
932 | lat_min_index = 0 |
---|
933 | else: |
---|
934 | lat_min_index = num.searchsorted(latitudes, minlat)-1 |
---|
935 | if lat_min_index < 0: |
---|
936 | lat_min_index = 0 |
---|
937 | |
---|
938 | if maxlat == None: |
---|
939 | lat_max_index = largest_lat_index #len(latitudes) |
---|
940 | else: |
---|
941 | lat_max_index = num.searchsorted(latitudes, maxlat) |
---|
942 | if lat_max_index > largest_lat_index: |
---|
943 | lat_max_index = largest_lat_index |
---|
944 | |
---|
945 | if minlon == None: |
---|
946 | lon_min_index = 0 |
---|
947 | else: |
---|
948 | lon_min_index = num.searchsorted(longitudes, minlon)-1 |
---|
949 | if lon_min_index < 0: |
---|
950 | lon_min_index = 0 |
---|
951 | |
---|
952 | if maxlon == None: |
---|
953 | lon_max_index = len(longitudes) |
---|
954 | else: |
---|
955 | lon_max_index = num.searchsorted(longitudes, maxlon) |
---|
956 | |
---|
957 | # Reversing the indexes, if the lat array is decending |
---|
958 | if lat_ascending is False: |
---|
959 | lat_min_index, lat_max_index = largest_lat_index - lat_max_index, \ |
---|
960 | largest_lat_index - lat_min_index |
---|
961 | lat_max_index = lat_max_index + 1 # taking into account how slicing works |
---|
962 | lon_max_index = lon_max_index + 1 # taking into account how slicing works |
---|
963 | |
---|
964 | return lat_min_index, lat_max_index, lon_min_index, lon_max_index |
---|
965 | |
---|
966 | |
---|