source: anuga_core/source/anuga/file_conversion/file_conversion.py @ 7742

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

Further split up file conversion to modularise shallow_water module.

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