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

Last change on this file since 8780 was 8780, checked in by steve, 11 years ago

Some changes to allow netcdf4 use

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 anuga.file.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
179    nan_ua = file_u.variables['UA'].missing_value
180    nan_va = file_v.variables['VA'].missing_value
181    if hasattr(file_e.variables[zname],'missing_value'):
182        nan_e  = file_e.variables[zname].missing_value
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.