1 | import os |
---|
2 | from struct import unpack |
---|
3 | from anuga.utilities.anuga_exceptions import ANUGAError |
---|
4 | from Numeric import Float, Int, zeros |
---|
5 | from Numeric import concatenate, array, Float, Int, Int32, resize, \ |
---|
6 | sometrue, searchsorted, zeros, allclose, around, reshape, \ |
---|
7 | transpose, sort, NewAxis, ArrayType, compress, take, arange, \ |
---|
8 | argmax, alltrue, shape, Float32, size |
---|
9 | |
---|
10 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
11 | from anuga.config import max_float |
---|
12 | from anuga.coordinate_transforms.geo_reference import Geo_reference, \ |
---|
13 | write_NetCDF_georeference, ensure_geo_reference |
---|
14 | """ |
---|
15 | The name of the urs file names must be; |
---|
16 | [basename_in]_velocity-z-mux2 |
---|
17 | [basename_in]_velocity-e-mux2 |
---|
18 | [basename_in]_waveheight-n-mux2 |
---|
19 | """ |
---|
20 | def read_binary_mux2(file_in,verbose=False): |
---|
21 | """ |
---|
22 | Program for demuxing small format (mux2) files. Based upon |
---|
23 | demux_compact_tgdata.c written by C. Thomas Geoscience Austrlia 2007 |
---|
24 | Author: John Jakeman |
---|
25 | |
---|
26 | Unlike c code this program only reads in one file (source) |
---|
27 | """ |
---|
28 | eps=0.00001 |
---|
29 | |
---|
30 | ###########E###### |
---|
31 | # Open mux2 File # |
---|
32 | ################## |
---|
33 | |
---|
34 | if (os.access(file_in, os.F_OK) == 0): |
---|
35 | msg = 'File %s does not exist or is not accessible' %file_name |
---|
36 | raise IOError, msg |
---|
37 | if file_in[-4:]!="mux2": |
---|
38 | msg ="This program operates only on multiplexed files in mux2 format\nAt least one input file name does not end with mux2" |
---|
39 | raise IOError, msg |
---|
40 | |
---|
41 | mux_file = open(file_in, 'rb') |
---|
42 | |
---|
43 | ###################### |
---|
44 | # Read in the header # |
---|
45 | ###################### |
---|
46 | |
---|
47 | # Number of points/stations |
---|
48 | (num_sites,)= unpack('i',mux_file.read(4)) |
---|
49 | |
---|
50 | geolat=zeros(num_sites,Float) |
---|
51 | geolon=zeros(num_sites,Float) |
---|
52 | depth=zeros(num_sites,Float) |
---|
53 | dt_old=0.0 |
---|
54 | nt_old=0.0 |
---|
55 | for i in range(num_sites): |
---|
56 | # location and water depth in geographic coordinates |
---|
57 | #-180/180 |
---|
58 | #-90/90 |
---|
59 | (geolat[i],)= unpack('f',mux_file.read(4)) |
---|
60 | |
---|
61 | #dt, float - time step, seconds |
---|
62 | (geolon[i],) = unpack('f', mux_file.read(4)) |
---|
63 | |
---|
64 | (mcolat,) = unpack('f', mux_file.read(4)) |
---|
65 | (mcolon,) = unpack('f', mux_file.read(4)) |
---|
66 | (ig,) = unpack('i', mux_file.read(4)) |
---|
67 | (ilon,) = unpack('i', mux_file.read(4)) #grid point location |
---|
68 | (ilat,) = unpack('i', mux_file.read(4)) |
---|
69 | (depth[i],) = unpack('f', mux_file.read(4)) |
---|
70 | (centerlat,) = unpack('f', mux_file.read(4)) |
---|
71 | (centerlon,) = unpack('f', mux_file.read(4)) |
---|
72 | (offset,) = unpack('f', mux_file.read(4)) |
---|
73 | (az,) = unpack('f', mux_file.read(4)) |
---|
74 | (baz,) = unpack('f', mux_file.read(4)) |
---|
75 | (dt,) = unpack('f', mux_file.read(4)) |
---|
76 | (nt,) = unpack('i', mux_file.read(4)) |
---|
77 | for j in range(4): # identifier |
---|
78 | (id,) = unpack('f', mux_file.read(4)) |
---|
79 | |
---|
80 | msg = "Bad data in the mux file." |
---|
81 | if num_sites < 0: |
---|
82 | mux_file.close() |
---|
83 | raise ANUGAError, msg |
---|
84 | if dt < 0: |
---|
85 | mux_file.close() |
---|
86 | raise ANUGAError, msg |
---|
87 | if nt < 0: |
---|
88 | mux_file.close() |
---|
89 | raise ANUGAError, msg |
---|
90 | |
---|
91 | if (verbose): |
---|
92 | print "num_stations:",num_sites |
---|
93 | print "station number:",i |
---|
94 | print "geolat:",geolat[i] |
---|
95 | print "geolon:",geolon[i] |
---|
96 | print "mcolat:",mcolat |
---|
97 | print "mcolo:",mcolon |
---|
98 | print "ig:",ig |
---|
99 | print "ilon:",ilon |
---|
100 | print "ilat:",ilat |
---|
101 | print "depth:",depth[i] |
---|
102 | print "centerlat:",centerlat |
---|
103 | print "centerlon:",centerlon |
---|
104 | print "offset:",offset |
---|
105 | print "az:",az |
---|
106 | print "baz:",baz |
---|
107 | print "tstep",dt |
---|
108 | print "num_tsteps",nt |
---|
109 | print |
---|
110 | |
---|
111 | if i>0: |
---|
112 | msg="Stations have different time step size" |
---|
113 | assert (dt==dt_old),msg |
---|
114 | msg="Stations have different time series length" |
---|
115 | assert (nt==nt_old),msg |
---|
116 | dt_old=dt |
---|
117 | nt_old=nt |
---|
118 | |
---|
119 | ################ |
---|
120 | # Read in Data # |
---|
121 | ################ |
---|
122 | |
---|
123 | # Make an array to hold the start and stop steps for each station for source |
---|
124 | first_tstep=zeros(num_sites,Int) |
---|
125 | last_tstep=zeros(num_sites,Int) |
---|
126 | # Read the start and stop steps for each site into the array for source 0 |
---|
127 | for i in range(num_sites): |
---|
128 | (first_tstep[i],)=unpack('i', mux_file.read(4)) |
---|
129 | for i in range(num_sites): |
---|
130 | (last_tstep[i],)=unpack('i', mux_file.read(4)) |
---|
131 | |
---|
132 | # Compute the size of the data block |
---|
133 | numDataTotal=0 |
---|
134 | numData=zeros(num_sites,Int) |
---|
135 | last_output_step=0 |
---|
136 | for i in range(num_sites): |
---|
137 | numDataTotal+= last_tstep[i]-first_tstep[i]+1 |
---|
138 | numData[i]= last_tstep[i]-first_tstep[i]+1 |
---|
139 | last_output_step=max(last_output_step,last_tstep[i]) |
---|
140 | numDataTotal+=last_output_step |
---|
141 | |
---|
142 | msg="Size of Data Block is negative" |
---|
143 | assert numDataTotal>=0,msg |
---|
144 | |
---|
145 | muxData = zeros(numDataTotal,Float) |
---|
146 | for i in range(numDataTotal): |
---|
147 | (muxData[i],)=unpack('f', mux_file.read(4)) |
---|
148 | |
---|
149 | ################## |
---|
150 | # Store Mux data # |
---|
151 | ################## |
---|
152 | |
---|
153 | # Make arrays of starting and finishing time steps for the tide gauges |
---|
154 | # and fill them from the file |
---|
155 | data=zeros((num_sites,nt),Float) |
---|
156 | maximum=zeros(num_sites,Float) |
---|
157 | |
---|
158 | # Find when first station starts recording |
---|
159 | min_tstep = min(first_tstep) #not necessary |
---|
160 | # Find when all stations have stoped recording |
---|
161 | max_tstep = max(last_tstep) |
---|
162 | |
---|
163 | times=zeros(max_tstep,Float) |
---|
164 | beg=dt |
---|
165 | |
---|
166 | for i in range(num_sites): |
---|
167 | if (verbose): print 'reading site '+str(i) |
---|
168 | istart=-1 |
---|
169 | istop=-1 |
---|
170 | offset=0 |
---|
171 | # Update start and stop timesteps for this gauge |
---|
172 | if istart==-1: |
---|
173 | istart=first_tstep[i] |
---|
174 | else: |
---|
175 | istart=min(first_tstep[i],istart) |
---|
176 | if istop==-1: |
---|
177 | istop=last_tstep[i] |
---|
178 | else: |
---|
179 | istop=min(last_tstep[i],istop) |
---|
180 | for t in range(max_tstep): |
---|
181 | last_t=t |
---|
182 | offset+=1 |
---|
183 | # Skip records from earlier tide gauges |
---|
184 | for j in range(i): |
---|
185 | if (t+1>=first_tstep[j]) and (t+1<=last_tstep[j]): |
---|
186 | offset+=1 |
---|
187 | # Deal with the tide gauge at hand |
---|
188 | if (t+1>=first_tstep[i]) and (t+1<=last_tstep[i]): |
---|
189 | # Gauge is recording at this time] |
---|
190 | data[i][t]=muxData[offset] |
---|
191 | offset+=1 |
---|
192 | elif (t+1<first_tstep[i]): |
---|
193 | # Gauge has not yet started recording |
---|
194 | data[i][t]=0.0 |
---|
195 | else: |
---|
196 | # Gauge has finished recording |
---|
197 | data[i][t]=0.0 |
---|
198 | break |
---|
199 | # Skip records from later tide gauges |
---|
200 | for j in range(i+1,num_sites): |
---|
201 | if (t+1>=first_tstep[j]) and (t+1<=last_tstep[j]): |
---|
202 | offset+=1 |
---|
203 | if i==0: |
---|
204 | times[t]=beg+t*dt |
---|
205 | |
---|
206 | if (last_t<max_tstep-1): |
---|
207 | # The loop was exited early because the gauge had finished recording |
---|
208 | for t in range(last_t+1,max_tstep): |
---|
209 | data[i][t]=0.0 #pad until all stations have stopped recording |
---|
210 | |
---|
211 | ####################################### |
---|
212 | # Compute the maxima for each station # |
---|
213 | ####################################### |
---|
214 | |
---|
215 | for t in range(max_tstep): |
---|
216 | if (data[i][t] > eps): |
---|
217 | maximum[i]=max(data[i][t],maximum[i]) |
---|
218 | if i==0: |
---|
219 | times[t]=beg+t*dt |
---|
220 | |
---|
221 | msg = "Mux file corrupted. No positive height found at station "+str(i) |
---|
222 | assert maximum[i]>0, msg |
---|
223 | |
---|
224 | return times, geolat, geolon, -depth, data |
---|
225 | |
---|
226 | |
---|
227 | class Write_sts: |
---|
228 | |
---|
229 | sts_quantities = ['stage'] |
---|
230 | |
---|
231 | |
---|
232 | RANGE = '_range' |
---|
233 | EXTREMA = ':extrema' |
---|
234 | |
---|
235 | def __init__(self): |
---|
236 | pass |
---|
237 | |
---|
238 | def store_header(self, |
---|
239 | outfile, |
---|
240 | times, |
---|
241 | number_of_points, |
---|
242 | description='Converted from URS mux2 format', |
---|
243 | sts_precision=Float32, |
---|
244 | verbose=False): |
---|
245 | """ |
---|
246 | outfile - the name of the file that will be written |
---|
247 | times - A list of the time slice times OR a start time |
---|
248 | Note, if a list is given the info will be made relative. |
---|
249 | number_of_points - the number of urs gauges sites |
---|
250 | """ |
---|
251 | |
---|
252 | outfile.institution = 'Geoscience Australia' |
---|
253 | outfile.description = description |
---|
254 | |
---|
255 | try: |
---|
256 | revision_number = get_revision_number() |
---|
257 | except: |
---|
258 | revision_number = None |
---|
259 | # Allow None to be stored as a string |
---|
260 | outfile.revision_number = str(revision_number) |
---|
261 | |
---|
262 | # times - A list or array of the time slice times OR a start time |
---|
263 | # Start time in seconds since the epoch (midnight 1/1/1970) |
---|
264 | |
---|
265 | # This is being used to seperate one number from a list. |
---|
266 | # what it is actually doing is sorting lists from numeric arrays. |
---|
267 | if type(times) is list or type(times) is ArrayType: |
---|
268 | number_of_times = len(times) |
---|
269 | times = ensure_numeric(times) |
---|
270 | if number_of_times == 0: |
---|
271 | starttime = 0 |
---|
272 | else: |
---|
273 | starttime = times[0] |
---|
274 | times = times - starttime #Store relative times |
---|
275 | else: |
---|
276 | number_of_times = 0 |
---|
277 | starttime = times |
---|
278 | |
---|
279 | outfile.starttime = starttime |
---|
280 | |
---|
281 | # Dimension definitions |
---|
282 | outfile.createDimension('number_of_points', number_of_points) |
---|
283 | outfile.createDimension('number_of_timesteps', number_of_times) |
---|
284 | outfile.createDimension('numbers_in_range', 2) |
---|
285 | |
---|
286 | # Variable definitions |
---|
287 | outfile.createVariable('x', sts_precision, ('number_of_points',)) |
---|
288 | outfile.createVariable('y', sts_precision, ('number_of_points',)) |
---|
289 | outfile.createVariable('elevation', sts_precision,('number_of_points',)) |
---|
290 | |
---|
291 | q = 'elevation' |
---|
292 | outfile.createVariable(q+Write_sts.RANGE, sts_precision, |
---|
293 | ('numbers_in_range',)) |
---|
294 | |
---|
295 | # Initialise ranges with small and large sentinels. |
---|
296 | # If this was in pure Python we could have used None sensibly |
---|
297 | outfile.variables[q+Write_sts.RANGE][0] = max_float # Min |
---|
298 | outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max |
---|
299 | |
---|
300 | outfile.createVariable('z', sts_precision, ('number_of_points',)) |
---|
301 | # Doing sts_precision instead of Float gives cast errors. |
---|
302 | outfile.createVariable('time', Float, ('number_of_timesteps',)) |
---|
303 | |
---|
304 | for q in Write_sts.sts_quantities: |
---|
305 | outfile.createVariable(q, sts_precision, |
---|
306 | ('number_of_timesteps', |
---|
307 | 'number_of_points')) |
---|
308 | outfile.createVariable(q+Write_sts.RANGE, sts_precision, |
---|
309 | ('numbers_in_range',)) |
---|
310 | # Initialise ranges with small and large sentinels. |
---|
311 | # If this was in pure Python we could have used None sensibly |
---|
312 | outfile.variables[q+Write_sts.RANGE][0] = max_float # Min |
---|
313 | outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max |
---|
314 | |
---|
315 | if type(times) is list or type(times) is ArrayType: |
---|
316 | outfile.variables['time'][:] = times #Store time relative |
---|
317 | |
---|
318 | if verbose: |
---|
319 | print '------------------------------------------------' |
---|
320 | print 'Statistics:' |
---|
321 | print ' t in [%f, %f], len(t) == %d'\ |
---|
322 | %(min(times.flat), max(times.flat), len(times.flat)) |
---|
323 | |
---|
324 | def store_points(self, |
---|
325 | outfile, |
---|
326 | points_utm, |
---|
327 | elevation, zone=None, new_origin=None, |
---|
328 | points_georeference=None, verbose=False): |
---|
329 | |
---|
330 | """ |
---|
331 | points_utm - currently a list or array of the points in UTM. |
---|
332 | points_georeference - the georeference of the points_utm |
---|
333 | |
---|
334 | How about passing new_origin and current_origin. |
---|
335 | If you get both, do a convertion from the old to the new. |
---|
336 | |
---|
337 | If you only get new_origin, the points are absolute, |
---|
338 | convert to relative |
---|
339 | |
---|
340 | if you only get the current_origin the points are relative, store |
---|
341 | as relative. |
---|
342 | |
---|
343 | if you get no georefs create a new georef based on the minimums of |
---|
344 | points_utm. (Another option would be to default to absolute) |
---|
345 | |
---|
346 | Yes, and this is done in another part of the code. |
---|
347 | Probably geospatial. |
---|
348 | |
---|
349 | If you don't supply either geo_refs, then supply a zone. If not |
---|
350 | the default zone will be used. |
---|
351 | |
---|
352 | precondition: |
---|
353 | header has been called. |
---|
354 | """ |
---|
355 | |
---|
356 | number_of_points = len(points_utm) |
---|
357 | points_utm = array(points_utm) |
---|
358 | |
---|
359 | # given the two geo_refs and the points, do the stuff |
---|
360 | # described in the method header |
---|
361 | points_georeference = ensure_geo_reference(points_georeference) |
---|
362 | new_origin = ensure_geo_reference(new_origin) |
---|
363 | |
---|
364 | if new_origin is None and points_georeference is not None: |
---|
365 | points = points_utm |
---|
366 | geo_ref = points_georeference |
---|
367 | else: |
---|
368 | if new_origin is None: |
---|
369 | new_origin = Geo_reference(zone,min(points_utm[:,0]), |
---|
370 | min(points_utm[:,1])) |
---|
371 | points = new_origin.change_points_geo_ref(points_utm, |
---|
372 | points_georeference) |
---|
373 | geo_ref = new_origin |
---|
374 | |
---|
375 | # At this stage I need a georef and points |
---|
376 | # the points are relative to the georef |
---|
377 | geo_ref.write_NetCDF(outfile) |
---|
378 | |
---|
379 | x = points[:,0] |
---|
380 | y = points[:,1] |
---|
381 | z = outfile.variables['z'][:] |
---|
382 | |
---|
383 | if verbose: |
---|
384 | print '------------------------------------------------' |
---|
385 | print 'More Statistics:' |
---|
386 | print ' Extent (/lon):' |
---|
387 | print ' x in [%f, %f], len(lat) == %d'\ |
---|
388 | %(min(x), max(x), |
---|
389 | len(x)) |
---|
390 | print ' y in [%f, %f], len(lon) == %d'\ |
---|
391 | %(min(y), max(y), |
---|
392 | len(y)) |
---|
393 | print ' z in [%f, %f], len(z) == %d'\ |
---|
394 | %(min(elevation), max(elevation), |
---|
395 | len(elevation)) |
---|
396 | print 'geo_ref: ',geo_ref |
---|
397 | print '------------------------------------------------' |
---|
398 | |
---|
399 | #z = resize(bath_grid,outfile.variables['z'][:].shape) |
---|
400 | #print "points[:,0]", points[:,0] |
---|
401 | outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner() |
---|
402 | outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner() |
---|
403 | outfile.variables['z'][:] = elevation |
---|
404 | outfile.variables['elevation'][:] = elevation #FIXME HACK4 |
---|
405 | |
---|
406 | q = 'elevation' |
---|
407 | # This updates the _range values |
---|
408 | outfile.variables[q+Write_sts.RANGE][0] = min(elevation) |
---|
409 | outfile.variables[q+Write_sts.RANGE][1] = max(elevation) |
---|
410 | |
---|
411 | def store_quantities(self, outfile, sts_precision=Float32, |
---|
412 | slice_index=None, time=None, |
---|
413 | verbose=False, **quant): |
---|
414 | |
---|
415 | """ |
---|
416 | Write the quantity info. |
---|
417 | |
---|
418 | **quant is extra keyword arguments passed in. These must be |
---|
419 | the sts quantities, currently; stage. |
---|
420 | |
---|
421 | if the time array is already been built, use the slice_index |
---|
422 | to specify the index. |
---|
423 | |
---|
424 | Otherwise, use time to increase the time dimension |
---|
425 | |
---|
426 | Maybe make this general, but the viewer assumes these quantities, |
---|
427 | so maybe we don't want it general - unless the viewer is general |
---|
428 | |
---|
429 | precondition: |
---|
430 | triangulation and |
---|
431 | header have been called. |
---|
432 | """ |
---|
433 | if time is not None: |
---|
434 | file_time = outfile.variables['time'] |
---|
435 | slice_index = len(file_time) |
---|
436 | file_time[slice_index] = time |
---|
437 | |
---|
438 | # Write the conserved quantities from Domain. |
---|
439 | # Typically stage, xmomentum, ymomentum |
---|
440 | # other quantities will be ignored, silently. |
---|
441 | # Also write the ranges: stage_range |
---|
442 | for q in Write_sts.sts_quantities: |
---|
443 | if not quant.has_key(q): |
---|
444 | msg = 'STS file can not write quantity %s' %q |
---|
445 | raise NewQuantity, msg |
---|
446 | else: |
---|
447 | q_values = quant[q] |
---|
448 | outfile.variables[q][slice_index] = \ |
---|
449 | q_values.astype(sts_precision) |
---|
450 | |
---|
451 | # This updates the _range values |
---|
452 | q_range = outfile.variables[q+Write_sts.RANGE][:] |
---|
453 | q_values_min = min(q_values) |
---|
454 | if q_values_min < q_range[0]: |
---|
455 | outfile.variables[q+Write_sts.RANGE][0] = q_values_min |
---|
456 | q_values_max = max(q_values) |
---|
457 | if q_values_max > q_range[1]: |
---|
458 | outfile.variables[q+Write_sts.RANGE][1] = q_values_max |
---|
459 | |
---|
460 | def urs2sts(basename_in, basename_out = None, verbose = False, origin = None, |
---|
461 | mean_stage=0.0,zscale=1.0, |
---|
462 | minlat = None, maxlat = None, |
---|
463 | minlon = None, maxlon = None): |
---|
464 | """Convert URS mux2 format for wave propagation to sts format |
---|
465 | |
---|
466 | Specify only basename_in and read files of the form |
---|
467 | out_waveheight-z-mux2 |
---|
468 | |
---|
469 | Also convert latitude and longitude to UTM. All coordinates are |
---|
470 | assumed to be given in the GDA94 datum |
---|
471 | |
---|
472 | origin is a 3-tuple with geo referenced |
---|
473 | UTM coordinates (zone, easting, northing) |
---|
474 | """ |
---|
475 | import os |
---|
476 | from Scientific.IO.NetCDF import NetCDFFile |
---|
477 | from Numeric import Float, Int, Int32, searchsorted, zeros, array |
---|
478 | from Numeric import allclose, around |
---|
479 | |
---|
480 | precision = Float |
---|
481 | |
---|
482 | msg = 'Must use latitudes and longitudes for minlat, maxlon etc' |
---|
483 | |
---|
484 | if minlat != None: |
---|
485 | assert -90 < minlat < 90 , msg |
---|
486 | if maxlat != None: |
---|
487 | assert -90 < maxlat < 90 , msg |
---|
488 | if minlat != None: |
---|
489 | assert maxlat > minlat |
---|
490 | if minlon != None: |
---|
491 | assert -180 < minlon < 180 , msg |
---|
492 | if maxlon != None: |
---|
493 | assert -180 < maxlon < 180 , msg |
---|
494 | if minlon != None: |
---|
495 | assert maxlon > minlon |
---|
496 | |
---|
497 | if basename_out is None: |
---|
498 | stsname = basename_in + '.sts' |
---|
499 | else: |
---|
500 | stsname = basename_out + '.sts' |
---|
501 | |
---|
502 | if (verbose): print 'reading mux2 file' |
---|
503 | times_urs, latitudes_urs, longitudes_urs, elevations, urs_stage = read_binary_mux2(basename_in,verbose) |
---|
504 | |
---|
505 | if (minlat is not None) and (minlon is not None) and (maxlat is not None) and (maxlon is not None): |
---|
506 | latitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),latitudes_urs) |
---|
507 | longitudes = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),longitudes_urs) |
---|
508 | times = compress((latitudes_urs>=minlat)&(latitudes_urs<=maxlat)&(longitudes_urs>=minlon)&(longitudes_urs<=maxlon),times_urs) |
---|
509 | |
---|
510 | |
---|
511 | number_of_points = latitudes.shape[0] |
---|
512 | number_of_times = times.shape[0] |
---|
513 | number_of_latitudes = latitudes.shape[0] |
---|
514 | number_of_longitudes = longitudes.shape[0] |
---|
515 | |
---|
516 | # NetCDF file definition |
---|
517 | outfile = NetCDFFile(stsname, 'w') |
---|
518 | |
---|
519 | description = 'Converted from URS mux2 files: %s'\ |
---|
520 | %(basename_in) |
---|
521 | |
---|
522 | # Create new file |
---|
523 | starttime = times[0] |
---|
524 | sts = Write_sts() |
---|
525 | sts.store_header(outfile, times,number_of_points, description=description, |
---|
526 | verbose=verbose,sts_precision=Float) |
---|
527 | |
---|
528 | # Store |
---|
529 | from anuga.coordinate_transforms.redfearn import redfearn |
---|
530 | x = zeros(number_of_points, Float) #Easting |
---|
531 | y = zeros(number_of_points, Float) #Northing |
---|
532 | |
---|
533 | # Check zone boundaries |
---|
534 | refzone, _, _ = redfearn(latitudes[0],longitudes[0]) |
---|
535 | |
---|
536 | for i in range(number_of_points): |
---|
537 | zone, easting, northing = redfearn(latitudes[i],longitudes[i]) |
---|
538 | x[i] = easting |
---|
539 | y[i] = northing |
---|
540 | #print zone,easting,northing |
---|
541 | |
---|
542 | if origin is None: |
---|
543 | origin = Geo_reference(refzone,min(x),min(y)) |
---|
544 | geo_ref = write_NetCDF_georeference(origin, outfile) |
---|
545 | |
---|
546 | z = elevations |
---|
547 | |
---|
548 | #print geo_ref.get_xllcorner() |
---|
549 | #print geo_ref.get_yllcorner() |
---|
550 | |
---|
551 | z = resize(z,outfile.variables['z'][:].shape) |
---|
552 | outfile.variables['x'][:] = x - geo_ref.get_xllcorner() |
---|
553 | outfile.variables['y'][:] = y - geo_ref.get_yllcorner() |
---|
554 | outfile.variables['z'][:] = z #FIXME HACK for bacwards compat. |
---|
555 | outfile.variables['elevation'][:] = z |
---|
556 | |
---|
557 | stage = outfile.variables['stage'] |
---|
558 | |
---|
559 | #print urs_stage.shape |
---|
560 | for j in range(len(times)): |
---|
561 | for i in range(number_of_points): |
---|
562 | w = zscale*urs_stage[i,j] + mean_stage |
---|
563 | stage[j,i] = w |
---|
564 | |
---|
565 | outfile.close() |
---|
566 | |
---|
567 | def read_sts(filename): |
---|
568 | from Scientific.IO.NetCDF import NetCDFFile |
---|
569 | from Numeric import asarray |
---|
570 | fid = NetCDFFile(filename+'.sts', 'r') #Open existing file for read |
---|
571 | time = fid.variables['time'] #Timesteps |
---|
572 | |
---|
573 | # Get the variables as Numeric arrays |
---|
574 | x = fid.variables['x'][:] #x-coordinates of vertices |
---|
575 | y = fid.variables['y'][:] #y-coordinates of vertices |
---|
576 | elevation = fid.variables['elevation'] #Elevation |
---|
577 | stage = fid.variables['stage'] #Water level |
---|
578 | |
---|
579 | starttime = fid.starttime[0] |
---|
580 | coordinates=transpose(asarray([x.tolist(),y.tolist()])) |
---|
581 | |
---|
582 | conserved_quantities = [] |
---|
583 | other_quantities = [] |
---|
584 | for quantity in fid.variables.keys(): |
---|
585 | dimensions = fid.variables[quantity].dimensions |
---|
586 | if 'number_of_timesteps' in dimensions: |
---|
587 | conserved_quantities.append(quantity) |
---|
588 | else: other_quantities.append(quantity) |
---|
589 | |
---|
590 | other_quantities.remove('x') |
---|
591 | other_quantities.remove('y') |
---|
592 | other_quantities.remove('z') |
---|
593 | try: |
---|
594 | other_quantities.remove('stage_range') |
---|
595 | other_quantities.remove('elevation_range') |
---|
596 | except: |
---|
597 | pass |
---|
598 | |
---|
599 | from pylab import plot,show |
---|
600 | plot(x/1000,y/1000,'o') |
---|
601 | show() |
---|
602 | |
---|
603 | |
---|
604 | #def plot_sts_gauges(): |
---|
605 | |
---|
606 | |
---|
607 | file_in = "big_out_waveheight-z-mux2" |
---|
608 | file_out="big" |
---|
609 | |
---|
610 | east=99 |
---|
611 | west=98 |
---|
612 | north=2 |
---|
613 | south=0 |
---|
614 | urs2sts(file_in,file_out,verbose=False,minlat=south,maxlat=north,minlon=west,maxlon=east) |
---|
615 | read_sts(file_out) |
---|