source: trunk/anuga_core/source/anuga/file_conversion/ferret2sww.py @ 7814

Last change on this file since 7814 was 7814, checked in by hudson, 14 years ago

New filename conventions for file conversion. Filenames must always be passed in with the correct extension.

File size: 15.3 KB
Line 
1"""
2    Convert a ferret file to an SWW file.
3"""
4# external modules
5import numpy as num
6
7
8# ANUGA modules
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_float
10from anuga.file.sww import Write_sww 
11from anuga.coordinate_transforms.geo_reference import Geo_reference, \
12     write_NetCDF_georeference                         
13import anuga.utilities.log as log
14
15#local modules
16from anuga.file_conversion.file_conversion import get_min_max_indices                           
17
18
19def ferret2sww(basename_in, name_out=None,
20               verbose=False,
21               minlat=None, maxlat=None,
22               minlon=None, maxlon=None,
23               mint=None, maxt=None, mean_stage=0,
24               origin=None, zscale=1,
25               fail_on_NaN=True,
26               NaN_filler=0,
27               elevation=None,
28               inverted_bathymetry=True
29               ): #FIXME: Bathymetry should be obtained
30                                  #from MOST somehow.
31                                  #Alternatively from elsewhere
32                                  #or, as a last resort,
33                                  #specified here.
34                                  #The value of -100 will work
35                                  #for the Wollongong tsunami
36                                  #scenario but is very hacky
37    """Convert MOST and 'Ferret' NetCDF format for wave propagation to
38    sww format native to abstract_2d_finite_volumes.
39
40    Specify only basename_in and read files of the form
41    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
42    relative height, x-velocity and y-velocity, respectively.
43
44    Also convert latitude and longitude to UTM. All coordinates are
45    assumed to be given in the GDA94 datum.
46
47    min's and max's: If omitted - full extend is used.
48    To include a value min may equal it, while max must exceed it.
49    Lat and lon are assuemd to be in decimal degrees
50
51    origin is a 3-tuple with geo referenced
52    UTM coordinates (zone, easting, northing)
53
54    nc format has values organised as HA[TIME, LATITUDE, LONGITUDE]
55    which means that longitude is the fastest
56    varying dimension (row major order, so to speak)
57
58    ferret2sww uses grid points as vertices in a triangular grid
59    counting vertices from lower left corner upwards, then right
60    """
61
62    from Scientific.IO.NetCDF import NetCDFFile
63
64    _assert_lat_long(minlat, maxlat, minlon, maxlon)
65
66    if name_out != None and name_out[-4:] != '.sww':
67        raise IOError('Output file %s should be of type .sww.' % name_out)
68
69    # Get NetCDF data
70    if verbose: log.critical('Reading files %s_*.nc' % basename_in)
71
72    # Wave amplitude (cm)
73    file_h = NetCDFFile(basename_in + '_ha.nc', netcdf_mode_r) 
74   
75    # Velocity (x) (cm/s)
76    file_u = NetCDFFile(basename_in + '_ua.nc', netcdf_mode_r)
77     
78    # Velocity (y) (cm/s)
79    file_v = NetCDFFile(basename_in + '_va.nc', netcdf_mode_r)
80   
81    # Elevation (z) (m)
82    file_e = NetCDFFile(basename_in + '_e.nc', netcdf_mode_r) 
83
84    if name_out is None:
85        swwname = basename_in + '.sww'
86    else:
87        swwname = name_out
88
89    # Get dimensions of file_h
90    for dimension in file_h.dimensions.keys():
91        if dimension[:3] == 'LON':
92            dim_h_longitude = dimension
93        if dimension[:3] == 'LAT':
94            dim_h_latitude = dimension
95        if dimension[:4] == 'TIME':
96            dim_h_time = dimension
97
98    times = file_h.variables[dim_h_time]
99    latitudes = file_h.variables[dim_h_latitude]
100    longitudes = file_h.variables[dim_h_longitude]
101
102    kmin, kmax, lmin, lmax = get_min_max_indices(latitudes[:],
103                                                  longitudes[:],
104                                                  minlat, maxlat,
105                                                  minlon, maxlon)
106    # get dimensions for file_e
107    for dimension in file_e.dimensions.keys():
108        if dimension[:3] == 'LON':
109            dim_e_longitude = dimension
110        if dimension[:3] == 'LAT':
111            dim_e_latitude = dimension
112
113    # get dimensions for file_u
114    for dimension in file_u.dimensions.keys():
115        if dimension[:3] == 'LON':
116            dim_u_longitude = dimension
117        if dimension[:3] == 'LAT':
118            dim_u_latitude = dimension
119
120    # get dimensions for file_v
121    for dimension in file_v.dimensions.keys():
122        if dimension[:3] == 'LON':
123            dim_v_longitude = dimension
124        if dimension[:3] == 'LAT':
125            dim_v_latitude = dimension
126
127    # Precision used by most for lat/lon is 4 or 5 decimals
128    e_lat = num.around(file_e.variables[dim_e_latitude][:], 5)
129    e_lon = num.around(file_e.variables[dim_e_longitude][:], 5)
130
131    # Check that files are compatible
132    assert num.allclose(latitudes, file_u.variables[dim_u_latitude])
133    assert num.allclose(latitudes, file_v.variables[dim_v_latitude])
134    assert num.allclose(latitudes, e_lat)
135
136    assert num.allclose(longitudes, file_u.variables[dim_u_longitude])
137    assert num.allclose(longitudes, file_v.variables[dim_v_longitude])
138    assert num.allclose(longitudes, e_lon)
139
140    if mint is None:
141        jmin = 0
142        mint = times[0]
143    else:
144        jmin = num.searchsorted(times, mint)
145       
146        # numpy.int32 didn't work in slicing of amplitude below
147        jmin = int(jmin)
148
149    if maxt is None:
150        jmax = len(times)
151        maxt = times[-1]
152    else:
153        jmax = num.searchsorted(times, maxt)
154       
155        # numpy.int32 didn't work in slicing of amplitude below
156        jmax = int(jmax)       
157
158    kmin, kmax, lmin, lmax = get_min_max_indices(latitudes[:],
159                                                  longitudes[:],
160                                                  minlat, maxlat,
161                                                  minlon, maxlon)
162
163
164    times = times[jmin:jmax]
165    latitudes = latitudes[kmin:kmax]
166    longitudes = longitudes[lmin:lmax]
167
168    if verbose: log.critical('cropping')
169
170    zname = 'ELEVATION'
171
172    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
173    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
174    vspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lat
175    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
176
177    # Get missing values
178    nan_ha = file_h.variables['HA'].missing_value[0]
179    nan_ua = file_u.variables['UA'].missing_value[0]
180    nan_va = file_v.variables['VA'].missing_value[0]
181    if hasattr(file_e.variables[zname],'missing_value'):
182        nan_e  = file_e.variables[zname].missing_value[0]
183    else:
184        nan_e = None
185
186    # Cleanup
187    missing = (amplitudes == nan_ha)
188    if num.sometrue (missing):
189        if fail_on_NaN:
190            msg = 'NetCDFFile %s contains missing values' \
191                  % basename_in + '_ha.nc'
192            raise DataMissingValuesError, msg
193        else:
194            amplitudes = amplitudes*(missing==0) + missing*NaN_filler
195
196    missing = (uspeed == nan_ua)
197    if num.sometrue (missing):
198        if fail_on_NaN:
199            msg = 'NetCDFFile %s contains missing values' \
200                  % basename_in + '_ua.nc'
201            raise DataMissingValuesError, msg
202        else:
203            uspeed = uspeed*(missing==0) + missing*NaN_filler
204
205    missing = (vspeed == nan_va)
206    if num.sometrue (missing):
207        if fail_on_NaN:
208            msg = 'NetCDFFile %s contains missing values' \
209                  % basename_in + '_va.nc'
210            raise DataMissingValuesError, msg
211        else:
212            vspeed = vspeed*(missing==0) + missing*NaN_filler
213
214    missing = (elevations == nan_e)
215    if num.sometrue (missing):
216        if fail_on_NaN:
217            msg = 'NetCDFFile %s contains missing values' \
218                  % basename_in + '_e.nc'
219            raise DataMissingValuesError, msg
220        else:
221            elevations = elevations*(missing==0) + missing*NaN_filler
222
223    number_of_times = times.shape[0]
224    number_of_latitudes = latitudes.shape[0]
225    number_of_longitudes = longitudes.shape[0]
226
227    assert amplitudes.shape[0] == number_of_times
228    assert amplitudes.shape[1] == number_of_latitudes
229    assert amplitudes.shape[2] == number_of_longitudes
230
231    if verbose:
232        _show_stats((latitudes, longitudes), times, amplitudes, \
233                    (uspeed, vspeed), elevations)
234
235    # print number_of_latitudes, number_of_longitudes
236    number_of_points = number_of_latitudes * number_of_longitudes
237    number_of_volumes = (number_of_latitudes-1) * (number_of_longitudes-1) * 2
238
239    file_h.close()
240    file_u.close()
241    file_v.close()
242    file_e.close()
243
244    # NetCDF file definition
245    outfile = NetCDFFile(swwname, netcdf_mode_w)
246
247    description = 'Converted from Ferret files: %s, %s, %s, %s' \
248                  % (basename_in + '_ha.nc',
249                     basename_in + '_ua.nc',
250                     basename_in + '_va.nc',
251                     basename_in + '_e.nc')
252
253    # Create new file
254    starttime = times[0]
255
256    sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
257    sww.store_header(outfile, times, number_of_volumes,
258                     number_of_points, description=description,
259                     verbose=verbose, sww_precision=netcdf_float)
260
261    # Store
262    from anuga.coordinate_transforms.redfearn import redfearn
263    x = num.zeros(number_of_points, num.float)  #Easting
264    y = num.zeros(number_of_points, num.float)  #Northing
265
266    if verbose:
267        log.critical('Making triangular grid')
268
269    # Check zone boundaries
270    refzone, _, _ = redfearn(latitudes[0], longitudes[0])
271
272    vertices = {}
273    i = 0
274    for k, lat in enumerate(latitudes):       # Y direction
275        for l, lon in enumerate(longitudes):  # X direction
276            vertices[l, k] = i
277
278            _, easting, northing = redfearn(lat, lon)
279
280            #msg = 'Zone boundary crossed at longitude =', lon
281            #assert zone == refzone, msg
282            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
283            x[i] = easting
284            y[i] = northing
285            i += 1
286
287    #Construct 2 triangles per 'rectangular' element
288    volumes = []
289    for l in range(number_of_longitudes-1):    # X direction
290        for k in range(number_of_latitudes-1): # Y direction
291            v1 = vertices[l, k+1]
292            v2 = vertices[l, k]
293            v3 = vertices[l+1, k+1]
294            v4 = vertices[l+1, k]
295
296            volumes.append([v1, v2, v3]) #Upper element
297            volumes.append([v4, v3, v2]) #Lower element
298
299    volumes = num.array(volumes, num.int)      #array default#
300
301    if origin is None:
302        origin = Geo_reference(refzone, min(x), min(y))
303    geo_ref = write_NetCDF_georeference(origin, outfile)
304
305    if elevation is not None:
306        z = elevation
307    else:
308        if inverted_bathymetry:
309            z = -1 * elevations
310        else:
311            z = elevations
312    #FIXME: z should be obtained from MOST and passed in here
313
314    #FIXME use the Write_sww instance(sww) to write this info
315    z = num.resize(z, outfile.variables['elevation'][:].shape)
316    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
317    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
318    #outfile.variables['z'][:] = z             #FIXME HACK for bacwards compat.
319    outfile.variables['elevation'][:] = z
320    outfile.variables['volumes'][:] = volumes.astype(num.int32) #For Opteron 64
321
322    #Time stepping
323    stage = outfile.variables['stage']
324    xmomentum = outfile.variables['xmomentum']
325    ymomentum = outfile.variables['ymomentum']
326
327    if verbose:
328        log.critical('Converting quantities')
329
330    n = len(times)
331    for j in range(n):
332        if verbose and j % ((n+10)/10) == 0:
333            log.critical('  Doing %d of %d' % (j, n))
334
335        i = 0
336        for k in range(number_of_latitudes):      # Y direction
337            for l in range(number_of_longitudes): # X direction
338                w = zscale * amplitudes[j, k, l] / 100 + mean_stage
339                stage[j, i] = w
340                h = w - z[i]
341                xmomentum[j, i] = uspeed[j, k, l]/100*h
342                ymomentum[j, i] = vspeed[j, k, l]/100*h
343                i += 1
344
345    #outfile.close()
346
347    #FIXME: Refactor using code from file_function.statistics
348    #Something like print swwstats(swwname)
349    if verbose:
350        time_info = times, starttime, mint, maxt
351        _show_sww_stats(outfile, swwname, geo_ref, time_info)
352
353    outfile.close()
354
355
356def _show_stats(latlong, times, amplitudes, speeds, elevations):
357    """ Print the statistics nicely to the log file """
358   
359    latitudes, longitudes = latlong
360    uspeed, vspeed = speeds
361   
362    log.critical('------------------------------------------------')
363    log.critical('Statistics:')
364    log.critical('  Extent (lat/lon):')
365    log.critical('    lat in [%f, %f], len(lat) == %d'
366                 % (num.min(latitudes), num.max(latitudes),
367                    len(latitudes.flat)))
368    log.critical('    lon in [%f, %f], len(lon) == %d'
369                 % (num.min(longitudes), num.max(longitudes),
370                    len(longitudes.flat)))
371    log.critical('    t in [%f, %f], len(t) == %d'
372                 % (num.min(times), num.max(times), len(times.flat)))
373
374    name = 'Amplitudes (ha) [cm]'
375    log.critical(%s in [%f, %f]'
376                 % (name, num.min(amplitudes), num.max(amplitudes)))
377
378    name = 'Speeds (ua) [cm/s]'
379    log.critical(%s in [%f, %f]'
380                 % (name, num.min(uspeed), num.max(uspeed)))
381
382    name = 'Speeds (va) [cm/s]'
383    log.critical(%s in [%f, %f]'
384                 % (name, num.min(vspeed), num.max(vspeed)))
385
386    name = 'Elevations (e) [m]'
387    log.critical(%s in [%f, %f]'
388                 % (name, num.min(elevations), num.max(elevations)))
389                 
390             
391
392def _show_sww_stats(outfile, swwname, geo_ref, time_info):         
393    """ Log SWW output stats. """
394    times, starttime, mint, maxt = time_info
395    x = outfile.variables['x'][:]
396    y = outfile.variables['y'][:]
397    log.critical('------------------------------------------------')
398    log.critical('Statistics of output file:')
399    log.critical('  Name: %s' %swwname)
400    log.critical('  Reference:')
401    log.critical('    Lower left corner: [%f, %f]'
402                 % (geo_ref.get_xllcorner(), geo_ref.get_yllcorner()))
403    log.critical(' Start time: %f' %starttime)
404    log.critical('    Min time: %f' %mint)
405    log.critical('    Max time: %f' %maxt)
406    log.critical('  Extent:')
407    log.critical('    x [m] in [%f, %f], len(x) == %d'
408                 % (num.min(x), num.max(x), len(x.flat)))
409    log.critical('    y [m] in [%f, %f], len(y) == %d'
410                 % (num.min(y), num.max(y), len(y.flat)))
411    log.critical('    t [s] in [%f, %f], len(t) == %d'
412                 % (min(times), max(times), len(times)))
413    log.critical('  Quantities [SI units]:')
414    for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
415        q = outfile.variables[name][:]    # .flatten()
416        log.critical('    %s in [%f, %f]' % \
417                        (name, num.min(q), num.max(q)))       
418                 
419
420def _assert_lat_long(minlat, maxlat, minlon, maxlon):
421    """Check latitudes and longitudes for validity."""
422   
423    msg = 'Must use latitudes and longitudes for minlat, maxlon etc'
424
425    if minlat != None:
426        assert -90 < minlat < 90 , msg
427    if maxlat != None:
428        assert -90 < maxlat < 90 , msg
429        if minlat != None:
430            assert maxlat > minlat
431    if minlon != None:
432        assert -180 < minlon < 180 , msg
433    if maxlon != None:
434        assert -180 < maxlon < 180 , msg
435        if minlon != None:
436            assert maxlon > minlon
Note: See TracBrowser for help on using the repository browser.