1 | """datamanager.py - input output for AnuGA |
---|
2 | |
---|
3 | |
---|
4 | This module takes care of reading and writing datafiles such as topograhies, |
---|
5 | model output, etc |
---|
6 | |
---|
7 | |
---|
8 | Formats used within AnuGA: |
---|
9 | |
---|
10 | .sww: Netcdf format for storing model output f(t,x,y) |
---|
11 | .tms: Netcdf format for storing time series f(t) |
---|
12 | |
---|
13 | .xya: ASCII format for storing arbitrary points and associated attributes |
---|
14 | .pts: NetCDF format for storing arbitrary points and associated attributes |
---|
15 | |
---|
16 | .asc: ASCII format of regular DEMs as output from ArcView |
---|
17 | .prj: Associated ArcView file giving more meta data for asc format |
---|
18 | .ers: ERMapper header format of regular DEMs for ArcView |
---|
19 | |
---|
20 | .dem: NetCDF representation of regular DEM data |
---|
21 | |
---|
22 | .tsh: ASCII format for storing meshes and associated boundary and region info |
---|
23 | .msh: NetCDF format for storing meshes and associated boundary and region info |
---|
24 | |
---|
25 | .nc: Native ferret NetCDF format |
---|
26 | .geo: Houdinis ascii geometry format (?) |
---|
27 | |
---|
28 | |
---|
29 | A typical dataflow can be described as follows |
---|
30 | |
---|
31 | Manually created files: |
---|
32 | ASC, PRJ: Digital elevation models (gridded) |
---|
33 | TSH: Triangular meshes (e.g. created from pmesh) |
---|
34 | NC Model outputs for use as boundary conditions (e.g from MOST) |
---|
35 | |
---|
36 | |
---|
37 | AUTOMATICALLY CREATED FILES: |
---|
38 | |
---|
39 | ASC, PRJ -> DEM -> PTS: Conversion of DEM's to native pts file |
---|
40 | |
---|
41 | NC -> SWW: Conversion of MOST bundary files to boundary sww |
---|
42 | |
---|
43 | PTS + TSH -> TSH with elevation: Least squares fit |
---|
44 | |
---|
45 | TSH -> SWW: Conversion of TSH to sww viewable using Swollen |
---|
46 | |
---|
47 | TSH + Boundary SWW -> SWW: Simluation using pyvolution |
---|
48 | |
---|
49 | """ |
---|
50 | |
---|
51 | from Numeric import concatenate |
---|
52 | |
---|
53 | from coordinate_transforms.geo_reference import Geo_reference |
---|
54 | |
---|
55 | def make_filename(s): |
---|
56 | """Transform argument string into a suitable filename |
---|
57 | """ |
---|
58 | |
---|
59 | s = s.strip() |
---|
60 | s = s.replace(' ', '_') |
---|
61 | s = s.replace('(', '') |
---|
62 | s = s.replace(')', '') |
---|
63 | s = s.replace('__', '_') |
---|
64 | |
---|
65 | return s |
---|
66 | |
---|
67 | |
---|
68 | def check_dir(path, verbose=None): |
---|
69 | """Check that specified path exists. |
---|
70 | If path does not exist it will be created if possible |
---|
71 | |
---|
72 | USAGE: |
---|
73 | checkdir(path, verbose): |
---|
74 | |
---|
75 | ARGUMENTS: |
---|
76 | path -- Directory |
---|
77 | verbose -- Flag verbose output (default: None) |
---|
78 | |
---|
79 | RETURN VALUE: |
---|
80 | Verified path including trailing separator |
---|
81 | |
---|
82 | """ |
---|
83 | |
---|
84 | import os, sys |
---|
85 | import os.path |
---|
86 | |
---|
87 | if sys.platform in ['nt', 'dos', 'win32', 'what else?']: |
---|
88 | unix = 0 |
---|
89 | else: |
---|
90 | unix = 1 |
---|
91 | |
---|
92 | |
---|
93 | if path[-1] != os.sep: |
---|
94 | path = path + os.sep # Add separator for directories |
---|
95 | |
---|
96 | path = os.path.expanduser(path) # Expand ~ or ~user in pathname |
---|
97 | if not (os.access(path,os.R_OK and os.W_OK) or path == ''): |
---|
98 | try: |
---|
99 | exitcode=os.mkdir(path) |
---|
100 | |
---|
101 | # Change access rights if possible |
---|
102 | # |
---|
103 | if unix: |
---|
104 | exitcode=os.system('chmod 775 '+path) |
---|
105 | else: |
---|
106 | pass # FIXME: What about acces rights under Windows? |
---|
107 | |
---|
108 | if verbose: print 'MESSAGE: Directory', path, 'created.' |
---|
109 | |
---|
110 | except: |
---|
111 | print 'WARNING: Directory', path, 'could not be created.' |
---|
112 | if unix: |
---|
113 | path = '/tmp/' |
---|
114 | else: |
---|
115 | path = 'C:' |
---|
116 | |
---|
117 | print 'Using directory %s instead' %path |
---|
118 | |
---|
119 | return(path) |
---|
120 | |
---|
121 | |
---|
122 | |
---|
123 | def del_dir(path): |
---|
124 | """Recursively delete directory path and all its contents |
---|
125 | """ |
---|
126 | |
---|
127 | import os |
---|
128 | |
---|
129 | if os.path.isdir(path): |
---|
130 | for file in os.listdir(path): |
---|
131 | X = os.path.join(path, file) |
---|
132 | |
---|
133 | |
---|
134 | if os.path.isdir(X) and not os.path.islink(X): |
---|
135 | del_dir(X) |
---|
136 | else: |
---|
137 | try: |
---|
138 | os.remove(X) |
---|
139 | except: |
---|
140 | print "Could not remove file %s" %X |
---|
141 | |
---|
142 | os.rmdir(path) |
---|
143 | |
---|
144 | |
---|
145 | |
---|
146 | def create_filename(datadir, filename, format, size=None, time=None): |
---|
147 | |
---|
148 | import os |
---|
149 | #from config import data_dir |
---|
150 | |
---|
151 | FN = check_dir(datadir) + filename |
---|
152 | |
---|
153 | if size is not None: |
---|
154 | FN += '_size%d' %size |
---|
155 | |
---|
156 | if time is not None: |
---|
157 | FN += '_time%.2f' %time |
---|
158 | |
---|
159 | FN += '.' + format |
---|
160 | return FN |
---|
161 | |
---|
162 | |
---|
163 | def get_files(datadir, filename, format, size): |
---|
164 | """Get all file (names) with gven name, size and format |
---|
165 | """ |
---|
166 | |
---|
167 | import glob |
---|
168 | |
---|
169 | import os |
---|
170 | #from config import data_dir |
---|
171 | |
---|
172 | dir = check_dir(datadir) |
---|
173 | |
---|
174 | pattern = dir + os.sep + filename + '_size=%d*.%s' %(size, format) |
---|
175 | return glob.glob(pattern) |
---|
176 | |
---|
177 | |
---|
178 | |
---|
179 | #Generic class for storing output to e.g. visualisation or checkpointing |
---|
180 | class Data_format: |
---|
181 | """Generic interface to data formats |
---|
182 | """ |
---|
183 | |
---|
184 | |
---|
185 | def __init__(self, domain, extension, mode = 'w'): |
---|
186 | assert mode in ['r', 'w', 'a'], '''Mode %s must be either:''' %mode +\ |
---|
187 | ''' 'w' (write)'''+\ |
---|
188 | ''' 'r' (read)''' +\ |
---|
189 | ''' 'a' (append)''' |
---|
190 | |
---|
191 | #Create filename |
---|
192 | #self.filename = create_filename(domain.get_datadir(), |
---|
193 | #domain.get_name(), extension, len(domain)) |
---|
194 | |
---|
195 | |
---|
196 | self.filename = create_filename(domain.get_datadir(), |
---|
197 | domain.get_name(), extension) |
---|
198 | |
---|
199 | #print 'F', self.filename |
---|
200 | self.timestep = 0 |
---|
201 | self.number_of_volumes = len(domain) |
---|
202 | self.domain = domain |
---|
203 | |
---|
204 | |
---|
205 | #FIXME: Should we have a general set_precision function? |
---|
206 | |
---|
207 | |
---|
208 | |
---|
209 | #Class for storing output to e.g. visualisation |
---|
210 | class Data_format_sww(Data_format): |
---|
211 | """Interface to native NetCDF format (.sww) for storing model output |
---|
212 | |
---|
213 | There are two kinds of data |
---|
214 | |
---|
215 | 1: Constant data: Vertex coordinates and field values. Stored once |
---|
216 | 2: Variable data: Conserved quantities. Stored once per timestep. |
---|
217 | |
---|
218 | All data is assumed to reside at vertex locations. |
---|
219 | """ |
---|
220 | |
---|
221 | |
---|
222 | def __init__(self, domain, mode = 'w',\ |
---|
223 | max_size = 2000000000, |
---|
224 | recursion = False): |
---|
225 | from Scientific.IO.NetCDF import NetCDFFile |
---|
226 | from Numeric import Int, Float, Float32 |
---|
227 | |
---|
228 | self.precision = Float32 #Use single precision |
---|
229 | if hasattr(domain, 'max_size'): |
---|
230 | self.max_size = domain.max_size #file size max is 2Gig |
---|
231 | else: |
---|
232 | self.max_size = max_size |
---|
233 | self.recursion = recursion |
---|
234 | self.mode = mode |
---|
235 | |
---|
236 | Data_format.__init__(self, domain, 'sww', mode) |
---|
237 | |
---|
238 | |
---|
239 | # NetCDF file definition |
---|
240 | fid = NetCDFFile(self.filename, mode) |
---|
241 | |
---|
242 | if mode == 'w': |
---|
243 | |
---|
244 | #Create new file |
---|
245 | fid.institution = 'Geoscience Australia' |
---|
246 | fid.description = 'Output from pyvolution suitable for plotting' |
---|
247 | |
---|
248 | if domain.smooth: |
---|
249 | fid.smoothing = 'Yes' |
---|
250 | else: |
---|
251 | fid.smoothing = 'No' |
---|
252 | |
---|
253 | fid.order = domain.default_order |
---|
254 | |
---|
255 | if hasattr(domain, 'texture'): |
---|
256 | fid.texture = domain.texture |
---|
257 | #else: |
---|
258 | # fid.texture = 'None' |
---|
259 | |
---|
260 | #Reference point |
---|
261 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
262 | #FIXME: Use Georef |
---|
263 | fid.starttime = domain.starttime |
---|
264 | |
---|
265 | # dimension definitions |
---|
266 | fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
267 | fid.createDimension('number_of_vertices', 3) |
---|
268 | |
---|
269 | if domain.smooth is True: |
---|
270 | fid.createDimension('number_of_points', len(domain.vertexlist)) |
---|
271 | else: |
---|
272 | fid.createDimension('number_of_points', 3*self.number_of_volumes) |
---|
273 | |
---|
274 | fid.createDimension('number_of_timesteps', None) #extensible |
---|
275 | |
---|
276 | # variable definitions |
---|
277 | fid.createVariable('x', self.precision, ('number_of_points',)) |
---|
278 | fid.createVariable('y', self.precision, ('number_of_points',)) |
---|
279 | fid.createVariable('elevation', self.precision, ('number_of_points',)) |
---|
280 | if domain.geo_reference is not None: |
---|
281 | domain.geo_reference.write_NetCDF(fid) |
---|
282 | |
---|
283 | #FIXME: Backwards compatibility |
---|
284 | fid.createVariable('z', self.precision, ('number_of_points',)) |
---|
285 | ################################# |
---|
286 | |
---|
287 | fid.createVariable('volumes', Int, ('number_of_volumes', |
---|
288 | 'number_of_vertices')) |
---|
289 | |
---|
290 | fid.createVariable('time', self.precision, |
---|
291 | ('number_of_timesteps',)) |
---|
292 | |
---|
293 | fid.createVariable('stage', self.precision, |
---|
294 | ('number_of_timesteps', |
---|
295 | 'number_of_points')) |
---|
296 | |
---|
297 | fid.createVariable('xmomentum', self.precision, |
---|
298 | ('number_of_timesteps', |
---|
299 | 'number_of_points')) |
---|
300 | |
---|
301 | fid.createVariable('ymomentum', self.precision, |
---|
302 | ('number_of_timesteps', |
---|
303 | 'number_of_points')) |
---|
304 | |
---|
305 | #Close |
---|
306 | fid.close() |
---|
307 | |
---|
308 | |
---|
309 | def store_connectivity(self): |
---|
310 | """Specialisation of store_connectivity for net CDF format |
---|
311 | |
---|
312 | Writes x,y,z coordinates of triangles constituting |
---|
313 | the bed elevation. |
---|
314 | """ |
---|
315 | |
---|
316 | from Scientific.IO.NetCDF import NetCDFFile |
---|
317 | |
---|
318 | from Numeric import concatenate, Int |
---|
319 | |
---|
320 | domain = self.domain |
---|
321 | |
---|
322 | #Get NetCDF |
---|
323 | fid = NetCDFFile(self.filename, 'a') #Open existing file for append |
---|
324 | |
---|
325 | # Get the variables |
---|
326 | x = fid.variables['x'] |
---|
327 | y = fid.variables['y'] |
---|
328 | z = fid.variables['elevation'] |
---|
329 | |
---|
330 | volumes = fid.variables['volumes'] |
---|
331 | |
---|
332 | # Get X, Y and bed elevation Z |
---|
333 | Q = domain.quantities['elevation'] |
---|
334 | X,Y,Z,V = Q.get_vertex_values(xy=True, |
---|
335 | precision = self.precision) |
---|
336 | |
---|
337 | |
---|
338 | |
---|
339 | x[:] = X.astype(self.precision) |
---|
340 | y[:] = Y.astype(self.precision) |
---|
341 | z[:] = Z.astype(self.precision) |
---|
342 | |
---|
343 | #FIXME: Backwards compatibility |
---|
344 | z = fid.variables['z'] |
---|
345 | z[:] = Z.astype(self.precision) |
---|
346 | ################################ |
---|
347 | |
---|
348 | volumes[:] = V.astype(volumes.typecode()) |
---|
349 | |
---|
350 | #Close |
---|
351 | fid.close() |
---|
352 | |
---|
353 | def store_timestep(self, names): |
---|
354 | """Store time and named quantities to file |
---|
355 | """ |
---|
356 | from Scientific.IO.NetCDF import NetCDFFile |
---|
357 | import types |
---|
358 | from time import sleep |
---|
359 | from os import stat |
---|
360 | |
---|
361 | minimum_allowed_depth = 0.001 |
---|
362 | #minimum_allowed_depth = 0.0 #FIXME pass in or read from domain |
---|
363 | from Numeric import choose |
---|
364 | |
---|
365 | #Get NetCDF |
---|
366 | retries = 0 |
---|
367 | file_open = False |
---|
368 | while not file_open and retries < 10: |
---|
369 | try: |
---|
370 | fid = NetCDFFile(self.filename, 'a') #Open existing file |
---|
371 | except IOError: |
---|
372 | #This could happen if someone was reading the file. |
---|
373 | #In that case, wait a while and try again |
---|
374 | msg = 'Warning (store_timestep): File %s could not be opened'\ |
---|
375 | %self.filename |
---|
376 | msg += ' - trying step %s again' %self.domain.time |
---|
377 | print msg |
---|
378 | retries += 1 |
---|
379 | sleep(1) |
---|
380 | else: |
---|
381 | file_open = True |
---|
382 | |
---|
383 | if not file_open: |
---|
384 | msg = 'File %s could not be opened for append' %self.filename |
---|
385 | raise msg |
---|
386 | |
---|
387 | |
---|
388 | |
---|
389 | #Check to see if the file is already too big: |
---|
390 | time = fid.variables['time'] |
---|
391 | i = len(time)+1 |
---|
392 | file_size = stat(self.filename)[6] |
---|
393 | file_size_increase = file_size/i |
---|
394 | if file_size + file_size_increase > self.max_size*(2**self.recursion): |
---|
395 | #in order to get the file name and start time correct, |
---|
396 | #I change the domian.filename and domain.starttime. |
---|
397 | #This is the only way to do this without changing |
---|
398 | #other modules (I think). |
---|
399 | |
---|
400 | #write a filename addon that won't break swollens reader |
---|
401 | #(10.sww is bad) |
---|
402 | filename_ext = '_time_%s'%self.domain.time |
---|
403 | filename_ext = filename_ext.replace('.', '_') |
---|
404 | #remember the old filename, then give domain a |
---|
405 | #name with the extension |
---|
406 | old_domain_filename = self.domain.filename |
---|
407 | if not self.recursion: |
---|
408 | self.domain.filename = self.domain.filename+filename_ext |
---|
409 | |
---|
410 | #change the domain starttime to the current time |
---|
411 | old_domain_starttime = self.domain.starttime |
---|
412 | self.domain.starttime = self.domain.time |
---|
413 | |
---|
414 | #build a new data_structure. |
---|
415 | next_data_structure=\ |
---|
416 | Data_format_sww(self.domain, mode=self.mode,\ |
---|
417 | max_size = self.max_size,\ |
---|
418 | recursion = self.recursion+1) |
---|
419 | if not self.recursion: |
---|
420 | print ' file_size = %s'%file_size |
---|
421 | print ' saving file to %s'%next_data_structure.filename |
---|
422 | #set up the new data_structure |
---|
423 | self.domain.writer = next_data_structure |
---|
424 | |
---|
425 | #FIXME - could be cleaner to use domain.store_timestep etc. |
---|
426 | next_data_structure.store_connectivity() |
---|
427 | next_data_structure.store_timestep(names) |
---|
428 | fid.sync() |
---|
429 | fid.close() |
---|
430 | |
---|
431 | #restore the old starttime and filename |
---|
432 | self.domain.starttime = old_domain_starttime |
---|
433 | self.domain.filename = old_domain_filename |
---|
434 | else: |
---|
435 | self.recursion = False |
---|
436 | domain = self.domain |
---|
437 | |
---|
438 | # Get the variables |
---|
439 | time = fid.variables['time'] |
---|
440 | stage = fid.variables['stage'] |
---|
441 | xmomentum = fid.variables['xmomentum'] |
---|
442 | ymomentum = fid.variables['ymomentum'] |
---|
443 | i = len(time) |
---|
444 | |
---|
445 | #Store time |
---|
446 | time[i] = self.domain.time |
---|
447 | |
---|
448 | |
---|
449 | if type(names) not in [types.ListType, types.TupleType]: |
---|
450 | names = [names] |
---|
451 | |
---|
452 | for name in names: |
---|
453 | # Get quantity |
---|
454 | Q = domain.quantities[name] |
---|
455 | A,V = Q.get_vertex_values(xy = False, |
---|
456 | precision = self.precision) |
---|
457 | |
---|
458 | #FIXME: Make this general (see below) |
---|
459 | if name == 'stage': |
---|
460 | z = fid.variables['elevation'] |
---|
461 | #print z[:] |
---|
462 | #print A-z[:] |
---|
463 | A = choose( A-z[:] >= minimum_allowed_depth, (z[:], A)) |
---|
464 | stage[i,:] = A.astype(self.precision) |
---|
465 | elif name == 'xmomentum': |
---|
466 | xmomentum[i,:] = A.astype(self.precision) |
---|
467 | elif name == 'ymomentum': |
---|
468 | ymomentum[i,:] = A.astype(self.precision) |
---|
469 | |
---|
470 | #As in.... |
---|
471 | #eval( name + '[i,:] = A.astype(self.precision)' ) |
---|
472 | #FIXME: But we need a UNIT test for that before refactoring |
---|
473 | |
---|
474 | |
---|
475 | |
---|
476 | #Flush and close |
---|
477 | fid.sync() |
---|
478 | fid.close() |
---|
479 | |
---|
480 | |
---|
481 | |
---|
482 | #Class for handling checkpoints data |
---|
483 | class Data_format_cpt(Data_format): |
---|
484 | """Interface to native NetCDF format (.cpt) |
---|
485 | """ |
---|
486 | |
---|
487 | |
---|
488 | def __init__(self, domain, mode = 'w'): |
---|
489 | from Scientific.IO.NetCDF import NetCDFFile |
---|
490 | from Numeric import Int, Float, Float |
---|
491 | |
---|
492 | self.precision = Float #Use full precision |
---|
493 | |
---|
494 | Data_format.__init__(self, domain, 'sww', mode) |
---|
495 | |
---|
496 | |
---|
497 | # NetCDF file definition |
---|
498 | fid = NetCDFFile(self.filename, mode) |
---|
499 | |
---|
500 | if mode == 'w': |
---|
501 | #Create new file |
---|
502 | fid.institution = 'Geoscience Australia' |
---|
503 | fid.description = 'Checkpoint data' |
---|
504 | #fid.smooth = domain.smooth |
---|
505 | fid.order = domain.default_order |
---|
506 | |
---|
507 | # dimension definitions |
---|
508 | fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
509 | fid.createDimension('number_of_vertices', 3) |
---|
510 | |
---|
511 | #Store info at all vertices (no smoothing) |
---|
512 | fid.createDimension('number_of_points', 3*self.number_of_volumes) |
---|
513 | fid.createDimension('number_of_timesteps', None) #extensible |
---|
514 | |
---|
515 | # variable definitions |
---|
516 | |
---|
517 | #Mesh |
---|
518 | fid.createVariable('x', self.precision, ('number_of_points',)) |
---|
519 | fid.createVariable('y', self.precision, ('number_of_points',)) |
---|
520 | |
---|
521 | |
---|
522 | fid.createVariable('volumes', Int, ('number_of_volumes', |
---|
523 | 'number_of_vertices')) |
---|
524 | |
---|
525 | fid.createVariable('time', self.precision, |
---|
526 | ('number_of_timesteps',)) |
---|
527 | |
---|
528 | #Allocate space for all quantities |
---|
529 | for name in domain.quantities.keys(): |
---|
530 | fid.createVariable(name, self.precision, |
---|
531 | ('number_of_timesteps', |
---|
532 | 'number_of_points')) |
---|
533 | |
---|
534 | #Close |
---|
535 | fid.close() |
---|
536 | |
---|
537 | |
---|
538 | def store_checkpoint(self): |
---|
539 | """ |
---|
540 | Write x,y coordinates of triangles. |
---|
541 | Write connectivity ( |
---|
542 | constituting |
---|
543 | the bed elevation. |
---|
544 | """ |
---|
545 | |
---|
546 | from Scientific.IO.NetCDF import NetCDFFile |
---|
547 | |
---|
548 | from Numeric import concatenate |
---|
549 | |
---|
550 | domain = self.domain |
---|
551 | |
---|
552 | #Get NetCDF |
---|
553 | fid = NetCDFFile(self.filename, 'a') #Open existing file for append |
---|
554 | |
---|
555 | # Get the variables |
---|
556 | x = fid.variables['x'] |
---|
557 | y = fid.variables['y'] |
---|
558 | |
---|
559 | volumes = fid.variables['volumes'] |
---|
560 | |
---|
561 | # Get X, Y and bed elevation Z |
---|
562 | Q = domain.quantities['elevation'] |
---|
563 | X,Y,Z,V = Q.get_vertex_values(xy=True, |
---|
564 | precision = self.precision) |
---|
565 | |
---|
566 | |
---|
567 | |
---|
568 | x[:] = X.astype(self.precision) |
---|
569 | y[:] = Y.astype(self.precision) |
---|
570 | z[:] = Z.astype(self.precision) |
---|
571 | |
---|
572 | volumes[:] = V |
---|
573 | |
---|
574 | #Close |
---|
575 | fid.close() |
---|
576 | |
---|
577 | |
---|
578 | def store_timestep(self, name): |
---|
579 | """Store time and named quantity to file |
---|
580 | """ |
---|
581 | from Scientific.IO.NetCDF import NetCDFFile |
---|
582 | from time import sleep |
---|
583 | |
---|
584 | #Get NetCDF |
---|
585 | retries = 0 |
---|
586 | file_open = False |
---|
587 | while not file_open and retries < 10: |
---|
588 | try: |
---|
589 | fid = NetCDFFile(self.filename, 'a') #Open existing file |
---|
590 | except IOError: |
---|
591 | #This could happen if someone was reading the file. |
---|
592 | #In that case, wait a while and try again |
---|
593 | msg = 'Warning (store_timestep): File %s could not be opened'\ |
---|
594 | %self.filename |
---|
595 | msg += ' - trying again' |
---|
596 | print msg |
---|
597 | retries += 1 |
---|
598 | sleep(1) |
---|
599 | else: |
---|
600 | file_open = True |
---|
601 | |
---|
602 | if not file_open: |
---|
603 | msg = 'File %s could not be opened for append' %self.filename |
---|
604 | raise msg |
---|
605 | |
---|
606 | |
---|
607 | domain = self.domain |
---|
608 | |
---|
609 | # Get the variables |
---|
610 | time = fid.variables['time'] |
---|
611 | stage = fid.variables['stage'] |
---|
612 | i = len(time) |
---|
613 | |
---|
614 | #Store stage |
---|
615 | time[i] = self.domain.time |
---|
616 | |
---|
617 | # Get quantity |
---|
618 | Q = domain.quantities[name] |
---|
619 | A,V = Q.get_vertex_values(xy=False, |
---|
620 | precision = self.precision) |
---|
621 | |
---|
622 | stage[i,:] = A.astype(self.precision) |
---|
623 | |
---|
624 | #Flush and close |
---|
625 | fid.sync() |
---|
626 | fid.close() |
---|
627 | |
---|
628 | |
---|
629 | |
---|
630 | |
---|
631 | |
---|
632 | #Function for storing xya output |
---|
633 | #FIXME Not done yet for this version |
---|
634 | class Data_format_xya(Data_format): |
---|
635 | """Generic interface to data formats |
---|
636 | """ |
---|
637 | |
---|
638 | |
---|
639 | def __init__(self, domain, mode = 'w'): |
---|
640 | from Scientific.IO.NetCDF import NetCDFFile |
---|
641 | from Numeric import Int, Float, Float32 |
---|
642 | |
---|
643 | self.precision = Float32 #Use single precision |
---|
644 | |
---|
645 | Data_format.__init__(self, domain, 'xya', mode) |
---|
646 | |
---|
647 | |
---|
648 | |
---|
649 | #FIXME -This is the old xya format |
---|
650 | def store_all(self): |
---|
651 | """Specialisation of store all for xya format |
---|
652 | |
---|
653 | Writes x,y,z coordinates of triangles constituting |
---|
654 | the bed elevation. |
---|
655 | """ |
---|
656 | |
---|
657 | from Numeric import concatenate |
---|
658 | |
---|
659 | domain = self.domain |
---|
660 | |
---|
661 | fd = open(self.filename, 'w') |
---|
662 | |
---|
663 | |
---|
664 | if domain.smooth is True: |
---|
665 | number_of_points = len(domain.vertexlist) |
---|
666 | else: |
---|
667 | number_of_points = 3*self.number_of_volumes |
---|
668 | |
---|
669 | numVertAttrib = 3 #Three attributes is what is assumed by the xya format |
---|
670 | |
---|
671 | fd.write(str(number_of_points) + " " + str(numVertAttrib) +\ |
---|
672 | " # <vertex #> <x> <y> [attributes]" + "\n") |
---|
673 | |
---|
674 | |
---|
675 | # Get X, Y, bed elevation and friction (index=0,1) |
---|
676 | X,Y,A,V = domain.get_vertex_values(xy=True, value_array='field_values', |
---|
677 | indices = (0,1), precision = self.precision) |
---|
678 | |
---|
679 | bed_eles = A[:,0] |
---|
680 | fricts = A[:,1] |
---|
681 | |
---|
682 | # Get stage (index=0) |
---|
683 | B,V = domain.get_vertex_values(xy=False, value_array='conserved_quantities', |
---|
684 | indices = (0,), precision = self.precision) |
---|
685 | |
---|
686 | stages = B[:,0] |
---|
687 | |
---|
688 | #<vertex #> <x> <y> [attributes] |
---|
689 | for x, y, bed_ele, stage, frict in map(None, X, Y, bed_eles, |
---|
690 | stages, fricts): |
---|
691 | |
---|
692 | s = '%.6f %.6f %.6f %.6f %.6f\n' %(x, y, bed_ele, stage, frict) |
---|
693 | fd.write(s) |
---|
694 | |
---|
695 | #close |
---|
696 | fd.close() |
---|
697 | |
---|
698 | |
---|
699 | def store_timestep(self, t, V0, V1, V2): |
---|
700 | """Store time, water heights (and momentums) to file |
---|
701 | """ |
---|
702 | pass |
---|
703 | |
---|
704 | |
---|
705 | #Auxiliary |
---|
706 | def write_obj(filename,x,y,z): |
---|
707 | """Store x,y,z vectors into filename (obj format) |
---|
708 | Vectors are assumed to have dimension (M,3) where |
---|
709 | M corresponds to the number elements. |
---|
710 | triangles are assumed to be disconnected |
---|
711 | |
---|
712 | The three numbers in each vector correspond to three vertices, |
---|
713 | |
---|
714 | e.g. the x coordinate of vertex 1 of element i is in x[i,1] |
---|
715 | |
---|
716 | """ |
---|
717 | #print 'Writing obj to %s' % filename |
---|
718 | |
---|
719 | import os.path |
---|
720 | |
---|
721 | root, ext = os.path.splitext(filename) |
---|
722 | if ext == '.obj': |
---|
723 | FN = filename |
---|
724 | else: |
---|
725 | FN = filename + '.obj' |
---|
726 | |
---|
727 | |
---|
728 | outfile = open(FN, 'wb') |
---|
729 | outfile.write("# Triangulation as an obj file\n") |
---|
730 | |
---|
731 | M, N = x.shape |
---|
732 | assert N==3 #Assuming three vertices per element |
---|
733 | |
---|
734 | for i in range(M): |
---|
735 | for j in range(N): |
---|
736 | outfile.write("v %f %f %f\n" % (x[i,j],y[i,j],z[i,j])) |
---|
737 | |
---|
738 | for i in range(M): |
---|
739 | base = i*N |
---|
740 | outfile.write("f %d %d %d\n" % (base+1,base+2,base+3)) |
---|
741 | |
---|
742 | outfile.close() |
---|
743 | |
---|
744 | |
---|
745 | ######################################################### |
---|
746 | #Conversion routines |
---|
747 | ######################################################## |
---|
748 | |
---|
749 | def sww2obj(basefilename, size): |
---|
750 | """Convert netcdf based data output to obj |
---|
751 | """ |
---|
752 | from Scientific.IO.NetCDF import NetCDFFile |
---|
753 | |
---|
754 | from Numeric import Float, zeros |
---|
755 | |
---|
756 | #Get NetCDF |
---|
757 | FN = create_filename('.', basefilename, 'sww', size) |
---|
758 | print 'Reading from ', FN |
---|
759 | fid = NetCDFFile(FN, 'r') #Open existing file for read |
---|
760 | |
---|
761 | |
---|
762 | # Get the variables |
---|
763 | x = fid.variables['x'] |
---|
764 | y = fid.variables['y'] |
---|
765 | z = fid.variables['elevation'] |
---|
766 | time = fid.variables['time'] |
---|
767 | stage = fid.variables['stage'] |
---|
768 | |
---|
769 | M = size #Number of lines |
---|
770 | xx = zeros((M,3), Float) |
---|
771 | yy = zeros((M,3), Float) |
---|
772 | zz = zeros((M,3), Float) |
---|
773 | |
---|
774 | for i in range(M): |
---|
775 | for j in range(3): |
---|
776 | xx[i,j] = x[i+j*M] |
---|
777 | yy[i,j] = y[i+j*M] |
---|
778 | zz[i,j] = z[i+j*M] |
---|
779 | |
---|
780 | #Write obj for bathymetry |
---|
781 | FN = create_filename('.', basefilename, 'obj', size) |
---|
782 | write_obj(FN,xx,yy,zz) |
---|
783 | |
---|
784 | |
---|
785 | #Now read all the data with variable information, combine with |
---|
786 | #x,y info and store as obj |
---|
787 | |
---|
788 | for k in range(len(time)): |
---|
789 | t = time[k] |
---|
790 | print 'Processing timestep %f' %t |
---|
791 | |
---|
792 | for i in range(M): |
---|
793 | for j in range(3): |
---|
794 | zz[i,j] = stage[k,i+j*M] |
---|
795 | |
---|
796 | |
---|
797 | #Write obj for variable data |
---|
798 | #FN = create_filename(basefilename, 'obj', size, time=t) |
---|
799 | FN = create_filename('.', basefilename[:5], 'obj', size, time=t) |
---|
800 | write_obj(FN,xx,yy,zz) |
---|
801 | |
---|
802 | |
---|
803 | def dat2obj(basefilename): |
---|
804 | """Convert line based data output to obj |
---|
805 | FIXME: Obsolete? |
---|
806 | """ |
---|
807 | |
---|
808 | import glob, os |
---|
809 | from config import data_dir |
---|
810 | |
---|
811 | |
---|
812 | #Get bathymetry and x,y's |
---|
813 | lines = open(data_dir+os.sep+basefilename+'_geometry.dat', 'r').readlines() |
---|
814 | |
---|
815 | from Numeric import zeros, Float |
---|
816 | |
---|
817 | M = len(lines) #Number of lines |
---|
818 | x = zeros((M,3), Float) |
---|
819 | y = zeros((M,3), Float) |
---|
820 | z = zeros((M,3), Float) |
---|
821 | |
---|
822 | ##i = 0 |
---|
823 | for i, line in enumerate(lines): |
---|
824 | tokens = line.split() |
---|
825 | values = map(float,tokens) |
---|
826 | |
---|
827 | for j in range(3): |
---|
828 | x[i,j] = values[j*3] |
---|
829 | y[i,j] = values[j*3+1] |
---|
830 | z[i,j] = values[j*3+2] |
---|
831 | |
---|
832 | ##i += 1 |
---|
833 | |
---|
834 | |
---|
835 | #Write obj for bathymetry |
---|
836 | write_obj(data_dir+os.sep+basefilename+'_geometry',x,y,z) |
---|
837 | |
---|
838 | |
---|
839 | #Now read all the data files with variable information, combine with |
---|
840 | #x,y info |
---|
841 | #and store as obj |
---|
842 | |
---|
843 | files = glob.glob(data_dir+os.sep+basefilename+'*.dat') |
---|
844 | |
---|
845 | for filename in files: |
---|
846 | print 'Processing %s' % filename |
---|
847 | |
---|
848 | lines = open(data_dir+os.sep+filename,'r').readlines() |
---|
849 | assert len(lines) == M |
---|
850 | root, ext = os.path.splitext(filename) |
---|
851 | |
---|
852 | #Get time from filename |
---|
853 | i0 = filename.find('_time=') |
---|
854 | if i0 == -1: |
---|
855 | #Skip bathymetry file |
---|
856 | continue |
---|
857 | |
---|
858 | i0 += 6 #Position where time starts |
---|
859 | i1 = filename.find('.dat') |
---|
860 | |
---|
861 | if i1 > i0: |
---|
862 | t = float(filename[i0:i1]) |
---|
863 | else: |
---|
864 | raise 'Hmmmm' |
---|
865 | |
---|
866 | |
---|
867 | |
---|
868 | ##i = 0 |
---|
869 | for i, line in enumerate(lines): |
---|
870 | tokens = line.split() |
---|
871 | values = map(float,tokens) |
---|
872 | |
---|
873 | for j in range(3): |
---|
874 | z[i,j] = values[j] |
---|
875 | |
---|
876 | ##i += 1 |
---|
877 | |
---|
878 | #Write obj for variable data |
---|
879 | write_obj(data_dir+os.sep+basefilename+'_time=%.4f' %t,x,y,z) |
---|
880 | |
---|
881 | |
---|
882 | def filter_netcdf(filename1, filename2, first=0, last=None, step = 1): |
---|
883 | """Read netcdf filename1, pick timesteps first:step:last and save to |
---|
884 | nettcdf file filename2 |
---|
885 | """ |
---|
886 | from Scientific.IO.NetCDF import NetCDFFile |
---|
887 | |
---|
888 | #Get NetCDF |
---|
889 | infile = NetCDFFile(filename1, 'r') #Open existing file for read |
---|
890 | outfile = NetCDFFile(filename2, 'w') #Open new file |
---|
891 | |
---|
892 | |
---|
893 | #Copy dimensions |
---|
894 | for d in infile.dimensions: |
---|
895 | outfile.createDimension(d, infile.dimensions[d]) |
---|
896 | |
---|
897 | for name in infile.variables: |
---|
898 | var = infile.variables[name] |
---|
899 | outfile.createVariable(name, var.typecode(), var.dimensions) |
---|
900 | |
---|
901 | |
---|
902 | #Copy the static variables |
---|
903 | for name in infile.variables: |
---|
904 | if name == 'time' or name == 'stage': |
---|
905 | pass |
---|
906 | else: |
---|
907 | #Copy |
---|
908 | outfile.variables[name][:] = infile.variables[name][:] |
---|
909 | |
---|
910 | #Copy selected timesteps |
---|
911 | time = infile.variables['time'] |
---|
912 | stage = infile.variables['stage'] |
---|
913 | |
---|
914 | newtime = outfile.variables['time'] |
---|
915 | newstage = outfile.variables['stage'] |
---|
916 | |
---|
917 | if last is None: |
---|
918 | last = len(time) |
---|
919 | |
---|
920 | selection = range(first, last, step) |
---|
921 | for i, j in enumerate(selection): |
---|
922 | print 'Copying timestep %d of %d (%f)' %(j, last-first, time[j]) |
---|
923 | newtime[i] = time[j] |
---|
924 | newstage[i,:] = stage[j,:] |
---|
925 | |
---|
926 | #Close |
---|
927 | infile.close() |
---|
928 | outfile.close() |
---|
929 | |
---|
930 | |
---|
931 | #Get data objects |
---|
932 | def get_dataobject(domain, mode='w'): |
---|
933 | """Return instance of class of given format using filename |
---|
934 | """ |
---|
935 | |
---|
936 | cls = eval('Data_format_%s' %domain.format) |
---|
937 | return cls(domain, mode) |
---|
938 | |
---|
939 | def xya2pts(basename_in, basename_out=None, verbose=False, |
---|
940 | #easting_min=None, easting_max=None, |
---|
941 | #northing_min=None, northing_max=None, |
---|
942 | stride = 1, |
---|
943 | attribute_name = 'elevation', |
---|
944 | z_func = None): |
---|
945 | """Read points data from ascii (.xya) |
---|
946 | |
---|
947 | Example: |
---|
948 | |
---|
949 | x(m) y(m) z(m) |
---|
950 | 0.00000e+00 0.00000e+00 1.3535000e-01 |
---|
951 | 0.00000e+00 1.40000e-02 1.3535000e-01 |
---|
952 | |
---|
953 | |
---|
954 | |
---|
955 | Convert to NetCDF pts format which is |
---|
956 | |
---|
957 | points: (Nx2) Float array |
---|
958 | elevation: N Float array |
---|
959 | |
---|
960 | Only lines that contain three numeric values are processed |
---|
961 | |
---|
962 | If z_func is specified, it will be applied to the third field |
---|
963 | """ |
---|
964 | |
---|
965 | import os |
---|
966 | #from Scientific.IO.NetCDF import NetCDFFile |
---|
967 | from Numeric import Float, arrayrange, concatenate |
---|
968 | |
---|
969 | root, ext = os.path.splitext(basename_in) |
---|
970 | |
---|
971 | if ext == '': ext = '.xya' |
---|
972 | |
---|
973 | #Get NetCDF |
---|
974 | infile = open(root + ext, 'r') #Open existing xya file for read |
---|
975 | |
---|
976 | if verbose: print 'Reading xya points from %s' %(root + ext) |
---|
977 | |
---|
978 | points = [] |
---|
979 | attribute = [] |
---|
980 | for i, line in enumerate(infile.readlines()): |
---|
981 | |
---|
982 | if i % stride != 0: continue |
---|
983 | |
---|
984 | fields = line.split() |
---|
985 | |
---|
986 | try: |
---|
987 | assert len(fields) == 3 |
---|
988 | except: |
---|
989 | print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) |
---|
990 | |
---|
991 | try: |
---|
992 | x = float( fields[0] ) |
---|
993 | y = float( fields[1] ) |
---|
994 | z = float( fields[2] ) |
---|
995 | except: |
---|
996 | continue |
---|
997 | |
---|
998 | points.append( [x, y] ) |
---|
999 | |
---|
1000 | if callable(z_func): |
---|
1001 | attribute.append(z_func(z)) |
---|
1002 | else: |
---|
1003 | attribute.append(z) |
---|
1004 | |
---|
1005 | |
---|
1006 | #Get output file |
---|
1007 | if basename_out == None: |
---|
1008 | ptsname = root + '.pts' |
---|
1009 | else: |
---|
1010 | ptsname = basename_out + '.pts' |
---|
1011 | |
---|
1012 | if verbose: print 'Store to NetCDF file %s' %ptsname |
---|
1013 | write_ptsfile(ptsname, points, attribute, attribute_name) |
---|
1014 | |
---|
1015 | |
---|
1016 | def write_ptsfile(ptsname, points, attribute, attribute_name = None): |
---|
1017 | """Write points and associated attribute to pts (NetCDF) format |
---|
1018 | """ |
---|
1019 | |
---|
1020 | from Numeric import Float |
---|
1021 | |
---|
1022 | if attribute_name is None: |
---|
1023 | attribute_name = 'attribute' |
---|
1024 | |
---|
1025 | |
---|
1026 | from Scientific.IO.NetCDF import NetCDFFile |
---|
1027 | |
---|
1028 | # NetCDF file definition |
---|
1029 | outfile = NetCDFFile(ptsname, 'w') |
---|
1030 | |
---|
1031 | |
---|
1032 | #Create new file |
---|
1033 | outfile.institution = 'Geoscience Australia' |
---|
1034 | outfile.description = 'NetCDF pts format for compact and '\ |
---|
1035 | 'portable storage of spatial point data' |
---|
1036 | |
---|
1037 | |
---|
1038 | #Georeferencing |
---|
1039 | from coordinate_transforms.geo_reference import Geo_reference |
---|
1040 | Geo_reference().write_NetCDF(outfile) |
---|
1041 | |
---|
1042 | |
---|
1043 | outfile.createDimension('number_of_points', len(points)) |
---|
1044 | outfile.createDimension('number_of_dimensions', 2) #This is 2d data |
---|
1045 | |
---|
1046 | # variable definitions |
---|
1047 | outfile.createVariable('points', Float, ('number_of_points', |
---|
1048 | 'number_of_dimensions')) |
---|
1049 | outfile.createVariable(attribute_name, Float, ('number_of_points',)) |
---|
1050 | |
---|
1051 | # Get handles to the variables |
---|
1052 | nc_points = outfile.variables['points'] |
---|
1053 | nc_attribute = outfile.variables[attribute_name] |
---|
1054 | |
---|
1055 | #Store data |
---|
1056 | nc_points[:, :] = points |
---|
1057 | nc_attribute[:] = attribute |
---|
1058 | |
---|
1059 | outfile.close() |
---|
1060 | |
---|
1061 | |
---|
1062 | def dem2pts(basename_in, basename_out=None, verbose=False, |
---|
1063 | easting_min=None, easting_max=None, |
---|
1064 | northing_min=None, northing_max=None): |
---|
1065 | """Read Digitial Elevation model from the following NetCDF format (.dem) |
---|
1066 | |
---|
1067 | Example: |
---|
1068 | |
---|
1069 | ncols 3121 |
---|
1070 | nrows 1800 |
---|
1071 | xllcorner 722000 |
---|
1072 | yllcorner 5893000 |
---|
1073 | cellsize 25 |
---|
1074 | NODATA_value -9999 |
---|
1075 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
1076 | |
---|
1077 | Convert to NetCDF pts format which is |
---|
1078 | |
---|
1079 | points: (Nx2) Float array |
---|
1080 | elevation: N Float array |
---|
1081 | """ |
---|
1082 | |
---|
1083 | #FIXME: Can this be written feasibly using write_pts? |
---|
1084 | |
---|
1085 | import os |
---|
1086 | from Scientific.IO.NetCDF import NetCDFFile |
---|
1087 | from Numeric import Float, zeros, reshape |
---|
1088 | |
---|
1089 | root = basename_in |
---|
1090 | |
---|
1091 | #Get NetCDF |
---|
1092 | infile = NetCDFFile(root + '.dem', 'r') #Open existing netcdf file for read |
---|
1093 | |
---|
1094 | if verbose: print 'Reading DEM from %s' %(root + '.dem') |
---|
1095 | |
---|
1096 | ncols = infile.ncols[0] |
---|
1097 | nrows = infile.nrows[0] |
---|
1098 | xllcorner = infile.xllcorner[0] #Easting of lower left corner |
---|
1099 | yllcorner = infile.yllcorner[0] #Northing of lower left corner |
---|
1100 | cellsize = infile.cellsize[0] |
---|
1101 | NODATA_value = infile.NODATA_value[0] |
---|
1102 | dem_elevation = infile.variables['elevation'] |
---|
1103 | |
---|
1104 | zone = infile.zone[0] |
---|
1105 | false_easting = infile.false_easting[0] |
---|
1106 | false_northing = infile.false_northing[0] |
---|
1107 | |
---|
1108 | #Text strings |
---|
1109 | projection = infile.projection |
---|
1110 | datum = infile.datum |
---|
1111 | units = infile.units |
---|
1112 | |
---|
1113 | |
---|
1114 | #Get output file |
---|
1115 | if basename_out == None: |
---|
1116 | ptsname = root + '.pts' |
---|
1117 | else: |
---|
1118 | ptsname = basename_out + '.pts' |
---|
1119 | |
---|
1120 | if verbose: print 'Store to NetCDF file %s' %ptsname |
---|
1121 | # NetCDF file definition |
---|
1122 | outfile = NetCDFFile(ptsname, 'w') |
---|
1123 | |
---|
1124 | #Create new file |
---|
1125 | outfile.institution = 'Geoscience Australia' |
---|
1126 | outfile.description = 'NetCDF pts format for compact and portable storage ' +\ |
---|
1127 | 'of spatial point data' |
---|
1128 | #assign default values |
---|
1129 | if easting_min is None: easting_min = xllcorner |
---|
1130 | if easting_max is None: easting_max = xllcorner + ncols*cellsize |
---|
1131 | if northing_min is None: northing_min = yllcorner |
---|
1132 | if northing_max is None: northing_max = yllcorner + nrows*cellsize |
---|
1133 | |
---|
1134 | #compute offsets to update georeferencing |
---|
1135 | easting_offset = xllcorner - easting_min |
---|
1136 | northing_offset = yllcorner - northing_min |
---|
1137 | |
---|
1138 | #Georeferencing |
---|
1139 | outfile.zone = zone |
---|
1140 | outfile.xllcorner = easting_min #Easting of lower left corner |
---|
1141 | outfile.yllcorner = northing_min #Northing of lower left corner |
---|
1142 | outfile.false_easting = false_easting |
---|
1143 | outfile.false_northing = false_northing |
---|
1144 | |
---|
1145 | outfile.projection = projection |
---|
1146 | outfile.datum = datum |
---|
1147 | outfile.units = units |
---|
1148 | |
---|
1149 | |
---|
1150 | #Grid info (FIXME: probably not going to be used, but heck) |
---|
1151 | outfile.ncols = ncols |
---|
1152 | outfile.nrows = nrows |
---|
1153 | |
---|
1154 | |
---|
1155 | # dimension definitions |
---|
1156 | nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize)) |
---|
1157 | ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize)) |
---|
1158 | outfile.createDimension('number_of_points', nrows_in_bounding_box*ncols_in_bounding_box) |
---|
1159 | outfile.createDimension('number_of_dimensions', 2) #This is 2d data |
---|
1160 | |
---|
1161 | # variable definitions |
---|
1162 | outfile.createVariable('points', Float, ('number_of_points', |
---|
1163 | 'number_of_dimensions')) |
---|
1164 | outfile.createVariable('elevation', Float, ('number_of_points',)) |
---|
1165 | |
---|
1166 | # Get handles to the variables |
---|
1167 | points = outfile.variables['points'] |
---|
1168 | elevation = outfile.variables['elevation'] |
---|
1169 | |
---|
1170 | dem_elevation_r = reshape(dem_elevation, (nrows, ncols)) |
---|
1171 | |
---|
1172 | #Store data |
---|
1173 | #FIXME: Could perhaps be faster using array operations (Fixed 27/7/05) |
---|
1174 | global_index = 0 |
---|
1175 | for i in range(nrows): |
---|
1176 | if verbose and i%((nrows+10)/10)==0: |
---|
1177 | print 'Processing row %d of %d' %(i, nrows) |
---|
1178 | |
---|
1179 | lower_index = global_index |
---|
1180 | tpoints = zeros((ncols_in_bounding_box, 2), Float) |
---|
1181 | telev = zeros(ncols_in_bounding_box, Float) |
---|
1182 | local_index = 0 |
---|
1183 | |
---|
1184 | y = (nrows-i)*cellsize + yllcorner |
---|
1185 | for j in range(ncols): |
---|
1186 | |
---|
1187 | x = j*cellsize + xllcorner |
---|
1188 | if easting_min <= x <= easting_max and \ |
---|
1189 | northing_min <= y <= northing_max: |
---|
1190 | tpoints[local_index, :] = [x-easting_min,y-northing_min] |
---|
1191 | telev[local_index] = dem_elevation_r[i, j] |
---|
1192 | global_index += 1 |
---|
1193 | local_index += 1 |
---|
1194 | |
---|
1195 | upper_index = global_index |
---|
1196 | |
---|
1197 | if upper_index == lower_index + ncols_in_bounding_box: |
---|
1198 | points[lower_index:upper_index, :] = tpoints |
---|
1199 | elevation[lower_index:upper_index] = telev |
---|
1200 | |
---|
1201 | assert global_index == nrows_in_bounding_box*ncols_in_bounding_box, 'index not equal to number of points' |
---|
1202 | |
---|
1203 | infile.close() |
---|
1204 | outfile.close() |
---|
1205 | |
---|
1206 | |
---|
1207 | |
---|
1208 | def sww2dem(basename_in, basename_out = None, |
---|
1209 | quantity = None, |
---|
1210 | timestep = None, |
---|
1211 | reduction = None, |
---|
1212 | cellsize = 10, |
---|
1213 | easting_min = None, |
---|
1214 | easting_max = None, |
---|
1215 | northing_min = None, |
---|
1216 | northing_max = None, |
---|
1217 | verbose = False, |
---|
1218 | origin = None, |
---|
1219 | datum = 'WGS84', |
---|
1220 | format = 'ers'): |
---|
1221 | |
---|
1222 | """Read SWW file and convert to Digitial Elevation model format (.asc or .ers) |
---|
1223 | |
---|
1224 | Example (ASC): |
---|
1225 | |
---|
1226 | ncols 3121 |
---|
1227 | nrows 1800 |
---|
1228 | xllcorner 722000 |
---|
1229 | yllcorner 5893000 |
---|
1230 | cellsize 25 |
---|
1231 | NODATA_value -9999 |
---|
1232 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
1233 | |
---|
1234 | Also write accompanying file with same basename_in but extension .prj |
---|
1235 | used to fix the UTM zone, datum, false northings and eastings. |
---|
1236 | |
---|
1237 | The prj format is assumed to be as |
---|
1238 | |
---|
1239 | Projection UTM |
---|
1240 | Zone 56 |
---|
1241 | Datum WGS84 |
---|
1242 | Zunits NO |
---|
1243 | Units METERS |
---|
1244 | Spheroid WGS84 |
---|
1245 | Xshift 0.0000000000 |
---|
1246 | Yshift 10000000.0000000000 |
---|
1247 | Parameters |
---|
1248 | |
---|
1249 | |
---|
1250 | if quantity is given, out values from quantity otherwise default to |
---|
1251 | elevation |
---|
1252 | |
---|
1253 | if timestep (an index) is given, output quantity at that timestep |
---|
1254 | |
---|
1255 | if reduction is given use that to reduce quantity over all timesteps. |
---|
1256 | |
---|
1257 | datum |
---|
1258 | |
---|
1259 | format can be either 'asc' or 'ers' |
---|
1260 | """ |
---|
1261 | |
---|
1262 | from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue |
---|
1263 | |
---|
1264 | msg = 'Format must be either asc or ers' |
---|
1265 | assert format.lower() in ['asc', 'ers'], msg |
---|
1266 | |
---|
1267 | false_easting = 500000 |
---|
1268 | false_northing = 10000000 |
---|
1269 | |
---|
1270 | if quantity is None: |
---|
1271 | quantity = 'elevation' |
---|
1272 | |
---|
1273 | if reduction is None: |
---|
1274 | reduction = max |
---|
1275 | |
---|
1276 | if basename_out is None: |
---|
1277 | basename_out = basename_in + '_%s' %quantity |
---|
1278 | |
---|
1279 | swwfile = basename_in + '.sww' |
---|
1280 | demfile = basename_out + '.' + format |
---|
1281 | # Note the use of a .ers extension is optional (write_ermapper_grid will |
---|
1282 | # deal with either option |
---|
1283 | |
---|
1284 | #Read sww file |
---|
1285 | if verbose: print 'Reading from %s' %swwfile |
---|
1286 | from Scientific.IO.NetCDF import NetCDFFile |
---|
1287 | fid = NetCDFFile(swwfile) |
---|
1288 | |
---|
1289 | #Get extent and reference |
---|
1290 | x = fid.variables['x'][:] |
---|
1291 | y = fid.variables['y'][:] |
---|
1292 | volumes = fid.variables['volumes'][:] |
---|
1293 | |
---|
1294 | number_of_timesteps = fid.dimensions['number_of_timesteps'] |
---|
1295 | number_of_points = fid.dimensions['number_of_points'] |
---|
1296 | if origin is None: |
---|
1297 | |
---|
1298 | #Get geo_reference |
---|
1299 | #sww files don't have to have a geo_ref |
---|
1300 | try: |
---|
1301 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
1302 | except AttributeError, e: |
---|
1303 | geo_reference = Geo_reference() #Default georef object |
---|
1304 | |
---|
1305 | xllcorner = geo_reference.get_xllcorner() |
---|
1306 | yllcorner = geo_reference.get_yllcorner() |
---|
1307 | zone = geo_reference.get_zone() |
---|
1308 | else: |
---|
1309 | zone = origin[0] |
---|
1310 | xllcorner = origin[1] |
---|
1311 | yllcorner = origin[2] |
---|
1312 | |
---|
1313 | |
---|
1314 | |
---|
1315 | #Get quantity and reduce if applicable |
---|
1316 | if verbose: print 'Reading quantity %s' %quantity |
---|
1317 | |
---|
1318 | if quantity.lower() == 'depth': |
---|
1319 | q = fid.variables['stage'][:] - fid.variables['elevation'][:] |
---|
1320 | else: |
---|
1321 | q = fid.variables[quantity][:] |
---|
1322 | |
---|
1323 | |
---|
1324 | if len(q.shape) == 2: |
---|
1325 | #q has a time component and needs to be reduced along |
---|
1326 | #the temporal dimension |
---|
1327 | if verbose: print 'Reducing quantity %s' %quantity |
---|
1328 | q_reduced = zeros( number_of_points, Float ) |
---|
1329 | |
---|
1330 | for k in range(number_of_points): |
---|
1331 | q_reduced[k] = reduction( q[:,k] ) |
---|
1332 | |
---|
1333 | q = q_reduced |
---|
1334 | |
---|
1335 | #Post condition: Now q has dimension: number_of_points |
---|
1336 | assert len(q.shape) == 1 |
---|
1337 | assert q.shape[0] == number_of_points |
---|
1338 | |
---|
1339 | |
---|
1340 | #Create grid and update xll/yll corner and x,y |
---|
1341 | |
---|
1342 | #Relative extent |
---|
1343 | if easting_min is None: |
---|
1344 | xmin = min(x) |
---|
1345 | else: |
---|
1346 | xmin = easting_min - xllcorner |
---|
1347 | |
---|
1348 | if easting_max is None: |
---|
1349 | xmax = max(x) |
---|
1350 | else: |
---|
1351 | xmax = easting_max - xllcorner |
---|
1352 | |
---|
1353 | if northing_min is None: |
---|
1354 | ymin = min(y) |
---|
1355 | else: |
---|
1356 | ymin = northing_min - yllcorner |
---|
1357 | |
---|
1358 | if northing_max is None: |
---|
1359 | ymax = max(y) |
---|
1360 | else: |
---|
1361 | ymax = northing_max - yllcorner |
---|
1362 | |
---|
1363 | |
---|
1364 | |
---|
1365 | if verbose: print 'Creating grid' |
---|
1366 | ncols = int((xmax-xmin)/cellsize)+1 |
---|
1367 | nrows = int((ymax-ymin)/cellsize)+1 |
---|
1368 | |
---|
1369 | |
---|
1370 | #New absolute reference and coordinates |
---|
1371 | newxllcorner = xmin+xllcorner |
---|
1372 | newyllcorner = ymin+yllcorner |
---|
1373 | |
---|
1374 | x = x+xllcorner-newxllcorner |
---|
1375 | y = y+yllcorner-newyllcorner |
---|
1376 | |
---|
1377 | vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) |
---|
1378 | assert len(vertex_points.shape) == 2 |
---|
1379 | |
---|
1380 | |
---|
1381 | from Numeric import zeros, Float |
---|
1382 | grid_points = zeros ( (ncols*nrows, 2), Float ) |
---|
1383 | |
---|
1384 | |
---|
1385 | for i in xrange(nrows): |
---|
1386 | if format.lower() == 'asc': |
---|
1387 | yg = i*cellsize |
---|
1388 | else: |
---|
1389 | #this will flip the order of the y values for ers |
---|
1390 | yg = (nrows-i)*cellsize |
---|
1391 | |
---|
1392 | for j in xrange(ncols): |
---|
1393 | xg = j*cellsize |
---|
1394 | k = i*ncols + j |
---|
1395 | |
---|
1396 | grid_points[k,0] = xg |
---|
1397 | grid_points[k,1] = yg |
---|
1398 | |
---|
1399 | #Interpolate |
---|
1400 | from least_squares import Interpolation |
---|
1401 | from util import inside_polygon |
---|
1402 | |
---|
1403 | #FIXME: This should be done with precrop = True (?), otherwise it'll |
---|
1404 | #take forever. With expand_search set to False, some grid points might |
---|
1405 | #miss out.... This will be addressed though Duncan's refactoring of least_squares |
---|
1406 | |
---|
1407 | interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0, |
---|
1408 | precrop = False, expand_search = False, |
---|
1409 | verbose = verbose) |
---|
1410 | |
---|
1411 | #Interpolate using quantity values |
---|
1412 | if verbose: print 'Interpolating' |
---|
1413 | grid_values = interp.interpolate(q).flat |
---|
1414 | |
---|
1415 | if format.lower() == 'ers': |
---|
1416 | # setup ERS header information |
---|
1417 | grid_values = reshape(grid_values,(nrows, ncols)) |
---|
1418 | NODATA_value = 0 |
---|
1419 | header = {} |
---|
1420 | header['datum'] = '"' + datum + '"' |
---|
1421 | # FIXME The use of hardwired UTM and zone number needs to be made optional |
---|
1422 | # FIXME Also need an automatic test for coordinate type (i.e. EN or LL) |
---|
1423 | header['projection'] = '"UTM-' + str(zone) + '"' |
---|
1424 | header['coordinatetype'] = 'EN' |
---|
1425 | if header['coordinatetype'] == 'LL': |
---|
1426 | header['longitude'] = str(newxllcorner) |
---|
1427 | header['latitude'] = str(newyllcorner) |
---|
1428 | elif header['coordinatetype'] == 'EN': |
---|
1429 | header['eastings'] = str(newxllcorner) |
---|
1430 | header['northings'] = str(newyllcorner) |
---|
1431 | header['nullcellvalue'] = str(NODATA_value) |
---|
1432 | header['xdimension'] = str(cellsize) |
---|
1433 | header['ydimension'] = str(cellsize) |
---|
1434 | header['value'] = '"' + quantity + '"' |
---|
1435 | |
---|
1436 | |
---|
1437 | #Write |
---|
1438 | if verbose: print 'Writing %s' %demfile |
---|
1439 | import ermapper_grids |
---|
1440 | ermapper_grids.write_ermapper_grid(demfile, grid_values, header) |
---|
1441 | |
---|
1442 | fid.close() |
---|
1443 | else: |
---|
1444 | #Write to Ascii format |
---|
1445 | |
---|
1446 | #Write prj file |
---|
1447 | prjfile = basename_out + '.prj' |
---|
1448 | |
---|
1449 | if verbose: print 'Writing %s' %prjfile |
---|
1450 | prjid = open(prjfile, 'w') |
---|
1451 | prjid.write('Projection %s\n' %'UTM') |
---|
1452 | prjid.write('Zone %d\n' %zone) |
---|
1453 | prjid.write('Datum %s\n' %datum) |
---|
1454 | prjid.write('Zunits NO\n') |
---|
1455 | prjid.write('Units METERS\n') |
---|
1456 | prjid.write('Spheroid %s\n' %datum) |
---|
1457 | prjid.write('Xshift %d\n' %false_easting) |
---|
1458 | prjid.write('Yshift %d\n' %false_northing) |
---|
1459 | prjid.write('Parameters\n') |
---|
1460 | prjid.close() |
---|
1461 | |
---|
1462 | |
---|
1463 | |
---|
1464 | if verbose: print 'Writing %s' %ascfile |
---|
1465 | NODATA_value = -9999 |
---|
1466 | |
---|
1467 | ascid = open(demfile, 'w') |
---|
1468 | |
---|
1469 | ascid.write('ncols %d\n' %ncols) |
---|
1470 | ascid.write('nrows %d\n' %nrows) |
---|
1471 | ascid.write('xllcorner %d\n' %newxllcorner) |
---|
1472 | ascid.write('yllcorner %d\n' %newyllcorner) |
---|
1473 | ascid.write('cellsize %f\n' %cellsize) |
---|
1474 | ascid.write('NODATA_value %d\n' %NODATA_value) |
---|
1475 | |
---|
1476 | |
---|
1477 | #Get bounding polygon from mesh |
---|
1478 | P = interp.mesh.get_boundary_polygon() |
---|
1479 | inside_indices = inside_polygon(grid_points, P) |
---|
1480 | |
---|
1481 | for i in range(nrows): |
---|
1482 | if verbose and i%((nrows+10)/10)==0: |
---|
1483 | print 'Doing row %d of %d' %(i, nrows) |
---|
1484 | |
---|
1485 | for j in range(ncols): |
---|
1486 | index = (nrows-i-1)*ncols+j |
---|
1487 | |
---|
1488 | if sometrue(inside_indices == index): |
---|
1489 | ascid.write('%f ' %grid_values[index]) |
---|
1490 | else: |
---|
1491 | ascid.write('%d ' %NODATA_value) |
---|
1492 | |
---|
1493 | ascid.write('\n') |
---|
1494 | |
---|
1495 | #Close |
---|
1496 | ascid.close() |
---|
1497 | fid.close() |
---|
1498 | |
---|
1499 | #Backwards compatibility |
---|
1500 | def sww2asc(basename_in, basename_out = None, |
---|
1501 | quantity = None, |
---|
1502 | timestep = None, |
---|
1503 | reduction = None, |
---|
1504 | cellsize = 10, |
---|
1505 | verbose = False, |
---|
1506 | origin = None): |
---|
1507 | print 'sww2asc will soon be obsoleted - please use sww2dem' |
---|
1508 | sww2dem(basename_in, |
---|
1509 | basename_out = basename_out, |
---|
1510 | quantity = quantity, |
---|
1511 | timestep = timestep, |
---|
1512 | reduction = reduction, |
---|
1513 | cellsize = cellsize, |
---|
1514 | verbose = verbose, |
---|
1515 | origin = origin, |
---|
1516 | datum = 'WGS84', |
---|
1517 | format = 'asc') |
---|
1518 | |
---|
1519 | def sww2ers(basename_in, basename_out = None, |
---|
1520 | quantity = None, |
---|
1521 | timestep = None, |
---|
1522 | reduction = None, |
---|
1523 | cellsize = 10, |
---|
1524 | verbose = False, |
---|
1525 | origin = None, |
---|
1526 | datum = 'WGS84'): |
---|
1527 | print 'sww2ers will soon be obsoleted - please use sww2dem' |
---|
1528 | sww2dem(basename_in, |
---|
1529 | basename_out = basename_out, |
---|
1530 | quantity = quantity, |
---|
1531 | timestep = timestep, |
---|
1532 | reduction = reduction, |
---|
1533 | cellsize = cellsize, |
---|
1534 | verbose = verbose, |
---|
1535 | origin = origin, |
---|
1536 | datum = datum, |
---|
1537 | format = 'ers') |
---|
1538 | ################################# END COMPATIBILITY ############## |
---|
1539 | |
---|
1540 | |
---|
1541 | |
---|
1542 | |
---|
1543 | def convert_dem_from_ascii2netcdf(basename_in, basename_out = None, |
---|
1544 | verbose=False): |
---|
1545 | """Read Digitial Elevation model from the following ASCII format (.asc) |
---|
1546 | |
---|
1547 | Example: |
---|
1548 | |
---|
1549 | ncols 3121 |
---|
1550 | nrows 1800 |
---|
1551 | xllcorner 722000 |
---|
1552 | yllcorner 5893000 |
---|
1553 | cellsize 25 |
---|
1554 | NODATA_value -9999 |
---|
1555 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
1556 | |
---|
1557 | Convert basename_in + '.asc' to NetCDF format (.dem) |
---|
1558 | mimicking the ASCII format closely. |
---|
1559 | |
---|
1560 | |
---|
1561 | An accompanying file with same basename_in but extension .prj must exist |
---|
1562 | and is used to fix the UTM zone, datum, false northings and eastings. |
---|
1563 | |
---|
1564 | The prj format is assumed to be as |
---|
1565 | |
---|
1566 | Projection UTM |
---|
1567 | Zone 56 |
---|
1568 | Datum WGS84 |
---|
1569 | Zunits NO |
---|
1570 | Units METERS |
---|
1571 | Spheroid WGS84 |
---|
1572 | Xshift 0.0000000000 |
---|
1573 | Yshift 10000000.0000000000 |
---|
1574 | Parameters |
---|
1575 | """ |
---|
1576 | |
---|
1577 | import os |
---|
1578 | from Scientific.IO.NetCDF import NetCDFFile |
---|
1579 | from Numeric import Float, array |
---|
1580 | |
---|
1581 | #root, ext = os.path.splitext(basename_in) |
---|
1582 | root = basename_in |
---|
1583 | |
---|
1584 | ########################################### |
---|
1585 | # Read Meta data |
---|
1586 | if verbose: print 'Reading METADATA from %s' %root + '.prj' |
---|
1587 | metadatafile = open(root + '.prj') |
---|
1588 | metalines = metadatafile.readlines() |
---|
1589 | metadatafile.close() |
---|
1590 | |
---|
1591 | L = metalines[0].strip().split() |
---|
1592 | assert L[0].strip().lower() == 'projection' |
---|
1593 | projection = L[1].strip() #TEXT |
---|
1594 | |
---|
1595 | L = metalines[1].strip().split() |
---|
1596 | assert L[0].strip().lower() == 'zone' |
---|
1597 | zone = int(L[1].strip()) |
---|
1598 | |
---|
1599 | L = metalines[2].strip().split() |
---|
1600 | assert L[0].strip().lower() == 'datum' |
---|
1601 | datum = L[1].strip() #TEXT |
---|
1602 | |
---|
1603 | L = metalines[3].strip().split() |
---|
1604 | assert L[0].strip().lower() == 'zunits' #IGNORE |
---|
1605 | zunits = L[1].strip() #TEXT |
---|
1606 | |
---|
1607 | L = metalines[4].strip().split() |
---|
1608 | assert L[0].strip().lower() == 'units' |
---|
1609 | units = L[1].strip() #TEXT |
---|
1610 | |
---|
1611 | L = metalines[5].strip().split() |
---|
1612 | assert L[0].strip().lower() == 'spheroid' #IGNORE |
---|
1613 | spheroid = L[1].strip() #TEXT |
---|
1614 | |
---|
1615 | L = metalines[6].strip().split() |
---|
1616 | assert L[0].strip().lower() == 'xshift' |
---|
1617 | false_easting = float(L[1].strip()) |
---|
1618 | |
---|
1619 | L = metalines[7].strip().split() |
---|
1620 | assert L[0].strip().lower() == 'yshift' |
---|
1621 | false_northing = float(L[1].strip()) |
---|
1622 | |
---|
1623 | #print false_easting, false_northing, zone, datum |
---|
1624 | |
---|
1625 | |
---|
1626 | ########################################### |
---|
1627 | #Read DEM data |
---|
1628 | |
---|
1629 | datafile = open(basename_in + '.asc') |
---|
1630 | |
---|
1631 | if verbose: print 'Reading DEM from %s' %(basename_in + '.asc') |
---|
1632 | lines = datafile.readlines() |
---|
1633 | datafile.close() |
---|
1634 | |
---|
1635 | if verbose: print 'Got', len(lines), ' lines' |
---|
1636 | |
---|
1637 | ncols = int(lines[0].split()[1].strip()) |
---|
1638 | nrows = int(lines[1].split()[1].strip()) |
---|
1639 | xllcorner = float(lines[2].split()[1].strip()) |
---|
1640 | yllcorner = float(lines[3].split()[1].strip()) |
---|
1641 | cellsize = float(lines[4].split()[1].strip()) |
---|
1642 | NODATA_value = int(lines[5].split()[1].strip()) |
---|
1643 | |
---|
1644 | assert len(lines) == nrows + 6 |
---|
1645 | |
---|
1646 | |
---|
1647 | ########################################## |
---|
1648 | |
---|
1649 | |
---|
1650 | if basename_out == None: |
---|
1651 | netcdfname = root + '.dem' |
---|
1652 | else: |
---|
1653 | netcdfname = basename_out + '.dem' |
---|
1654 | |
---|
1655 | if verbose: print 'Store to NetCDF file %s' %netcdfname |
---|
1656 | # NetCDF file definition |
---|
1657 | fid = NetCDFFile(netcdfname, 'w') |
---|
1658 | |
---|
1659 | #Create new file |
---|
1660 | fid.institution = 'Geoscience Australia' |
---|
1661 | fid.description = 'NetCDF DEM format for compact and portable storage ' +\ |
---|
1662 | 'of spatial point data' |
---|
1663 | |
---|
1664 | fid.ncols = ncols |
---|
1665 | fid.nrows = nrows |
---|
1666 | fid.xllcorner = xllcorner |
---|
1667 | fid.yllcorner = yllcorner |
---|
1668 | fid.cellsize = cellsize |
---|
1669 | fid.NODATA_value = NODATA_value |
---|
1670 | |
---|
1671 | fid.zone = zone |
---|
1672 | fid.false_easting = false_easting |
---|
1673 | fid.false_northing = false_northing |
---|
1674 | fid.projection = projection |
---|
1675 | fid.datum = datum |
---|
1676 | fid.units = units |
---|
1677 | |
---|
1678 | |
---|
1679 | # dimension definitions |
---|
1680 | fid.createDimension('number_of_rows', nrows) |
---|
1681 | fid.createDimension('number_of_columns', ncols) |
---|
1682 | |
---|
1683 | # variable definitions |
---|
1684 | fid.createVariable('elevation', Float, ('number_of_rows', |
---|
1685 | 'number_of_columns')) |
---|
1686 | |
---|
1687 | # Get handles to the variables |
---|
1688 | elevation = fid.variables['elevation'] |
---|
1689 | |
---|
1690 | #Store data |
---|
1691 | n = len(lines[6:]) |
---|
1692 | for i, line in enumerate(lines[6:]): |
---|
1693 | fields = line.split() |
---|
1694 | if verbose and i%((n+10)/10)==0: |
---|
1695 | print 'Processing row %d of %d' %(i, nrows) |
---|
1696 | |
---|
1697 | elevation[i, :] = array([float(x) for x in fields]) |
---|
1698 | |
---|
1699 | fid.close() |
---|
1700 | |
---|
1701 | |
---|
1702 | |
---|
1703 | |
---|
1704 | |
---|
1705 | def ferret2sww(basename_in, basename_out = None, |
---|
1706 | verbose = False, |
---|
1707 | minlat = None, maxlat = None, |
---|
1708 | minlon = None, maxlon = None, |
---|
1709 | mint = None, maxt = None, mean_stage = 0, |
---|
1710 | origin = None, zscale = 1, |
---|
1711 | fail_on_NaN = True, |
---|
1712 | NaN_filler = 0, |
---|
1713 | elevation = None, |
---|
1714 | inverted_bathymetry = False |
---|
1715 | ): #FIXME: Bathymetry should be obtained |
---|
1716 | #from MOST somehow. |
---|
1717 | #Alternatively from elsewhere |
---|
1718 | #or, as a last resort, |
---|
1719 | #specified here. |
---|
1720 | #The value of -100 will work |
---|
1721 | #for the Wollongong tsunami |
---|
1722 | #scenario but is very hacky |
---|
1723 | """Convert 'Ferret' NetCDF format for wave propagation to |
---|
1724 | sww format native to pyvolution. |
---|
1725 | |
---|
1726 | Specify only basename_in and read files of the form |
---|
1727 | basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing |
---|
1728 | relative height, x-velocity and y-velocity, respectively. |
---|
1729 | |
---|
1730 | Also convert latitude and longitude to UTM. All coordinates are |
---|
1731 | assumed to be given in the GDA94 datum. |
---|
1732 | |
---|
1733 | min's and max's: If omitted - full extend is used. |
---|
1734 | To include a value min may equal it, while max must exceed it. |
---|
1735 | Lat and lon are assuemd to be in decimal degrees |
---|
1736 | |
---|
1737 | origin is a 3-tuple with geo referenced |
---|
1738 | UTM coordinates (zone, easting, northing) |
---|
1739 | |
---|
1740 | nc format has values organised as HA[TIME, LATITUDE, LONGITUDE] |
---|
1741 | which means that longitude is the fastest |
---|
1742 | varying dimension (row major order, so to speak) |
---|
1743 | |
---|
1744 | ferret2sww uses grid points as vertices in a triangular grid |
---|
1745 | counting vertices from lower left corner upwards, then right |
---|
1746 | """ |
---|
1747 | |
---|
1748 | import os |
---|
1749 | from Scientific.IO.NetCDF import NetCDFFile |
---|
1750 | from Numeric import Float, Int, Int32, searchsorted, zeros, array |
---|
1751 | from Numeric import allclose, around |
---|
1752 | |
---|
1753 | precision = Float |
---|
1754 | |
---|
1755 | |
---|
1756 | #Get NetCDF data |
---|
1757 | if verbose: print 'Reading files %s_*.nc' %basename_in |
---|
1758 | file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm) |
---|
1759 | file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s) |
---|
1760 | file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s) |
---|
1761 | file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m) |
---|
1762 | |
---|
1763 | if basename_out is None: |
---|
1764 | swwname = basename_in + '.sww' |
---|
1765 | else: |
---|
1766 | swwname = basename_out + '.sww' |
---|
1767 | |
---|
1768 | times = file_h.variables['TIME'] |
---|
1769 | latitudes = file_h.variables['LAT'] |
---|
1770 | longitudes = file_h.variables['LON'] |
---|
1771 | |
---|
1772 | |
---|
1773 | |
---|
1774 | #Precision used by most for lat/lon is 4 or 5 decimals |
---|
1775 | e_lat = around(file_e.variables['LAT'][:], 5) |
---|
1776 | e_lon = around(file_e.variables['LON'][:], 5) |
---|
1777 | |
---|
1778 | #Check that files are compatible |
---|
1779 | assert allclose(latitudes, file_u.variables['LAT']) |
---|
1780 | assert allclose(latitudes, file_v.variables['LAT']) |
---|
1781 | assert allclose(latitudes, e_lat) |
---|
1782 | |
---|
1783 | assert allclose(longitudes, file_u.variables['LON']) |
---|
1784 | assert allclose(longitudes, file_v.variables['LON']) |
---|
1785 | assert allclose(longitudes, e_lon) |
---|
1786 | |
---|
1787 | |
---|
1788 | |
---|
1789 | if mint == None: |
---|
1790 | jmin = 0 |
---|
1791 | else: |
---|
1792 | jmin = searchsorted(times, mint) |
---|
1793 | |
---|
1794 | if maxt == None: |
---|
1795 | jmax=len(times) |
---|
1796 | else: |
---|
1797 | jmax = searchsorted(times, maxt) |
---|
1798 | |
---|
1799 | if minlat == None: |
---|
1800 | kmin=0 |
---|
1801 | else: |
---|
1802 | kmin = searchsorted(latitudes, minlat) |
---|
1803 | |
---|
1804 | if maxlat == None: |
---|
1805 | kmax = len(latitudes) |
---|
1806 | else: |
---|
1807 | kmax = searchsorted(latitudes, maxlat) |
---|
1808 | |
---|
1809 | if minlon == None: |
---|
1810 | lmin=0 |
---|
1811 | else: |
---|
1812 | lmin = searchsorted(longitudes, minlon) |
---|
1813 | |
---|
1814 | if maxlon == None: |
---|
1815 | lmax = len(longitudes) |
---|
1816 | else: |
---|
1817 | lmax = searchsorted(longitudes, maxlon) |
---|
1818 | |
---|
1819 | |
---|
1820 | |
---|
1821 | times = times[jmin:jmax] |
---|
1822 | latitudes = latitudes[kmin:kmax] |
---|
1823 | longitudes = longitudes[lmin:lmax] |
---|
1824 | |
---|
1825 | |
---|
1826 | if verbose: print 'cropping' |
---|
1827 | zname = 'ELEVATION' |
---|
1828 | |
---|
1829 | |
---|
1830 | amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax] |
---|
1831 | uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon |
---|
1832 | vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat |
---|
1833 | elevations = file_e.variables[zname][kmin:kmax, lmin:lmax] |
---|
1834 | |
---|
1835 | # if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]: |
---|
1836 | # elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax] |
---|
1837 | # elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]: |
---|
1838 | # from Numeric import asarray |
---|
1839 | # elevations=elevations.tolist() |
---|
1840 | # elevations.reverse() |
---|
1841 | # elevations=asarray(elevations) |
---|
1842 | # else: |
---|
1843 | # from Numeric import asarray |
---|
1844 | # elevations=elevations.tolist() |
---|
1845 | # elevations.reverse() |
---|
1846 | # elevations=asarray(elevations) |
---|
1847 | # 'print hmmm' |
---|
1848 | |
---|
1849 | |
---|
1850 | |
---|
1851 | #Get missing values |
---|
1852 | nan_ha = file_h.variables['HA'].missing_value[0] |
---|
1853 | nan_ua = file_u.variables['UA'].missing_value[0] |
---|
1854 | nan_va = file_v.variables['VA'].missing_value[0] |
---|
1855 | if hasattr(file_e.variables[zname],'missing_value'): |
---|
1856 | nan_e = file_e.variables[zname].missing_value[0] |
---|
1857 | else: |
---|
1858 | nan_e = None |
---|
1859 | |
---|
1860 | #Cleanup |
---|
1861 | from Numeric import sometrue |
---|
1862 | |
---|
1863 | missing = (amplitudes == nan_ha) |
---|
1864 | if sometrue (missing): |
---|
1865 | if fail_on_NaN: |
---|
1866 | msg = 'NetCDFFile %s contains missing values'\ |
---|
1867 | %(basename_in+'_ha.nc') |
---|
1868 | raise msg |
---|
1869 | else: |
---|
1870 | amplitudes = amplitudes*(missing==0) + missing*NaN_filler |
---|
1871 | |
---|
1872 | missing = (uspeed == nan_ua) |
---|
1873 | if sometrue (missing): |
---|
1874 | if fail_on_NaN: |
---|
1875 | msg = 'NetCDFFile %s contains missing values'\ |
---|
1876 | %(basename_in+'_ua.nc') |
---|
1877 | raise msg |
---|
1878 | else: |
---|
1879 | uspeed = uspeed*(missing==0) + missing*NaN_filler |
---|
1880 | |
---|
1881 | missing = (vspeed == nan_va) |
---|
1882 | if sometrue (missing): |
---|
1883 | if fail_on_NaN: |
---|
1884 | msg = 'NetCDFFile %s contains missing values'\ |
---|
1885 | %(basename_in+'_va.nc') |
---|
1886 | raise msg |
---|
1887 | else: |
---|
1888 | vspeed = vspeed*(missing==0) + missing*NaN_filler |
---|
1889 | |
---|
1890 | |
---|
1891 | missing = (elevations == nan_e) |
---|
1892 | if sometrue (missing): |
---|
1893 | if fail_on_NaN: |
---|
1894 | msg = 'NetCDFFile %s contains missing values'\ |
---|
1895 | %(basename_in+'_e.nc') |
---|
1896 | raise msg |
---|
1897 | else: |
---|
1898 | elevations = elevations*(missing==0) + missing*NaN_filler |
---|
1899 | |
---|
1900 | ####### |
---|
1901 | |
---|
1902 | |
---|
1903 | |
---|
1904 | number_of_times = times.shape[0] |
---|
1905 | number_of_latitudes = latitudes.shape[0] |
---|
1906 | number_of_longitudes = longitudes.shape[0] |
---|
1907 | |
---|
1908 | assert amplitudes.shape[0] == number_of_times |
---|
1909 | assert amplitudes.shape[1] == number_of_latitudes |
---|
1910 | assert amplitudes.shape[2] == number_of_longitudes |
---|
1911 | |
---|
1912 | if verbose: |
---|
1913 | print '------------------------------------------------' |
---|
1914 | print 'Statistics:' |
---|
1915 | print ' Extent (lat/lon):' |
---|
1916 | print ' lat in [%f, %f], len(lat) == %d'\ |
---|
1917 | %(min(latitudes.flat), max(latitudes.flat), |
---|
1918 | len(latitudes.flat)) |
---|
1919 | print ' lon in [%f, %f], len(lon) == %d'\ |
---|
1920 | %(min(longitudes.flat), max(longitudes.flat), |
---|
1921 | len(longitudes.flat)) |
---|
1922 | print ' t in [%f, %f], len(t) == %d'\ |
---|
1923 | %(min(times.flat), max(times.flat), len(times.flat)) |
---|
1924 | |
---|
1925 | q = amplitudes.flat |
---|
1926 | name = 'Amplitudes (ha) [cm]' |
---|
1927 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
1928 | |
---|
1929 | q = uspeed.flat |
---|
1930 | name = 'Speeds (ua) [cm/s]' |
---|
1931 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
1932 | |
---|
1933 | q = vspeed.flat |
---|
1934 | name = 'Speeds (va) [cm/s]' |
---|
1935 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
1936 | |
---|
1937 | q = elevations.flat |
---|
1938 | name = 'Elevations (e) [m]' |
---|
1939 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
1940 | |
---|
1941 | |
---|
1942 | #print number_of_latitudes, number_of_longitudes |
---|
1943 | number_of_points = number_of_latitudes*number_of_longitudes |
---|
1944 | number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2 |
---|
1945 | |
---|
1946 | |
---|
1947 | file_h.close() |
---|
1948 | file_u.close() |
---|
1949 | file_v.close() |
---|
1950 | file_e.close() |
---|
1951 | |
---|
1952 | |
---|
1953 | # NetCDF file definition |
---|
1954 | outfile = NetCDFFile(swwname, 'w') |
---|
1955 | |
---|
1956 | #Create new file |
---|
1957 | outfile.institution = 'Geoscience Australia' |
---|
1958 | outfile.description = 'Converted from Ferret files: %s, %s, %s, %s'\ |
---|
1959 | %(basename_in + '_ha.nc', |
---|
1960 | basename_in + '_ua.nc', |
---|
1961 | basename_in + '_va.nc', |
---|
1962 | basename_in + '_e.nc') |
---|
1963 | |
---|
1964 | |
---|
1965 | #For sww compatibility |
---|
1966 | outfile.smoothing = 'Yes' |
---|
1967 | outfile.order = 1 |
---|
1968 | |
---|
1969 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
1970 | outfile.starttime = starttime = times[0] |
---|
1971 | times = times - starttime #Store relative times |
---|
1972 | |
---|
1973 | # dimension definitions |
---|
1974 | outfile.createDimension('number_of_volumes', number_of_volumes) |
---|
1975 | |
---|
1976 | outfile.createDimension('number_of_vertices', 3) |
---|
1977 | outfile.createDimension('number_of_points', number_of_points) |
---|
1978 | |
---|
1979 | |
---|
1980 | #outfile.createDimension('number_of_timesteps', len(times)) |
---|
1981 | outfile.createDimension('number_of_timesteps', len(times)) |
---|
1982 | |
---|
1983 | # variable definitions |
---|
1984 | outfile.createVariable('x', precision, ('number_of_points',)) |
---|
1985 | outfile.createVariable('y', precision, ('number_of_points',)) |
---|
1986 | outfile.createVariable('elevation', precision, ('number_of_points',)) |
---|
1987 | |
---|
1988 | #FIXME: Backwards compatibility |
---|
1989 | outfile.createVariable('z', precision, ('number_of_points',)) |
---|
1990 | ################################# |
---|
1991 | |
---|
1992 | outfile.createVariable('volumes', Int, ('number_of_volumes', |
---|
1993 | 'number_of_vertices')) |
---|
1994 | |
---|
1995 | outfile.createVariable('time', precision, |
---|
1996 | ('number_of_timesteps',)) |
---|
1997 | |
---|
1998 | outfile.createVariable('stage', precision, |
---|
1999 | ('number_of_timesteps', |
---|
2000 | 'number_of_points')) |
---|
2001 | |
---|
2002 | outfile.createVariable('xmomentum', precision, |
---|
2003 | ('number_of_timesteps', |
---|
2004 | 'number_of_points')) |
---|
2005 | |
---|
2006 | outfile.createVariable('ymomentum', precision, |
---|
2007 | ('number_of_timesteps', |
---|
2008 | 'number_of_points')) |
---|
2009 | |
---|
2010 | |
---|
2011 | #Store |
---|
2012 | from coordinate_transforms.redfearn import redfearn |
---|
2013 | x = zeros(number_of_points, Float) #Easting |
---|
2014 | y = zeros(number_of_points, Float) #Northing |
---|
2015 | |
---|
2016 | |
---|
2017 | if verbose: print 'Making triangular grid' |
---|
2018 | #Check zone boundaries |
---|
2019 | refzone, _, _ = redfearn(latitudes[0],longitudes[0]) |
---|
2020 | |
---|
2021 | vertices = {} |
---|
2022 | i = 0 |
---|
2023 | for k, lat in enumerate(latitudes): #Y direction |
---|
2024 | for l, lon in enumerate(longitudes): #X direction |
---|
2025 | |
---|
2026 | vertices[l,k] = i |
---|
2027 | |
---|
2028 | zone, easting, northing = redfearn(lat,lon) |
---|
2029 | |
---|
2030 | msg = 'Zone boundary crossed at longitude =', lon |
---|
2031 | #assert zone == refzone, msg |
---|
2032 | #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing) |
---|
2033 | x[i] = easting |
---|
2034 | y[i] = northing |
---|
2035 | i += 1 |
---|
2036 | |
---|
2037 | |
---|
2038 | #Construct 2 triangles per 'rectangular' element |
---|
2039 | volumes = [] |
---|
2040 | for l in range(number_of_longitudes-1): #X direction |
---|
2041 | for k in range(number_of_latitudes-1): #Y direction |
---|
2042 | v1 = vertices[l,k+1] |
---|
2043 | v2 = vertices[l,k] |
---|
2044 | v3 = vertices[l+1,k+1] |
---|
2045 | v4 = vertices[l+1,k] |
---|
2046 | |
---|
2047 | volumes.append([v1,v2,v3]) #Upper element |
---|
2048 | volumes.append([v4,v3,v2]) #Lower element |
---|
2049 | |
---|
2050 | volumes = array(volumes) |
---|
2051 | |
---|
2052 | if origin == None: |
---|
2053 | zone = refzone |
---|
2054 | xllcorner = min(x) |
---|
2055 | yllcorner = min(y) |
---|
2056 | else: |
---|
2057 | zone = origin[0] |
---|
2058 | xllcorner = origin[1] |
---|
2059 | yllcorner = origin[2] |
---|
2060 | |
---|
2061 | |
---|
2062 | outfile.xllcorner = xllcorner |
---|
2063 | outfile.yllcorner = yllcorner |
---|
2064 | outfile.zone = zone |
---|
2065 | |
---|
2066 | |
---|
2067 | if elevation is not None: |
---|
2068 | z = elevation |
---|
2069 | else: |
---|
2070 | if inverted_bathymetry: |
---|
2071 | z = -1*elevations |
---|
2072 | else: |
---|
2073 | z = elevations |
---|
2074 | #FIXME: z should be obtained from MOST and passed in here |
---|
2075 | |
---|
2076 | from Numeric import resize |
---|
2077 | z = resize(z,outfile.variables['z'][:].shape) |
---|
2078 | outfile.variables['x'][:] = x - xllcorner |
---|
2079 | outfile.variables['y'][:] = y - yllcorner |
---|
2080 | outfile.variables['z'][:] = z |
---|
2081 | outfile.variables['elevation'][:] = z #FIXME HACK |
---|
2082 | outfile.variables['time'][:] = times #Store time relative |
---|
2083 | outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64 |
---|
2084 | |
---|
2085 | |
---|
2086 | |
---|
2087 | #Time stepping |
---|
2088 | stage = outfile.variables['stage'] |
---|
2089 | xmomentum = outfile.variables['xmomentum'] |
---|
2090 | ymomentum = outfile.variables['ymomentum'] |
---|
2091 | |
---|
2092 | if verbose: print 'Converting quantities' |
---|
2093 | n = len(times) |
---|
2094 | for j in range(n): |
---|
2095 | if verbose and j%((n+10)/10)==0: print ' Doing %d of %d' %(j, n) |
---|
2096 | i = 0 |
---|
2097 | for k in range(number_of_latitudes): #Y direction |
---|
2098 | for l in range(number_of_longitudes): #X direction |
---|
2099 | w = zscale*amplitudes[j,k,l]/100 + mean_stage |
---|
2100 | stage[j,i] = w |
---|
2101 | h = w - z[i] |
---|
2102 | xmomentum[j,i] = uspeed[j,k,l]/100*h |
---|
2103 | ymomentum[j,i] = vspeed[j,k,l]/100*h |
---|
2104 | i += 1 |
---|
2105 | |
---|
2106 | |
---|
2107 | if verbose: |
---|
2108 | x = outfile.variables['x'][:] |
---|
2109 | y = outfile.variables['y'][:] |
---|
2110 | print '------------------------------------------------' |
---|
2111 | print 'Statistics of output file:' |
---|
2112 | print ' Name: %s' %swwname |
---|
2113 | print ' Reference:' |
---|
2114 | print ' Lower left corner: [%f, %f]'\ |
---|
2115 | %(xllcorner, yllcorner) |
---|
2116 | print ' Start time: %f' %starttime |
---|
2117 | print ' Extent:' |
---|
2118 | print ' x [m] in [%f, %f], len(x) == %d'\ |
---|
2119 | %(min(x.flat), max(x.flat), len(x.flat)) |
---|
2120 | print ' y [m] in [%f, %f], len(y) == %d'\ |
---|
2121 | %(min(y.flat), max(y.flat), len(y.flat)) |
---|
2122 | print ' t [s] in [%f, %f], len(t) == %d'\ |
---|
2123 | %(min(times), max(times), len(times)) |
---|
2124 | print ' Quantities [SI units]:' |
---|
2125 | for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']: |
---|
2126 | q = outfile.variables[name][:].flat |
---|
2127 | print ' %s in [%f, %f]' %(name, min(q), max(q)) |
---|
2128 | |
---|
2129 | |
---|
2130 | |
---|
2131 | |
---|
2132 | outfile.close() |
---|
2133 | |
---|
2134 | |
---|
2135 | |
---|
2136 | |
---|
2137 | def timefile2netcdf(filename, quantity_names = None): |
---|
2138 | """Template for converting typical text files with time series to |
---|
2139 | NetCDF tms file. |
---|
2140 | |
---|
2141 | |
---|
2142 | The file format is assumed to be either two fields separated by a comma: |
---|
2143 | |
---|
2144 | time [DD/MM/YY hh:mm:ss], value0 value1 value2 ... |
---|
2145 | |
---|
2146 | E.g |
---|
2147 | |
---|
2148 | 31/08/04 00:00:00, 1.328223 0 0 |
---|
2149 | 31/08/04 00:15:00, 1.292912 0 0 |
---|
2150 | |
---|
2151 | will provide a time dependent function f(t) with three attributes |
---|
2152 | |
---|
2153 | filename is assumed to be the rootname with extenisons .txt and .sww |
---|
2154 | """ |
---|
2155 | |
---|
2156 | import time, calendar |
---|
2157 | from Numeric import array |
---|
2158 | from config import time_format |
---|
2159 | from util import ensure_numeric |
---|
2160 | |
---|
2161 | |
---|
2162 | fid = open(filename + '.txt') |
---|
2163 | line = fid.readline() |
---|
2164 | fid.close() |
---|
2165 | |
---|
2166 | fields = line.split(',') |
---|
2167 | msg = 'File %s must have the format date, value0 value1 value2 ...' |
---|
2168 | assert len(fields) == 2, msg |
---|
2169 | |
---|
2170 | try: |
---|
2171 | starttime = calendar.timegm(time.strptime(fields[0], time_format)) |
---|
2172 | except ValueError: |
---|
2173 | msg = 'First field in file %s must be' %filename |
---|
2174 | msg += ' date-time with format %s.\n' %time_format |
---|
2175 | msg += 'I got %s instead.' %fields[0] |
---|
2176 | raise msg |
---|
2177 | |
---|
2178 | |
---|
2179 | #Split values |
---|
2180 | values = [] |
---|
2181 | for value in fields[1].split(): |
---|
2182 | values.append(float(value)) |
---|
2183 | |
---|
2184 | q = ensure_numeric(values) |
---|
2185 | |
---|
2186 | msg = 'ERROR: File must contain at least one independent value' |
---|
2187 | assert len(q.shape) == 1, msg |
---|
2188 | |
---|
2189 | |
---|
2190 | |
---|
2191 | #Read times proper |
---|
2192 | from Numeric import zeros, Float, alltrue |
---|
2193 | from config import time_format |
---|
2194 | import time, calendar |
---|
2195 | |
---|
2196 | fid = open(filename + '.txt') |
---|
2197 | lines = fid.readlines() |
---|
2198 | fid.close() |
---|
2199 | |
---|
2200 | N = len(lines) |
---|
2201 | d = len(q) |
---|
2202 | |
---|
2203 | T = zeros(N, Float) #Time |
---|
2204 | Q = zeros((N, d), Float) #Values |
---|
2205 | |
---|
2206 | for i, line in enumerate(lines): |
---|
2207 | fields = line.split(',') |
---|
2208 | realtime = calendar.timegm(time.strptime(fields[0], time_format)) |
---|
2209 | |
---|
2210 | T[i] = realtime - starttime |
---|
2211 | |
---|
2212 | for j, value in enumerate(fields[1].split()): |
---|
2213 | Q[i, j] = float(value) |
---|
2214 | |
---|
2215 | msg = 'File %s must list time as a monotonuosly ' %filename |
---|
2216 | msg += 'increasing sequence' |
---|
2217 | assert alltrue( T[1:] - T[:-1] > 0 ), msg |
---|
2218 | |
---|
2219 | |
---|
2220 | #Create NetCDF file |
---|
2221 | from Scientific.IO.NetCDF import NetCDFFile |
---|
2222 | |
---|
2223 | fid = NetCDFFile(filename + '.tms', 'w') |
---|
2224 | |
---|
2225 | |
---|
2226 | fid.institution = 'Geoscience Australia' |
---|
2227 | fid.description = 'Time series' |
---|
2228 | |
---|
2229 | |
---|
2230 | #Reference point |
---|
2231 | #Start time in seconds since the epoch (midnight 1/1/1970) |
---|
2232 | #FIXME: Use Georef |
---|
2233 | fid.starttime = starttime |
---|
2234 | |
---|
2235 | |
---|
2236 | # dimension definitions |
---|
2237 | #fid.createDimension('number_of_volumes', self.number_of_volumes) |
---|
2238 | #fid.createDimension('number_of_vertices', 3) |
---|
2239 | |
---|
2240 | |
---|
2241 | fid.createDimension('number_of_timesteps', len(T)) |
---|
2242 | |
---|
2243 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
---|
2244 | |
---|
2245 | fid.variables['time'][:] = T |
---|
2246 | |
---|
2247 | for i in range(Q.shape[1]): |
---|
2248 | try: |
---|
2249 | name = quantity_names[i] |
---|
2250 | except: |
---|
2251 | name = 'Attribute%d'%i |
---|
2252 | |
---|
2253 | fid.createVariable(name, Float, ('number_of_timesteps',)) |
---|
2254 | fid.variables[name][:] = Q[:,i] |
---|
2255 | |
---|
2256 | fid.close() |
---|
2257 | |
---|
2258 | |
---|
2259 | def extent_sww(file_name): |
---|
2260 | """ |
---|
2261 | Read in an sww file. |
---|
2262 | |
---|
2263 | Input; |
---|
2264 | file_name - the sww file |
---|
2265 | |
---|
2266 | Output; |
---|
2267 | z - Vector of bed elevation |
---|
2268 | volumes - Array. Each row has 3 values, representing |
---|
2269 | the vertices that define the volume |
---|
2270 | time - Vector of the times where there is stage information |
---|
2271 | stage - array with respect to time and vertices (x,y) |
---|
2272 | """ |
---|
2273 | |
---|
2274 | |
---|
2275 | from Scientific.IO.NetCDF import NetCDFFile |
---|
2276 | |
---|
2277 | #Check contents |
---|
2278 | #Get NetCDF |
---|
2279 | fid = NetCDFFile(file_name, 'r') |
---|
2280 | |
---|
2281 | # Get the variables |
---|
2282 | x = fid.variables['x'][:] |
---|
2283 | y = fid.variables['y'][:] |
---|
2284 | stage = fid.variables['stage'][:] |
---|
2285 | #print "stage",stage |
---|
2286 | #print "stage.shap",stage.shape |
---|
2287 | #print "min(stage.flat), mpythonax(stage.flat)",min(stage.flat), max(stage.flat) |
---|
2288 | #print "min(stage)",min(stage) |
---|
2289 | |
---|
2290 | fid.close() |
---|
2291 | |
---|
2292 | return [min(x),max(x),min(y),max(y),min(stage.flat),max(stage.flat)] |
---|
2293 | |
---|
2294 | |
---|
2295 | def sww2domain(filename,boundary=None,t=None,\ |
---|
2296 | fail_if_NaN=True,NaN_filler=0\ |
---|
2297 | ,verbose = True,very_verbose = False): |
---|
2298 | """ |
---|
2299 | Usage: domain = sww2domain('file.sww',t=time (default = last time in file)) |
---|
2300 | |
---|
2301 | Boundary is not recommended if domian.smooth is not selected, as it |
---|
2302 | uses unique coordinates, but not unique boundaries. This means that |
---|
2303 | the boundary file will not be compatable with the coordiantes, and will |
---|
2304 | give a different final boundary, or crash. |
---|
2305 | """ |
---|
2306 | NaN=9.969209968386869e+036 |
---|
2307 | #initialise NaN. |
---|
2308 | |
---|
2309 | from Scientific.IO.NetCDF import NetCDFFile |
---|
2310 | from shallow_water import Domain |
---|
2311 | from Numeric import asarray, transpose, resize |
---|
2312 | |
---|
2313 | if verbose: print 'Reading from ', filename |
---|
2314 | fid = NetCDFFile(filename, 'r') #Open existing file for read |
---|
2315 | time = fid.variables['time'] #Timesteps |
---|
2316 | if t is None: |
---|
2317 | t = time[-1] |
---|
2318 | time_interp = get_time_interp(time,t) |
---|
2319 | |
---|
2320 | # Get the variables as Numeric arrays |
---|
2321 | x = fid.variables['x'][:] #x-coordinates of vertices |
---|
2322 | y = fid.variables['y'][:] #y-coordinates of vertices |
---|
2323 | elevation = fid.variables['elevation'] #Elevation |
---|
2324 | stage = fid.variables['stage'] #Water level |
---|
2325 | xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction |
---|
2326 | ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction |
---|
2327 | |
---|
2328 | starttime = fid.starttime[0] |
---|
2329 | volumes = fid.variables['volumes'][:] #Connectivity |
---|
2330 | coordinates=transpose(asarray([x.tolist(),y.tolist()])) |
---|
2331 | |
---|
2332 | conserved_quantities = [] |
---|
2333 | interpolated_quantities = {} |
---|
2334 | other_quantities = [] |
---|
2335 | |
---|
2336 | # get geo_reference |
---|
2337 | #sww files don't have to have a geo_ref |
---|
2338 | try: |
---|
2339 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
2340 | except: #AttributeError, e: |
---|
2341 | geo_reference = None |
---|
2342 | |
---|
2343 | if verbose: print ' getting quantities' |
---|
2344 | for quantity in fid.variables.keys(): |
---|
2345 | dimensions = fid.variables[quantity].dimensions |
---|
2346 | if 'number_of_timesteps' in dimensions: |
---|
2347 | conserved_quantities.append(quantity) |
---|
2348 | interpolated_quantities[quantity]=\ |
---|
2349 | interpolated_quantity(fid.variables[quantity][:],time_interp) |
---|
2350 | else: other_quantities.append(quantity) |
---|
2351 | |
---|
2352 | other_quantities.remove('x') |
---|
2353 | other_quantities.remove('y') |
---|
2354 | other_quantities.remove('z') |
---|
2355 | other_quantities.remove('volumes') |
---|
2356 | |
---|
2357 | conserved_quantities.remove('time') |
---|
2358 | |
---|
2359 | if verbose: print ' building domain' |
---|
2360 | # From domain.Domain: |
---|
2361 | # domain = Domain(coordinates, volumes,\ |
---|
2362 | # conserved_quantities = conserved_quantities,\ |
---|
2363 | # other_quantities = other_quantities,zone=zone,\ |
---|
2364 | # xllcorner=xllcorner, yllcorner=yllcorner) |
---|
2365 | |
---|
2366 | # From shallow_water.Domain: |
---|
2367 | coordinates=coordinates.tolist() |
---|
2368 | volumes=volumes.tolist() |
---|
2369 | #FIXME:should this be in mesh?(peter row) |
---|
2370 | if fid.smoothing == 'Yes': unique = False |
---|
2371 | else: unique = True |
---|
2372 | if unique: |
---|
2373 | coordinates,volumes,boundary=weed(coordinates,volumes,boundary) |
---|
2374 | |
---|
2375 | |
---|
2376 | domain = Domain(coordinates, volumes, boundary) |
---|
2377 | |
---|
2378 | if not boundary is None: |
---|
2379 | domain.boundary = boundary |
---|
2380 | |
---|
2381 | domain.geo_reference = geo_reference |
---|
2382 | |
---|
2383 | domain.starttime=float(starttime)+float(t) |
---|
2384 | domain.time=0.0 |
---|
2385 | |
---|
2386 | for quantity in other_quantities: |
---|
2387 | try: |
---|
2388 | NaN = fid.variables[quantity].missing_value |
---|
2389 | except: |
---|
2390 | pass #quantity has no missing_value number |
---|
2391 | X = fid.variables[quantity][:] |
---|
2392 | if very_verbose: |
---|
2393 | print ' ',quantity |
---|
2394 | print ' NaN =',NaN |
---|
2395 | print ' max(X)' |
---|
2396 | print ' ',max(X) |
---|
2397 | print ' max(X)==NaN' |
---|
2398 | print ' ',max(X)==NaN |
---|
2399 | print '' |
---|
2400 | if (max(X)==NaN) or (min(X)==NaN): |
---|
2401 | if fail_if_NaN: |
---|
2402 | msg = 'quantity "%s" contains no_data entry'%quantity |
---|
2403 | raise msg |
---|
2404 | else: |
---|
2405 | data = (X<>NaN) |
---|
2406 | X = (X*data)+(data==0)*NaN_filler |
---|
2407 | if unique: |
---|
2408 | X = resize(X,(len(X)/3,3)) |
---|
2409 | domain.set_quantity(quantity,X) |
---|
2410 | # |
---|
2411 | for quantity in conserved_quantities: |
---|
2412 | try: |
---|
2413 | NaN = fid.variables[quantity].missing_value |
---|
2414 | except: |
---|
2415 | pass #quantity has no missing_value number |
---|
2416 | X = interpolated_quantities[quantity] |
---|
2417 | if very_verbose: |
---|
2418 | print ' ',quantity |
---|
2419 | print ' NaN =',NaN |
---|
2420 | print ' max(X)' |
---|
2421 | print ' ',max(X) |
---|
2422 | print ' max(X)==NaN' |
---|
2423 | print ' ',max(X)==NaN |
---|
2424 | print '' |
---|
2425 | if (max(X)==NaN) or (min(X)==NaN): |
---|
2426 | if fail_if_NaN: |
---|
2427 | msg = 'quantity "%s" contains no_data entry'%quantity |
---|
2428 | raise msg |
---|
2429 | else: |
---|
2430 | data = (X<>NaN) |
---|
2431 | X = (X*data)+(data==0)*NaN_filler |
---|
2432 | if unique: |
---|
2433 | X = resize(X,(X.shape[0]/3,3)) |
---|
2434 | domain.set_quantity(quantity,X) |
---|
2435 | fid.close() |
---|
2436 | return domain |
---|
2437 | |
---|
2438 | def interpolated_quantity(saved_quantity,time_interp): |
---|
2439 | |
---|
2440 | #given an index and ratio, interpolate quantity with respect to time. |
---|
2441 | index,ratio = time_interp |
---|
2442 | Q = saved_quantity |
---|
2443 | if ratio > 0: |
---|
2444 | q = (1-ratio)*Q[index]+ ratio*Q[index+1] |
---|
2445 | else: |
---|
2446 | q = Q[index] |
---|
2447 | #Return vector of interpolated values |
---|
2448 | return q |
---|
2449 | |
---|
2450 | def get_time_interp(time,t=None): |
---|
2451 | #Finds the ratio and index for time interpolation. |
---|
2452 | #It is borrowed from previous pyvolution code. |
---|
2453 | if t is None: |
---|
2454 | t=time[-1] |
---|
2455 | index = -1 |
---|
2456 | ratio = 0. |
---|
2457 | else: |
---|
2458 | T = time |
---|
2459 | tau = t |
---|
2460 | index=0 |
---|
2461 | msg = 'Time interval derived from file %s [%s:%s]'\ |
---|
2462 | %('FIXMEfilename', T[0], T[-1]) |
---|
2463 | msg += ' does not match model time: %s' %tau |
---|
2464 | if tau < time[0]: raise msg |
---|
2465 | if tau > time[-1]: raise msg |
---|
2466 | while tau > time[index]: index += 1 |
---|
2467 | while tau < time[index]: index -= 1 |
---|
2468 | if tau == time[index]: |
---|
2469 | #Protect against case where tau == time[-1] (last time) |
---|
2470 | # - also works in general when tau == time[i] |
---|
2471 | ratio = 0 |
---|
2472 | else: |
---|
2473 | #t is now between index and index+1 |
---|
2474 | ratio = (tau - time[index])/(time[index+1] - time[index]) |
---|
2475 | return (index,ratio) |
---|
2476 | |
---|
2477 | |
---|
2478 | def weed(coordinates,volumes,boundary = None): |
---|
2479 | if type(coordinates)=='array': |
---|
2480 | coordinates = coordinates.tolist() |
---|
2481 | if type(volumes)=='array': |
---|
2482 | volumes = volumes.tolist() |
---|
2483 | |
---|
2484 | unique = False |
---|
2485 | point_dict = {} |
---|
2486 | same_point = {} |
---|
2487 | for i in range(len(coordinates)): |
---|
2488 | point = tuple(coordinates[i]) |
---|
2489 | if point_dict.has_key(point): |
---|
2490 | unique = True |
---|
2491 | same_point[i]=point |
---|
2492 | #to change all point i references to point j |
---|
2493 | else: |
---|
2494 | point_dict[point]=i |
---|
2495 | same_point[i]=point |
---|
2496 | |
---|
2497 | coordinates = [] |
---|
2498 | i = 0 |
---|
2499 | for point in point_dict.keys(): |
---|
2500 | point = tuple(point) |
---|
2501 | coordinates.append(list(point)) |
---|
2502 | point_dict[point]=i |
---|
2503 | i+=1 |
---|
2504 | |
---|
2505 | |
---|
2506 | for volume in volumes: |
---|
2507 | for i in range(len(volume)): |
---|
2508 | index = volume[i] |
---|
2509 | if index>-1: |
---|
2510 | volume[i]=point_dict[same_point[index]] |
---|
2511 | |
---|
2512 | new_boundary = {} |
---|
2513 | if not boundary is None: |
---|
2514 | for segment in boundary.keys(): |
---|
2515 | point0 = point_dict[same_point[segment[0]]] |
---|
2516 | point1 = point_dict[same_point[segment[1]]] |
---|
2517 | label = boundary[segment] |
---|
2518 | #FIXME should the bounday attributes be concaterated |
---|
2519 | #('exterior, pond') or replaced ('pond')(peter row) |
---|
2520 | |
---|
2521 | if new_boundary.has_key((point0,point1)): |
---|
2522 | new_boundary[(point0,point1)]=new_boundary[(point0,point1)]#\ |
---|
2523 | #+','+label |
---|
2524 | |
---|
2525 | elif new_boundary.has_key((point1,point0)): |
---|
2526 | new_boundary[(point1,point0)]=new_boundary[(point1,point0)]#\ |
---|
2527 | #+','+label |
---|
2528 | else: new_boundary[(point0,point1)]=label |
---|
2529 | |
---|
2530 | boundary = new_boundary |
---|
2531 | |
---|
2532 | return coordinates,volumes,boundary |
---|
2533 | |
---|
2534 | |
---|
2535 | |
---|
2536 | |
---|
2537 | |
---|
2538 | |
---|
2539 | |
---|
2540 | |
---|
2541 | |
---|
2542 | |
---|
2543 | ############################################### |
---|
2544 | #OBSOLETE STUFF |
---|
2545 | #Native checkpoint format. |
---|
2546 | #Information needed to recreate a state is preserved |
---|
2547 | #FIXME: Rethink and maybe use netcdf format |
---|
2548 | def cpt_variable_writer(filename, t, v0, v1, v2): |
---|
2549 | """Store all conserved quantities to file |
---|
2550 | """ |
---|
2551 | |
---|
2552 | M, N = v0.shape |
---|
2553 | |
---|
2554 | FN = create_filename(filename, 'cpt', M, t) |
---|
2555 | #print 'Writing to %s' %FN |
---|
2556 | |
---|
2557 | fid = open(FN, 'w') |
---|
2558 | for i in range(M): |
---|
2559 | for j in range(N): |
---|
2560 | fid.write('%.16e ' %v0[i,j]) |
---|
2561 | for j in range(N): |
---|
2562 | fid.write('%.16e ' %v1[i,j]) |
---|
2563 | for j in range(N): |
---|
2564 | fid.write('%.16e ' %v2[i,j]) |
---|
2565 | |
---|
2566 | fid.write('\n') |
---|
2567 | fid.close() |
---|
2568 | |
---|
2569 | |
---|
2570 | def cpt_variable_reader(filename, t, v0, v1, v2): |
---|
2571 | """Store all conserved quantities to file |
---|
2572 | """ |
---|
2573 | |
---|
2574 | M, N = v0.shape |
---|
2575 | |
---|
2576 | FN = create_filename(filename, 'cpt', M, t) |
---|
2577 | #print 'Reading from %s' %FN |
---|
2578 | |
---|
2579 | fid = open(FN) |
---|
2580 | |
---|
2581 | |
---|
2582 | for i in range(M): |
---|
2583 | values = fid.readline().split() #Get one line |
---|
2584 | |
---|
2585 | for j in range(N): |
---|
2586 | v0[i,j] = float(values[j]) |
---|
2587 | v1[i,j] = float(values[3+j]) |
---|
2588 | v2[i,j] = float(values[6+j]) |
---|
2589 | |
---|
2590 | fid.close() |
---|
2591 | |
---|
2592 | def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2): |
---|
2593 | """Writes x,y,z,z,z coordinates of triangles constituting the bed |
---|
2594 | elevation. |
---|
2595 | FIXME: Not in use pt |
---|
2596 | """ |
---|
2597 | |
---|
2598 | M, N = v0.shape |
---|
2599 | |
---|
2600 | |
---|
2601 | print X0 |
---|
2602 | import sys; sys.exit() |
---|
2603 | FN = create_filename(filename, 'cpt', M) |
---|
2604 | print 'Writing to %s' %FN |
---|
2605 | |
---|
2606 | fid = open(FN, 'w') |
---|
2607 | for i in range(M): |
---|
2608 | for j in range(2): |
---|
2609 | fid.write('%.16e ' %X0[i,j]) #x, y |
---|
2610 | for j in range(N): |
---|
2611 | fid.write('%.16e ' %v0[i,j]) #z,z,z, |
---|
2612 | |
---|
2613 | for j in range(2): |
---|
2614 | fid.write('%.16e ' %X1[i,j]) #x, y |
---|
2615 | for j in range(N): |
---|
2616 | fid.write('%.16e ' %v1[i,j]) |
---|
2617 | |
---|
2618 | for j in range(2): |
---|
2619 | fid.write('%.16e ' %X2[i,j]) #x, y |
---|
2620 | for j in range(N): |
---|
2621 | fid.write('%.16e ' %v2[i,j]) |
---|
2622 | |
---|
2623 | fid.write('\n') |
---|
2624 | fid.close() |
---|
2625 | |
---|
2626 | |
---|
2627 | |
---|
2628 | #Function for storing out to e.g. visualisation |
---|
2629 | #FIXME: Do we want this? |
---|
2630 | #FIXME: Not done yet for this version |
---|
2631 | def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2): |
---|
2632 | """Writes x,y,z coordinates of triangles constituting the bed elevation. |
---|
2633 | """ |
---|
2634 | |
---|
2635 | M, N = v0.shape |
---|
2636 | |
---|
2637 | FN = create_filename(filename, 'dat', M) |
---|
2638 | #print 'Writing to %s' %FN |
---|
2639 | |
---|
2640 | fid = open(FN, 'w') |
---|
2641 | for i in range(M): |
---|
2642 | for j in range(2): |
---|
2643 | fid.write('%f ' %X0[i,j]) #x, y |
---|
2644 | fid.write('%f ' %v0[i,0]) #z |
---|
2645 | |
---|
2646 | for j in range(2): |
---|
2647 | fid.write('%f ' %X1[i,j]) #x, y |
---|
2648 | fid.write('%f ' %v1[i,0]) #z |
---|
2649 | |
---|
2650 | for j in range(2): |
---|
2651 | fid.write('%f ' %X2[i,j]) #x, y |
---|
2652 | fid.write('%f ' %v2[i,0]) #z |
---|
2653 | |
---|
2654 | fid.write('\n') |
---|
2655 | fid.close() |
---|
2656 | |
---|
2657 | |
---|
2658 | |
---|
2659 | def dat_variable_writer(filename, t, v0, v1, v2): |
---|
2660 | """Store water height to file |
---|
2661 | """ |
---|
2662 | |
---|
2663 | M, N = v0.shape |
---|
2664 | |
---|
2665 | FN = create_filename(filename, 'dat', M, t) |
---|
2666 | #print 'Writing to %s' %FN |
---|
2667 | |
---|
2668 | fid = open(FN, 'w') |
---|
2669 | for i in range(M): |
---|
2670 | fid.write('%.4f ' %v0[i,0]) |
---|
2671 | fid.write('%.4f ' %v1[i,0]) |
---|
2672 | fid.write('%.4f ' %v2[i,0]) |
---|
2673 | |
---|
2674 | fid.write('\n') |
---|
2675 | fid.close() |
---|
2676 | |
---|
2677 | |
---|
2678 | def read_sww(filename): |
---|
2679 | """Read sww Net CDF file containing Shallow Water Wave simulation |
---|
2680 | |
---|
2681 | The integer array volumes is of shape Nx3 where N is the number of |
---|
2682 | triangles in the mesh. |
---|
2683 | |
---|
2684 | Each entry in volumes is an index into the x,y arrays (the location). |
---|
2685 | |
---|
2686 | Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions |
---|
2687 | number_of_timesteps, number_of_points. |
---|
2688 | |
---|
2689 | The momentum is not always stored. |
---|
2690 | |
---|
2691 | """ |
---|
2692 | from Scientific.IO.NetCDF import NetCDFFile |
---|
2693 | print 'Reading from ', filename |
---|
2694 | fid = NetCDFFile(filename, 'r') #Open existing file for read |
---|
2695 | #latitude, longitude |
---|
2696 | # Get the variables as Numeric arrays |
---|
2697 | x = fid.variables['x'] #x-coordinates of vertices |
---|
2698 | y = fid.variables['y'] #y-coordinates of vertices |
---|
2699 | z = fid.variables['elevation'] #Elevation |
---|
2700 | time = fid.variables['time'] #Timesteps |
---|
2701 | stage = fid.variables['stage'] #Water level |
---|
2702 | #xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction |
---|
2703 | #ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction |
---|
2704 | |
---|
2705 | volumes = fid.variables['volumes'] #Connectivity |
---|
2706 | |
---|
2707 | |
---|
2708 | def decimate_dem(basename_in, stencil, cellsize_new, basename_out=None, |
---|
2709 | verbose=False): |
---|
2710 | """Read Digitial Elevation model from the following NetCDF format (.dem) |
---|
2711 | |
---|
2712 | Example: |
---|
2713 | |
---|
2714 | ncols 3121 |
---|
2715 | nrows 1800 |
---|
2716 | xllcorner 722000 |
---|
2717 | yllcorner 5893000 |
---|
2718 | cellsize 25 |
---|
2719 | NODATA_value -9999 |
---|
2720 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
2721 | |
---|
2722 | Decimate data to cellsize_new using stencil and write to NetCDF dem format. |
---|
2723 | """ |
---|
2724 | |
---|
2725 | import os |
---|
2726 | from Scientific.IO.NetCDF import NetCDFFile |
---|
2727 | from Numeric import Float, zeros, sum, reshape, equal |
---|
2728 | |
---|
2729 | root = basename_in |
---|
2730 | inname = root + '.dem' |
---|
2731 | |
---|
2732 | #Open existing netcdf file to read |
---|
2733 | infile = NetCDFFile(inname, 'r') |
---|
2734 | if verbose: print 'Reading DEM from %s' %inname |
---|
2735 | |
---|
2736 | #Read metadata |
---|
2737 | ncols = infile.ncols[0] |
---|
2738 | nrows = infile.nrows[0] |
---|
2739 | xllcorner = infile.xllcorner[0] |
---|
2740 | yllcorner = infile.yllcorner[0] |
---|
2741 | cellsize = infile.cellsize[0] |
---|
2742 | NODATA_value = infile.NODATA_value[0] |
---|
2743 | zone = infile.zone[0] |
---|
2744 | false_easting = infile.false_easting[0] |
---|
2745 | false_northing = infile.false_northing[0] |
---|
2746 | projection = infile.projection |
---|
2747 | datum = infile.datum |
---|
2748 | units = infile.units |
---|
2749 | |
---|
2750 | dem_elevation = infile.variables['elevation'] |
---|
2751 | |
---|
2752 | #Get output file name |
---|
2753 | if basename_out == None: |
---|
2754 | outname = root + '_' + repr(cellsize_new) + '.dem' |
---|
2755 | else: |
---|
2756 | outname = basename_out + '.dem' |
---|
2757 | |
---|
2758 | if verbose: print 'Write decimated NetCDF file to %s' %outname |
---|
2759 | |
---|
2760 | #Determine some dimensions for decimated grid |
---|
2761 | (nrows_stencil, ncols_stencil) = stencil.shape |
---|
2762 | x_offset = ncols_stencil / 2 |
---|
2763 | y_offset = nrows_stencil / 2 |
---|
2764 | cellsize_ratio = int(cellsize_new / cellsize) |
---|
2765 | ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio |
---|
2766 | nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio |
---|
2767 | |
---|
2768 | #Open netcdf file for output |
---|
2769 | outfile = NetCDFFile(outname, 'w') |
---|
2770 | |
---|
2771 | #Create new file |
---|
2772 | outfile.institution = 'Geoscience Australia' |
---|
2773 | outfile.description = 'NetCDF DEM format for compact and portable storage ' +\ |
---|
2774 | 'of spatial point data' |
---|
2775 | #Georeferencing |
---|
2776 | outfile.zone = zone |
---|
2777 | outfile.projection = projection |
---|
2778 | outfile.datum = datum |
---|
2779 | outfile.units = units |
---|
2780 | |
---|
2781 | outfile.cellsize = cellsize_new |
---|
2782 | outfile.NODATA_value = NODATA_value |
---|
2783 | outfile.false_easting = false_easting |
---|
2784 | outfile.false_northing = false_northing |
---|
2785 | |
---|
2786 | outfile.xllcorner = xllcorner + (x_offset * cellsize) |
---|
2787 | outfile.yllcorner = yllcorner + (y_offset * cellsize) |
---|
2788 | outfile.ncols = ncols_new |
---|
2789 | outfile.nrows = nrows_new |
---|
2790 | |
---|
2791 | # dimension definition |
---|
2792 | outfile.createDimension('number_of_points', nrows_new*ncols_new) |
---|
2793 | |
---|
2794 | # variable definition |
---|
2795 | outfile.createVariable('elevation', Float, ('number_of_points',)) |
---|
2796 | |
---|
2797 | # Get handle to the variable |
---|
2798 | elevation = outfile.variables['elevation'] |
---|
2799 | |
---|
2800 | dem_elevation_r = reshape(dem_elevation, (nrows, ncols)) |
---|
2801 | |
---|
2802 | #Store data |
---|
2803 | global_index = 0 |
---|
2804 | for i in range(nrows_new): |
---|
2805 | if verbose: print 'Processing row %d of %d' %(i, nrows_new) |
---|
2806 | lower_index = global_index |
---|
2807 | telev = zeros(ncols_new, Float) |
---|
2808 | local_index = 0 |
---|
2809 | trow = i * cellsize_ratio |
---|
2810 | |
---|
2811 | for j in range(ncols_new): |
---|
2812 | tcol = j * cellsize_ratio |
---|
2813 | tmp = dem_elevation_r[trow:trow+nrows_stencil, tcol:tcol+ncols_stencil] |
---|
2814 | |
---|
2815 | #if dem contains 1 or more NODATA_values set value in |
---|
2816 | #decimated dem to NODATA_value, else compute decimated |
---|
2817 | #value using stencil |
---|
2818 | if sum(sum(equal(tmp, NODATA_value))) > 0: |
---|
2819 | telev[local_index] = NODATA_value |
---|
2820 | else: |
---|
2821 | telev[local_index] = sum(sum(tmp * stencil)) |
---|
2822 | |
---|
2823 | global_index += 1 |
---|
2824 | local_index += 1 |
---|
2825 | |
---|
2826 | upper_index = global_index |
---|
2827 | |
---|
2828 | elevation[lower_index:upper_index] = telev |
---|
2829 | |
---|
2830 | assert global_index == nrows_new*ncols_new, 'index not equal to number of points' |
---|
2831 | |
---|
2832 | infile.close() |
---|
2833 | outfile.close() |
---|
2834 | |
---|
2835 | |
---|
2836 | |
---|
2837 | def sww2asc_obsolete(basename_in, basename_out = None, |
---|
2838 | quantity = None, |
---|
2839 | timestep = None, |
---|
2840 | reduction = None, |
---|
2841 | cellsize = 10, |
---|
2842 | verbose = False, |
---|
2843 | origin = None): |
---|
2844 | """Read SWW file and convert to Digitial Elevation model format (.asc) |
---|
2845 | |
---|
2846 | Example: |
---|
2847 | |
---|
2848 | ncols 3121 |
---|
2849 | nrows 1800 |
---|
2850 | xllcorner 722000 |
---|
2851 | yllcorner 5893000 |
---|
2852 | cellsize 25 |
---|
2853 | NODATA_value -9999 |
---|
2854 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
2855 | |
---|
2856 | Also write accompanying file with same basename_in but extension .prj |
---|
2857 | used to fix the UTM zone, datum, false northings and eastings. |
---|
2858 | |
---|
2859 | The prj format is assumed to be as |
---|
2860 | |
---|
2861 | Projection UTM |
---|
2862 | Zone 56 |
---|
2863 | Datum WGS84 |
---|
2864 | Zunits NO |
---|
2865 | Units METERS |
---|
2866 | Spheroid WGS84 |
---|
2867 | Xshift 0.0000000000 |
---|
2868 | Yshift 10000000.0000000000 |
---|
2869 | Parameters |
---|
2870 | |
---|
2871 | |
---|
2872 | if quantity is given, out values from quantity otherwise default to |
---|
2873 | elevation |
---|
2874 | |
---|
2875 | if timestep (an index) is given, output quantity at that timestep |
---|
2876 | |
---|
2877 | if reduction is given use that to reduce quantity over all timesteps. |
---|
2878 | |
---|
2879 | """ |
---|
2880 | from Numeric import array, Float, concatenate, NewAxis, zeros,\ |
---|
2881 | sometrue |
---|
2882 | |
---|
2883 | |
---|
2884 | #FIXME: Should be variable |
---|
2885 | datum = 'WGS84' |
---|
2886 | false_easting = 500000 |
---|
2887 | false_northing = 10000000 |
---|
2888 | |
---|
2889 | if quantity is None: |
---|
2890 | quantity = 'elevation' |
---|
2891 | |
---|
2892 | if reduction is None: |
---|
2893 | reduction = max |
---|
2894 | |
---|
2895 | if basename_out is None: |
---|
2896 | basename_out = basename_in + '_%s' %quantity |
---|
2897 | |
---|
2898 | swwfile = basename_in + '.sww' |
---|
2899 | ascfile = basename_out + '.asc' |
---|
2900 | prjfile = basename_out + '.prj' |
---|
2901 | |
---|
2902 | |
---|
2903 | if verbose: print 'Reading from %s' %swwfile |
---|
2904 | #Read sww file |
---|
2905 | from Scientific.IO.NetCDF import NetCDFFile |
---|
2906 | fid = NetCDFFile(swwfile) |
---|
2907 | |
---|
2908 | #Get extent and reference |
---|
2909 | x = fid.variables['x'][:] |
---|
2910 | y = fid.variables['y'][:] |
---|
2911 | volumes = fid.variables['volumes'][:] |
---|
2912 | |
---|
2913 | ymin = min(y); ymax = max(y) |
---|
2914 | xmin = min(x); xmax = max(x) |
---|
2915 | |
---|
2916 | number_of_timesteps = fid.dimensions['number_of_timesteps'] |
---|
2917 | number_of_points = fid.dimensions['number_of_points'] |
---|
2918 | if origin is None: |
---|
2919 | |
---|
2920 | #Get geo_reference |
---|
2921 | #sww files don't have to have a geo_ref |
---|
2922 | try: |
---|
2923 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
2924 | except AttributeError, e: |
---|
2925 | geo_reference = Geo_reference() #Default georef object |
---|
2926 | |
---|
2927 | xllcorner = geo_reference.get_xllcorner() |
---|
2928 | yllcorner = geo_reference.get_yllcorner() |
---|
2929 | zone = geo_reference.get_zone() |
---|
2930 | else: |
---|
2931 | zone = origin[0] |
---|
2932 | xllcorner = origin[1] |
---|
2933 | yllcorner = origin[2] |
---|
2934 | |
---|
2935 | |
---|
2936 | #Get quantity and reduce if applicable |
---|
2937 | if verbose: print 'Reading quantity %s' %quantity |
---|
2938 | |
---|
2939 | if quantity.lower() == 'depth': |
---|
2940 | q = fid.variables['stage'][:] - fid.variables['elevation'][:] |
---|
2941 | else: |
---|
2942 | q = fid.variables[quantity][:] |
---|
2943 | |
---|
2944 | |
---|
2945 | if len(q.shape) == 2: |
---|
2946 | if verbose: print 'Reducing quantity %s' %quantity |
---|
2947 | q_reduced = zeros( number_of_points, Float ) |
---|
2948 | |
---|
2949 | for k in range(number_of_points): |
---|
2950 | q_reduced[k] = reduction( q[:,k] ) |
---|
2951 | |
---|
2952 | q = q_reduced |
---|
2953 | |
---|
2954 | #Now q has dimension: number_of_points |
---|
2955 | |
---|
2956 | #Create grid and update xll/yll corner and x,y |
---|
2957 | if verbose: print 'Creating grid' |
---|
2958 | ncols = int((xmax-xmin)/cellsize)+1 |
---|
2959 | nrows = int((ymax-ymin)/cellsize)+1 |
---|
2960 | |
---|
2961 | newxllcorner = xmin+xllcorner |
---|
2962 | newyllcorner = ymin+yllcorner |
---|
2963 | |
---|
2964 | x = x+xllcorner-newxllcorner |
---|
2965 | y = y+yllcorner-newyllcorner |
---|
2966 | |
---|
2967 | vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) |
---|
2968 | assert len(vertex_points.shape) == 2 |
---|
2969 | |
---|
2970 | |
---|
2971 | from Numeric import zeros, Float |
---|
2972 | grid_points = zeros ( (ncols*nrows, 2), Float ) |
---|
2973 | |
---|
2974 | |
---|
2975 | for i in xrange(nrows): |
---|
2976 | yg = i*cellsize |
---|
2977 | for j in xrange(ncols): |
---|
2978 | xg = j*cellsize |
---|
2979 | k = i*ncols + j |
---|
2980 | |
---|
2981 | grid_points[k,0] = xg |
---|
2982 | grid_points[k,1] = yg |
---|
2983 | |
---|
2984 | #Interpolate |
---|
2985 | from least_squares import Interpolation |
---|
2986 | from util import inside_polygon |
---|
2987 | |
---|
2988 | #FIXME: This should be done with precrop = True, otherwise it'll |
---|
2989 | #take forever. With expand_search set to False, some grid points might |
---|
2990 | #miss out.... |
---|
2991 | |
---|
2992 | interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0, |
---|
2993 | precrop = False, expand_search = False, |
---|
2994 | verbose = verbose) |
---|
2995 | |
---|
2996 | #Interpolate using quantity values |
---|
2997 | if verbose: print 'Interpolating' |
---|
2998 | grid_values = interp.interpolate(q).flat |
---|
2999 | |
---|
3000 | #Write |
---|
3001 | #Write prj file |
---|
3002 | if verbose: print 'Writing %s' %prjfile |
---|
3003 | prjid = open(prjfile, 'w') |
---|
3004 | prjid.write('Projection %s\n' %'UTM') |
---|
3005 | prjid.write('Zone %d\n' %zone) |
---|
3006 | prjid.write('Datum %s\n' %datum) |
---|
3007 | prjid.write('Zunits NO\n') |
---|
3008 | prjid.write('Units METERS\n') |
---|
3009 | prjid.write('Spheroid %s\n' %datum) |
---|
3010 | prjid.write('Xshift %d\n' %false_easting) |
---|
3011 | prjid.write('Yshift %d\n' %false_northing) |
---|
3012 | prjid.write('Parameters\n') |
---|
3013 | prjid.close() |
---|
3014 | |
---|
3015 | |
---|
3016 | |
---|
3017 | if verbose: print 'Writing %s' %ascfile |
---|
3018 | NODATA_value = -9999 |
---|
3019 | |
---|
3020 | ascid = open(ascfile, 'w') |
---|
3021 | |
---|
3022 | ascid.write('ncols %d\n' %ncols) |
---|
3023 | ascid.write('nrows %d\n' %nrows) |
---|
3024 | ascid.write('xllcorner %d\n' %newxllcorner) |
---|
3025 | ascid.write('yllcorner %d\n' %newyllcorner) |
---|
3026 | ascid.write('cellsize %f\n' %cellsize) |
---|
3027 | ascid.write('NODATA_value %d\n' %NODATA_value) |
---|
3028 | |
---|
3029 | |
---|
3030 | #Get bounding polygon from mesh |
---|
3031 | P = interp.mesh.get_boundary_polygon() |
---|
3032 | inside_indices = inside_polygon(grid_points, P) |
---|
3033 | |
---|
3034 | for i in range(nrows): |
---|
3035 | if verbose and i%((nrows+10)/10)==0: |
---|
3036 | print 'Doing row %d of %d' %(i, nrows) |
---|
3037 | |
---|
3038 | for j in range(ncols): |
---|
3039 | index = (nrows-i-1)*ncols+j |
---|
3040 | |
---|
3041 | if sometrue(inside_indices == index): |
---|
3042 | ascid.write('%f ' %grid_values[index]) |
---|
3043 | else: |
---|
3044 | ascid.write('%d ' %NODATA_value) |
---|
3045 | |
---|
3046 | ascid.write('\n') |
---|
3047 | |
---|
3048 | #Close |
---|
3049 | ascid.close() |
---|
3050 | fid.close() |
---|
3051 | |
---|
3052 | def sww2ers_obsolete(basename_in, basename_out = None, |
---|
3053 | quantity = None, |
---|
3054 | timestep = None, |
---|
3055 | reduction = None, |
---|
3056 | cellsize = 10, |
---|
3057 | verbose = False, |
---|
3058 | origin = None, |
---|
3059 | datum = 'WGS84'): |
---|
3060 | |
---|
3061 | """Read SWW file and convert to Digitial Elevation model format (.asc) |
---|
3062 | |
---|
3063 | Example: |
---|
3064 | |
---|
3065 | ncols 3121 |
---|
3066 | nrows 1800 |
---|
3067 | xllcorner 722000 |
---|
3068 | yllcorner 5893000 |
---|
3069 | cellsize 25 |
---|
3070 | NODATA_value -9999 |
---|
3071 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
3072 | |
---|
3073 | Also write accompanying file with same basename_in but extension .prj |
---|
3074 | used to fix the UTM zone, datum, false northings and eastings. |
---|
3075 | |
---|
3076 | The prj format is assumed to be as |
---|
3077 | |
---|
3078 | Projection UTM |
---|
3079 | Zone 56 |
---|
3080 | Datum WGS84 |
---|
3081 | Zunits NO |
---|
3082 | Units METERS |
---|
3083 | Spheroid WGS84 |
---|
3084 | Xshift 0.0000000000 |
---|
3085 | Yshift 10000000.0000000000 |
---|
3086 | Parameters |
---|
3087 | |
---|
3088 | |
---|
3089 | if quantity is given, out values from quantity otherwise default to |
---|
3090 | elevation |
---|
3091 | |
---|
3092 | if timestep (an index) is given, output quantity at that timestep |
---|
3093 | |
---|
3094 | if reduction is given use that to reduce quantity over all timesteps. |
---|
3095 | |
---|
3096 | """ |
---|
3097 | from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue |
---|
3098 | import ermapper_grids |
---|
3099 | |
---|
3100 | header = {} |
---|
3101 | false_easting = 500000 |
---|
3102 | false_northing = 10000000 |
---|
3103 | NODATA_value = 0 |
---|
3104 | |
---|
3105 | if quantity is None: |
---|
3106 | quantity = 'elevation' |
---|
3107 | |
---|
3108 | if reduction is None: |
---|
3109 | reduction = max |
---|
3110 | |
---|
3111 | if basename_out is None: |
---|
3112 | basename_out = basename_in + '_%s' %quantity |
---|
3113 | |
---|
3114 | swwfile = basename_in + '.sww' |
---|
3115 | # Note the use of a .ers extension is optional (write_ermapper_grid will |
---|
3116 | # deal with either option |
---|
3117 | ersfile = basename_out |
---|
3118 | ## prjfile = basename_out + '.prj' |
---|
3119 | |
---|
3120 | |
---|
3121 | if verbose: print 'Reading from %s' %swwfile |
---|
3122 | #Read sww file |
---|
3123 | from Scientific.IO.NetCDF import NetCDFFile |
---|
3124 | fid = NetCDFFile(swwfile) |
---|
3125 | |
---|
3126 | #Get extent and reference |
---|
3127 | x = fid.variables['x'][:] |
---|
3128 | y = fid.variables['y'][:] |
---|
3129 | volumes = fid.variables['volumes'][:] |
---|
3130 | |
---|
3131 | ymin = min(y); ymax = max(y) |
---|
3132 | xmin = min(x); xmax = max(x) |
---|
3133 | |
---|
3134 | number_of_timesteps = fid.dimensions['number_of_timesteps'] |
---|
3135 | number_of_points = fid.dimensions['number_of_points'] |
---|
3136 | if origin is None: |
---|
3137 | |
---|
3138 | #Get geo_reference |
---|
3139 | #sww files don't have to have a geo_ref |
---|
3140 | try: |
---|
3141 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
3142 | except AttributeError, e: |
---|
3143 | geo_reference = Geo_reference() #Default georef object |
---|
3144 | |
---|
3145 | xllcorner = geo_reference.get_xllcorner() |
---|
3146 | yllcorner = geo_reference.get_yllcorner() |
---|
3147 | zone = geo_reference.get_zone() |
---|
3148 | else: |
---|
3149 | zone = origin[0] |
---|
3150 | xllcorner = origin[1] |
---|
3151 | yllcorner = origin[2] |
---|
3152 | |
---|
3153 | |
---|
3154 | #Get quantity and reduce if applicable |
---|
3155 | if verbose: print 'Reading quantity %s' %quantity |
---|
3156 | |
---|
3157 | if quantity.lower() == 'depth': |
---|
3158 | q = fid.variables['stage'][:] - fid.variables['elevation'][:] |
---|
3159 | else: |
---|
3160 | q = fid.variables[quantity][:] |
---|
3161 | |
---|
3162 | |
---|
3163 | if len(q.shape) == 2: |
---|
3164 | if verbose: print 'Reducing quantity %s' %quantity |
---|
3165 | q_reduced = zeros( number_of_points, Float ) |
---|
3166 | |
---|
3167 | for k in range(number_of_points): |
---|
3168 | q_reduced[k] = reduction( q[:,k] ) |
---|
3169 | |
---|
3170 | q = q_reduced |
---|
3171 | |
---|
3172 | #Now q has dimension: number_of_points |
---|
3173 | |
---|
3174 | #Create grid and update xll/yll corner and x,y |
---|
3175 | if verbose: print 'Creating grid' |
---|
3176 | ncols = int((xmax-xmin)/cellsize)+1 |
---|
3177 | nrows = int((ymax-ymin)/cellsize)+1 |
---|
3178 | |
---|
3179 | newxllcorner = xmin+xllcorner |
---|
3180 | newyllcorner = ymin+yllcorner |
---|
3181 | |
---|
3182 | x = x+xllcorner-newxllcorner |
---|
3183 | y = y+yllcorner-newyllcorner |
---|
3184 | |
---|
3185 | vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) |
---|
3186 | assert len(vertex_points.shape) == 2 |
---|
3187 | |
---|
3188 | |
---|
3189 | from Numeric import zeros, Float |
---|
3190 | grid_points = zeros ( (ncols*nrows, 2), Float ) |
---|
3191 | |
---|
3192 | |
---|
3193 | for i in xrange(nrows): |
---|
3194 | yg = (nrows-i)*cellsize # this will flip the order of the y values |
---|
3195 | for j in xrange(ncols): |
---|
3196 | xg = j*cellsize |
---|
3197 | k = i*ncols + j |
---|
3198 | |
---|
3199 | grid_points[k,0] = xg |
---|
3200 | grid_points[k,1] = yg |
---|
3201 | |
---|
3202 | #Interpolate |
---|
3203 | from least_squares import Interpolation |
---|
3204 | from util import inside_polygon |
---|
3205 | |
---|
3206 | #FIXME: This should be done with precrop = True (?), otherwise it'll |
---|
3207 | #take forever. With expand_search set to False, some grid points might |
---|
3208 | #miss out.... |
---|
3209 | |
---|
3210 | interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0, |
---|
3211 | precrop = False, expand_search = False, |
---|
3212 | verbose = verbose) |
---|
3213 | |
---|
3214 | #Interpolate using quantity values |
---|
3215 | if verbose: print 'Interpolating' |
---|
3216 | grid_values = interp.interpolate(q).flat |
---|
3217 | grid_values = reshape(grid_values,(nrows, ncols)) |
---|
3218 | |
---|
3219 | |
---|
3220 | # setup header information |
---|
3221 | header['datum'] = '"' + datum + '"' |
---|
3222 | # FIXME The use of hardwired UTM and zone number needs to be made optional |
---|
3223 | # FIXME Also need an automatic test for coordinate type (i.e. EN or LL) |
---|
3224 | header['projection'] = '"UTM-' + str(zone) + '"' |
---|
3225 | header['coordinatetype'] = 'EN' |
---|
3226 | if header['coordinatetype'] == 'LL': |
---|
3227 | header['longitude'] = str(newxllcorner) |
---|
3228 | header['latitude'] = str(newyllcorner) |
---|
3229 | elif header['coordinatetype'] == 'EN': |
---|
3230 | header['eastings'] = str(newxllcorner) |
---|
3231 | header['northings'] = str(newyllcorner) |
---|
3232 | header['nullcellvalue'] = str(NODATA_value) |
---|
3233 | header['xdimension'] = str(cellsize) |
---|
3234 | header['ydimension'] = str(cellsize) |
---|
3235 | header['value'] = '"' + quantity + '"' |
---|
3236 | |
---|
3237 | |
---|
3238 | #Write |
---|
3239 | if verbose: print 'Writing %s' %ersfile |
---|
3240 | ermapper_grids.write_ermapper_grid(ersfile,grid_values, header) |
---|
3241 | |
---|
3242 | fid.close() |
---|
3243 | |
---|
3244 | ## ascid = open(ascfile, 'w') |
---|
3245 | ## |
---|
3246 | ## ascid.write('ncols %d\n' %ncols) |
---|
3247 | ## ascid.write('nrows %d\n' %nrows) |
---|
3248 | ## ascid.write('xllcorner %d\n' %newxllcorner) |
---|
3249 | ## ascid.write('yllcorner %d\n' %newyllcorner) |
---|
3250 | ## ascid.write('cellsize %f\n' %cellsize) |
---|
3251 | ## ascid.write('NODATA_value %d\n' %NODATA_value) |
---|
3252 | ## |
---|
3253 | ## |
---|
3254 | ## #Get bounding polygon from mesh |
---|
3255 | ## P = interp.mesh.get_boundary_polygon() |
---|
3256 | ## inside_indices = inside_polygon(grid_points, P) |
---|
3257 | ## |
---|
3258 | ## for i in range(nrows): |
---|
3259 | ## if verbose and i%((nrows+10)/10)==0: |
---|
3260 | ## print 'Doing row %d of %d' %(i, nrows) |
---|
3261 | ## |
---|
3262 | ## for j in range(ncols): |
---|
3263 | ## index = (nrows-i-1)*ncols+j |
---|
3264 | ## |
---|
3265 | ## if sometrue(inside_indices == index): |
---|
3266 | ## ascid.write('%f ' %grid_values[index]) |
---|
3267 | ## else: |
---|
3268 | ## ascid.write('%d ' %NODATA_value) |
---|
3269 | ## |
---|
3270 | ## ascid.write('\n') |
---|
3271 | ## |
---|
3272 | ## #Close |
---|
3273 | ## ascid.close() |
---|
3274 | ## fid.close() |
---|
3275 | |
---|