source: trunk/anuga_work/development/boxingday08/urs.py @ 7884

Last change on this file since 7884 was 5265, checked in by jakeman, 16 years ago

added file to read and write spatial timeseries (.sts) files

File size: 21.5 KB
Line 
1import os
2from struct import unpack
3from anuga.utilities.anuga_exceptions import ANUGAError
4from Numeric import Float, Int, zeros
5from 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
10from anuga.utilities.numerical_tools import ensure_numeric
11from anuga.config import max_float
12from anuga.coordinate_transforms.geo_reference import Geo_reference, \
13write_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"""
20def 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
227class 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
460def 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
567def 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
607file_in = "big_out_waveheight-z-mux2"
608file_out="big"
609
610east=99
611west=98
612north=2
613south=0
614urs2sts(file_in,file_out,verbose=False,minlat=south,maxlat=north,minlon=west,maxlon=east)
615read_sts(file_out)
Note: See TracBrowser for help on using the repository browser.