source: trunk/anuga_core/source/anuga/file/test_mux.py @ 9213

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

Some changes to allow netcdf4 use
File size: 60.8 KB
Line 
1import unittest
2import tempfile
3import numpy as num
4import os
5from struct import pack, unpack
6from anuga.file.netcdf import NetCDFFile
7
8from anuga.utilities.numerical_tools import ensure_numeric
9from anuga.coordinate_transforms.redfearn import redfearn
10from anuga.coordinate_transforms.geo_reference import Geo_reference
11
12from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
13                            NORTH_VELOCITY_LABEL
14
15from mux import WAVEHEIGHT_MUX2_LABEL, EAST_VELOCITY_MUX2_LABEL, \
16                NORTH_VELOCITY_MUX2_LABEL
17               
18from mux import read_mux2_py
19from anuga.file_conversion.urs2sts import urs2sts
20from anuga.file.urs import Read_urs
21
22class Test_Mux(unittest.TestCase):
23    def setUp(self):
24        pass
25
26    def tearDown(self):
27        pass
28
29    def write_mux(self, lat_long_points, time_step_count, time_step,
30                  depth=None, ha=None, ua=None, va=None):
31        """
32        This will write 3 non-gridded mux files, for testing.
33        If no quantities are passed in,
34        na and va quantities will be the Easting values.
35        Depth and ua will be the Northing value.
36       
37        The mux file format has south as positive so
38        this function will swap the sign for va. 
39        """
40
41        #print "lat_long_points", lat_long_points
42        #print "time_step_count",time_step_count
43        #print "time_step",
44
45       
46        points_num = len(lat_long_points)
47        lonlatdeps = []
48        quantities = ['HA','UA','VA']
49       
50        mux_names = [WAVEHEIGHT_MUX_LABEL,
51                     EAST_VELOCITY_LABEL,
52                     NORTH_VELOCITY_LABEL]
53        quantities_init = [[],[],[]]
54        # urs binary is latitude fastest
55        for point in lat_long_points:
56            lat = point[0]
57            lon = point[1]
58            _ , e, n = redfearn(lat, lon)
59            if depth is None:
60                this_depth = n
61            else:
62                this_depth = depth
63            if ha is None:
64                this_ha = e
65            else:
66                this_ha = ha
67            if ua is None:
68                this_ua = n
69            else:
70                this_ua = ua
71            if va is None:
72                this_va = e   
73            else:
74                this_va = va         
75            lonlatdeps.append([lon, lat, this_depth])
76            quantities_init[0].append(this_ha) # HA
77            quantities_init[1].append(this_ua) # UA
78            quantities_init[2].append(this_va) # VA
79               
80        file_handle, base_name = tempfile.mkstemp("")
81        os.close(file_handle)
82        os.remove(base_name)
83
84        files = []       
85        for i, q in enumerate(quantities): 
86            quantities_init[i] = ensure_numeric(quantities_init[i])
87            #print "HA_init", HA_init
88            q_time = num.zeros((time_step_count, points_num), num.float)
89            for time in range(time_step_count):
90                q_time[time,:] = quantities_init[i] #* time * 4
91           
92            #Write C files
93            columns = 3 # long, lat , depth
94            file = base_name + mux_names[i]
95            #print "base_name file",file
96            f = open(file, 'wb')
97            files.append(file)
98            f.write(pack('i',points_num))
99            f.write(pack('i',time_step_count))
100            f.write(pack('f',time_step))
101
102            #write lat/long info
103            for lonlatdep in lonlatdeps:
104                for float in lonlatdep:
105                    f.write(pack('f',float))
106                   
107            # Write quantity info
108            for time in  range(time_step_count):
109                for point_i in range(points_num):
110                    f.write(pack('f',q_time[time,point_i]))
111                    #print " mux_names[i]", mux_names[i]
112                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
113            f.close()
114        return base_name, files
115       
116   
117    def delete_mux(self, files):
118        for file in files:
119            try:
120                os.remove(file)
121            except:
122                pass
123       
124    def write_mux2(self, lat_long_points, time_step_count, time_step,
125                   first_tstep, last_tstep,
126                   depth=None, ha=None, ua=None, va=None):
127        """
128        This will write 3 non-gridded mux files, for testing.
129        If no quantities are passed in,
130        na and va quantities will be the Easting values.
131        Depth and ua will be the Northing value.
132        """
133        #print "lat_long_points", lat_long_points
134        #print "time_step_count",time_step_count
135        #print "time_step",
136
137        #irrelevant header information
138        ig=ilon=ilat=0
139        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
140
141        points_num = len(lat_long_points)
142        latlondeps = []
143        quantities = ['HA','UA','VA']
144
145        mux_names = [WAVEHEIGHT_MUX2_LABEL,
146                     EAST_VELOCITY_MUX2_LABEL,
147                     NORTH_VELOCITY_MUX2_LABEL]
148
149        msg='first_tstep and last_step arrays must have same length as number of points'
150        assert len(first_tstep)==points_num,msg
151        assert len(last_tstep)==points_num,msg
152
153        if depth is not None:
154            depth=ensure_numeric(depth)
155            assert len(depth)==points_num
156        if ha is not None:
157            ha=ensure_numeric(ha)
158            assert ha.shape==(points_num,time_step_count)
159        if ua is not None:
160            ua=ensure_numeric(ua)
161            assert ua.shape==(points_num,time_step_count)
162        if va is not None:
163            va=ensure_numeric(va)
164            assert va.shape==(points_num,time_step_count)
165
166        quantities_init = [[],[],[]]
167        # urs binary is latitude fastest
168        for i,point in enumerate(lat_long_points):
169            lat = point[0]
170            lon = point[1]
171            _ , e, n = redfearn(lat, lon)
172            if depth is None:
173                this_depth = n
174            else:
175                this_depth = depth[i]
176            latlondeps.append([lat, lon, this_depth])
177
178            if ha is None:
179                this_ha = e
180                quantities_init[0].append(num.ones(time_step_count,num.float)*this_ha) # HA
181            else:
182                quantities_init[0].append(ha[i])
183            if ua is None:
184                this_ua = n
185                quantities_init[1].append(num.ones(time_step_count,num.float)*this_ua) # UA
186            else:
187                quantities_init[1].append(ua[i])
188            if va is None:
189                this_va = e
190                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
191            else:
192                quantities_init[2].append(-va[i]) # South is negative in MUX
193
194        file_handle, base_name = tempfile.mkstemp("write_mux2")
195        os.close(file_handle)
196        os.remove(base_name)
197
198        files = []       
199        for i, q in enumerate(quantities):
200            q_time = num.zeros((time_step_count, points_num), num.float)
201            quantities_init[i] = ensure_numeric(quantities_init[i])
202            for time in range(time_step_count):
203                #print i, q, time, quantities_init[i][:,time]
204                q_time[time,:] = quantities_init[i][:,time]
205                #print i, q, time, q_time[time, :]               
206
207            #Write C files
208            columns = 3 # long, lat , depth
209            file = base_name + mux_names[i]
210           
211            #print 'base_name file', file
212            f = open(file, 'wb')
213            files.append(file)
214
215            f.write(pack('i',points_num))
216            #write mux 2 header
217            for latlondep in latlondeps:
218                f.write(pack('f',latlondep[0]))
219                f.write(pack('f',latlondep[1]))
220                f.write(pack('f',mcolat))
221                f.write(pack('f',mcolon))
222                f.write(pack('i',ig))
223                f.write(pack('i',ilon))
224                f.write(pack('i',ilat))
225                f.write(pack('f',latlondep[2]))
226                f.write(pack('f',centerlat))
227                f.write(pack('f',centerlon))
228                f.write(pack('f',offset))
229                f.write(pack('f',az))
230                f.write(pack('f',baz))
231                f.write(pack('f',time_step))
232                f.write(pack('i',time_step_count))
233                for j in range(4): # identifier
234                    f.write(pack('f',id))   
235
236            #first_tstep=1
237            #last_tstep=time_step_count
238            for i,latlondep in enumerate(latlondeps):
239                f.write(pack('i',first_tstep[i]))
240            for i,latlondep in enumerate(latlondeps):
241                f.write(pack('i',last_tstep[i]))
242
243            # Find when first station starts recording
244            min_tstep = min(first_tstep)
245            # Find when all stations have stopped recording
246            max_tstep = max(last_tstep)
247
248            #for time in  range(time_step_count):
249            for time in range(min_tstep-1,max_tstep):
250                f.write(pack('f',time*time_step))               
251                for point_i in range(points_num):
252                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
253                        #print 'writing', time, point_i, q_time[time, point_i]
254                        f.write(pack('f', q_time[time, point_i]))
255            f.close()
256
257        return base_name, files
258
259
260    def test_urs2sts_read_mux2_pyI(self):
261        """test_urs2sts_read_mux2_pyI(self):
262        Constant stage,momentum at each gauge
263        """
264        tide = 1
265        time_step_count = 3
266        time_step = 2
267        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
268        n=len(lat_long_points)
269        first_tstep=num.ones(n,num.int)
270        last_tstep=time_step_count*num.ones(n,num.int)
271        depth=20*num.ones(n,num.float)
272        ha=2*num.ones((n,time_step_count),num.float)
273        ua=5*num.ones((n,time_step_count),num.float)
274        va=-10*num.ones((n,time_step_count),num.float)
275        #-ve added to take into account mux file format where south is positive.
276        base_name, files = self.write_mux2(lat_long_points,
277                                      time_step_count, time_step,
278                                      first_tstep, last_tstep,
279                                      depth=depth,
280                                      ha=ha,
281                                      ua=ua,
282                                      va=va)
283
284        weights=num.ones(1, num.float)
285        #ensure that files are indeed mux2 files
286        times, latitudes, longitudes, elevation, stage, starttime = read_mux2_py([files[0]],weights)
287        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
288        msg='ha and ua have different gauge meta data'
289        assert num.allclose(times,ua_times) and num.allclose(latitudes,ua_latitudes) and num.allclose(longitudes,ua_longitudes) and num.allclose(elevation,ua_elevation) and num.allclose(starttime,starttime_ua), msg
290        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity, starttime_va=read_mux2_py([files[2]],weights)
291        msg='ha and va have different gauge meta data'
292        assert num.allclose(times,va_times) and num.allclose(latitudes,va_latitudes) and num.allclose(longitudes,va_longitudes) and num.allclose(elevation,va_elevation) and num.allclose(starttime,starttime_va), msg
293
294        self.delete_mux(files)
295
296        msg='time array has incorrect length'
297        assert times.shape[0]==time_step_count,msg
298       
299        msg = 'time array is incorrect'
300        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
301        assert num.allclose(times,time_step*num.arange(time_step_count)), msg
302       
303        msg='Incorrect gauge positions returned'
304        for i,point in enumerate(lat_long_points):
305            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
306
307        msg='Incorrect gauge depths returned'
308        assert num.allclose(elevation,-depth),msg
309        msg='incorrect gauge height time series returned'
310        assert num.allclose(stage,ha)
311        msg='incorrect gauge ua time series returned'
312        assert num.allclose(xvelocity,ua)
313        msg='incorrect gauge va time series returned'
314        assert num.allclose(yvelocity, -va)
315
316    def test_urs2sts_read_mux2_pyII(self):
317        """Spatially varing stage
318        """
319        tide = 1
320        time_step_count = 3
321        time_step = 2
322        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
323        n=len(lat_long_points)
324        first_tstep=num.ones(n,num.int)
325        last_tstep=(time_step_count)*num.ones(n,num.int)
326        depth=20*num.ones(n,num.float)
327        ha=2*num.ones((n,time_step_count),num.float)
328        ha[0]=num.arange(0,time_step_count)+1
329        ha[1]=time_step_count-num.arange(1,time_step_count+1)
330        ha[1]=num.arange(time_step_count,2*time_step_count)
331        ha[2]=num.arange(2*time_step_count,3*time_step_count)
332        ha[3]=num.arange(3*time_step_count,4*time_step_count)
333        ua=5*num.ones((n,time_step_count),num.float)
334        va=-10*num.ones((n,time_step_count),num.float)
335        #-ve added to take into account mux file format where south is positive.
336        base_name, files = self.write_mux2(lat_long_points,
337                                      time_step_count, time_step,
338                                      first_tstep, last_tstep,
339                                      depth=depth,
340                                      ha=ha,
341                                      ua=ua,
342                                      va=va)
343
344        weights=num.ones(1, num.float)
345        #ensure that files are indeed mux2 files
346        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
347        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
348        msg='ha and ua have different gauge meta data'
349        assert num.allclose(times,ua_times) and num.allclose(latitudes,ua_latitudes) and num.allclose(longitudes,ua_longitudes) and num.allclose(elevation,ua_elevation) and num.allclose(starttime,starttime_ua), msg
350        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
351        msg='ha and va have different gauge meta data'
352        assert num.allclose(times,va_times) and num.allclose(latitudes,va_latitudes) and num.allclose(longitudes,va_longitudes) and num.allclose(elevation,va_elevation) and num.allclose(starttime,starttime_va), msg
353
354
355        self.delete_mux(files)
356
357        msg='time array has incorrect length'
358        #assert times.shape[0]==time_step_count,msg
359        msg = 'time array is incorrect'
360        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
361        msg='Incorrect gauge positions returned'
362        for i,point in enumerate(lat_long_points):
363            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
364
365        msg='Incorrect gauge depths returned'
366        assert num.allclose(elevation, -depth),msg
367        msg='incorrect gauge height time series returned'
368        assert num.allclose(stage, ha)
369        msg='incorrect gauge ua time series returned'
370        assert num.allclose(xvelocity, ua)
371        msg='incorrect gauge va time series returned'
372        assert num.allclose(yvelocity, -va) # South is positive in MUX
373
374    def test_urs2sts_read_mux2_pyIII(self):
375        """Varying start and finish times
376        """
377        tide = 1
378        time_step_count = 3
379        time_step = 2
380        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
381        n=len(lat_long_points)
382        first_tstep=num.ones(n,num.int)
383        first_tstep[0]+=1
384        first_tstep[2]+=1
385        last_tstep=(time_step_count)*num.ones(n,num.int)
386        last_tstep[0]-=1
387
388        depth=20*num.ones(n,num.float)
389        ha=2*num.ones((n,time_step_count),num.float)
390        ha[0]=num.arange(0,time_step_count)
391        ha[1]=num.arange(time_step_count,2*time_step_count)
392        ha[2]=num.arange(2*time_step_count,3*time_step_count)
393        ha[3]=num.arange(3*time_step_count,4*time_step_count)
394        ua=5*num.ones((n,time_step_count),num.float)
395        va=-10*num.ones((n,time_step_count),num.float)
396        #-ve added to take into account mux file format where south is positive.
397        base_name, files = self.write_mux2(lat_long_points,
398                                      time_step_count, time_step,
399                                      first_tstep, last_tstep,
400                                      depth=depth,
401                                      ha=ha,
402                                      ua=ua,
403                                      va=va)
404
405        weights=num.ones(1, num.float)
406        #ensure that files are indeed mux2 files
407        times, latitudes, longitudes, elevation, stage, starttime=read_mux2_py([files[0]],weights)
408        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity, starttime_ua=read_mux2_py([files[1]],weights)
409        msg='ha and ua have different gauge meta data'
410        assert num.allclose(times,ua_times) and num.allclose(latitudes,ua_latitudes) and num.allclose(longitudes,ua_longitudes) and num.allclose(elevation,ua_elevation) and num.allclose(starttime,starttime_ua), msg
411        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
412        msg='ha and va have different gauge meta data'
413        assert num.allclose(times,va_times) and num.allclose(latitudes,va_latitudes) and num.allclose(longitudes,va_longitudes) and num.allclose(elevation,va_elevation) and num.allclose(starttime,starttime_va), msg
414
415        self.delete_mux(files)
416
417        msg='time array has incorrect length'
418        #assert times.shape[0]==time_step_count,msg
419        msg = 'time array is incorrect'
420        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
421        msg='Incorrect gauge positions returned'
422        for i,point in enumerate(lat_long_points):
423            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
424
425
426        # Set original data used to write mux file to be zero when gauges are
427        #not recdoring
428        ha[0][0]=0.0
429        ha[0][time_step_count-1]=0.0;
430        ha[2][0]=0.0;
431        ua[0][0]=0.0
432        ua[0][time_step_count-1]=0.0;
433        ua[2][0]=0.0;
434        va[0][0]=0.0
435        va[0][time_step_count-1]=0.0;
436        va[2][0]=0.0;
437        msg='Incorrect gauge depths returned'
438        assert num.allclose(elevation,-depth),msg
439        msg='incorrect gauge height time series returned'
440        assert num.allclose(stage,ha)
441        msg='incorrect gauge ua time series returned'
442        assert num.allclose(xvelocity,ua)
443        msg='incorrect gauge va time series returned'
444        assert num.allclose(yvelocity, -va) # South is positive in mux
445       
446
447       
448    def test_read_mux_platform_problem1(self):
449        """test_read_mux_platform_problem1
450       
451        This is to test a situation where read_mux returned
452        wrong values Win32
453
454        This test passes on Windows but test_read_mux_platform_problem2
455        does not
456        """
457       
458        from urs_ext import read_mux2
459       
460        verbose = False
461               
462        tide = 1.5
463        time_step_count = 10
464        time_step = 0.2
465        times_ref = num.arange(0, time_step_count*time_step, time_step)
466
467        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
468        n = len(lat_long_points)
469       
470        # Create different timeseries starting and ending at different times
471        first_tstep=num.ones(n, num.int)
472        first_tstep[0]+=2   # Point 0 starts at 2
473        first_tstep[1]+=4   # Point 1 starts at 4       
474        first_tstep[2]+=3   # Point 2 starts at 3
475       
476        last_tstep=(time_step_count)*num.ones(n,num.int)
477        last_tstep[0]-=1    # Point 0 ends 1 step early
478        last_tstep[1]-=2    # Point 1 ends 2 steps early               
479        last_tstep[4]-=3    # Point 4 ends 3 steps early       
480       
481        # Create varying elevation data (positive values for seafloor)
482        gauge_depth=20*num.ones(n,num.float)
483        for i in range(n):
484            gauge_depth[i] += i**2
485           
486        # Create data to be written to first mux file       
487        ha0=2*num.ones((n,time_step_count),num.float)
488        ha0[0]=num.arange(0,time_step_count)
489        ha0[1]=num.arange(time_step_count,2*time_step_count)
490        ha0[2]=num.arange(2*time_step_count,3*time_step_count)
491        ha0[3]=num.arange(3*time_step_count,4*time_step_count)
492        ua0=5*num.ones((n,time_step_count),num.float)
493        va0=-10*num.ones((n,time_step_count),num.float)
494
495        # Ensure data used to write mux file to be zero when gauges are
496        # not recording
497        for i in range(n):
498             # For each point
499             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
500                 # For timesteps before and after recording range
501                 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0                                 
502       
503        # Write first mux file to be combined by urs2sts
504        base_nameI, filesI = self.write_mux2(lat_long_points,
505                                             time_step_count, time_step,
506                                             first_tstep, last_tstep,
507                                             depth=gauge_depth,
508                                             ha=ha0,
509                                             ua=ua0,
510                                             va=va0)
511
512        # Create ordering file
513        permutation = ensure_numeric([4,0,2])
514
515        _, ordering_filename = tempfile.mkstemp('')
516        order_fid = open(ordering_filename, 'w') 
517        order_fid.write('index, longitude, latitude\n')
518        for index in permutation:
519            order_fid.write('%d, %f, %f\n' %(index, 
520                                             lat_long_points[index][1], 
521                                             lat_long_points[index][0]))
522        order_fid.close()
523       
524       
525
526        # -------------------------------------
527        # Now read files back and check values
528        weights = ensure_numeric([1.0])
529
530        # For each quantity read the associated list of source mux2 file with
531        # extention associated with that quantity
532        file_params=-1*num.ones(3,num.float) #[nsta,dt,nt]
533        OFFSET = 5
534
535        for j, file in enumerate(filesI):
536            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
537
538            number_of_selected_stations = data.shape[0]
539
540            # Index where data ends and parameters begin
541            parameters_index = data.shape[1]-OFFSET         
542         
543            for i in range(number_of_selected_stations):
544                if j == 0: assert num.allclose(data[i][:parameters_index], ha0[permutation[i], :])
545                if j == 1: assert num.allclose(data[i][:parameters_index], ua0[permutation[i], :])
546                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
547       
548        self.delete_mux(filesI)
549
550
551       
552       
553    def test_read_mux_platform_problem2(self):
554        """test_read_mux_platform_problem2
555       
556        This is to test a situation where read_mux returned
557        wrong values Win32
558
559        This test does not pass on Windows but test_read_mux_platform_problem1
560        does
561        """
562       
563        from urs_ext import read_mux2
564       
565        from anuga.config import single_precision as epsilon       
566       
567        verbose = False
568               
569        tide = 1.5
570        time_step_count = 10
571        time_step = 0.2
572       
573        times_ref = num.arange(0, time_step_count*time_step, time_step)
574       
575        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
576                           (-21.,115.), (-22., 117.)]
577        n = len(lat_long_points)
578       
579        # Create different timeseries starting and ending at different times
580        first_tstep=num.ones(n,num.int)
581        first_tstep[0]+=2   # Point 0 starts at 2
582        first_tstep[1]+=4   # Point 1 starts at 4       
583        first_tstep[2]+=3   # Point 2 starts at 3
584       
585        last_tstep=(time_step_count)*num.ones(n,num.int)
586        last_tstep[0]-=1    # Point 0 ends 1 step early
587        last_tstep[1]-=2    # Point 1 ends 2 steps early               
588        last_tstep[4]-=3    # Point 4 ends 3 steps early       
589       
590        # Create varying elevation data (positive values for seafloor)
591        gauge_depth=20*num.ones(n,num.float)
592        for i in range(n):
593            gauge_depth[i] += i**2
594           
595        # Create data to be written to second mux file       
596        ha1=num.ones((n,time_step_count),num.float)
597        ha1[0]=num.sin(times_ref)
598        ha1[1]=2*num.sin(times_ref - 3)
599        ha1[2]=5*num.sin(4*times_ref)
600        ha1[3]=num.sin(times_ref)
601        ha1[4]=num.sin(2*times_ref-0.7)
602               
603        ua1=num.zeros((n,time_step_count),num.float)
604        ua1[0]=3*num.cos(times_ref)       
605        ua1[1]=2*num.sin(times_ref-0.7)   
606        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
607        ua1[4]=2*num.ones(time_step_count)
608       
609        va1=num.zeros((n,time_step_count),num.float)
610        va1[0]=2*num.cos(times_ref-0.87)       
611        va1[1]=3*num.ones(time_step_count)
612        va1[3]=2*num.sin(times_ref-0.71)       
613       
614        # Ensure data used to write mux file to be zero when gauges are
615        # not recording
616        for i in range(n):
617             # For each point
618             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
619                 # For timesteps before and after recording range
620                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 
621
622
623        #print 'Second station to be written to MUX'
624        #print 'ha', ha1[0,:]
625        #print 'ua', ua1[0,:]
626        #print 'va', va1[0,:]
627       
628        # Write second mux file to be combined by urs2sts
629        base_nameII, filesII = self.write_mux2(lat_long_points,
630                                               time_step_count, time_step,
631                                               first_tstep, last_tstep,
632                                               depth=gauge_depth,
633                                               ha=ha1,
634                                               ua=ua1,
635                                               va=va1)
636
637
638
639
640        # Read mux file back and verify it's correcness
641
642        ####################################################
643        # FIXME (Ole): This is where the test should
644        # verify that the MUX files are correct.
645
646        #JJ: It appears as though
647        #that certain quantities are not being stored with enough precision
648        #inn muxfile or more likely that they are being cast into a
649        #lower precision when read in using read_mux2 Time step and q_time
650        # are equal but only to approx 1e-7
651        ####################################################
652
653        #define information as it should be stored in mus2 files
654        points_num=len(lat_long_points)
655        depth=gauge_depth
656        ha=ha1
657        ua=ua1
658        va=va1
659       
660        quantities = ['HA','UA','VA']
661        mux_names = [WAVEHEIGHT_MUX2_LABEL,
662                     EAST_VELOCITY_MUX2_LABEL,
663                     NORTH_VELOCITY_MUX2_LABEL]
664        quantities_init = [[],[],[]]
665        latlondeps = []
666        #irrelevant header information
667        ig=ilon=ilat=0
668        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
669        # urs binary is latitude fastest
670        for i,point in enumerate(lat_long_points):
671            lat = point[0]
672            lon = point[1]
673            _ , e, n = redfearn(lat, lon)
674            if depth is None:
675                this_depth = n
676            else:
677                this_depth = depth[i]
678            latlondeps.append([lat, lon, this_depth])
679
680            if ha is None:
681                this_ha = e
682                quantities_init[0].append(num.ones(time_step_count,num.float)*this_ha) # HA
683            else:
684                quantities_init[0].append(ha[i])
685            if ua is None:
686                this_ua = n
687                quantities_init[1].append(num.ones(time_step_count,num.float)*this_ua) # UA
688            else:
689                quantities_init[1].append(ua[i])
690            if va is None:
691                this_va = e
692                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
693            else:
694                quantities_init[2].append(va[i])
695
696        for i, q in enumerate(quantities):
697            #print
698            #print i, q
699           
700            q_time = num.zeros((time_step_count, points_num), num.float)
701            quantities_init[i] = ensure_numeric(quantities_init[i])
702            for time in range(time_step_count):
703                #print i, q, time, quantities_init[i][:,time]
704                q_time[time,:] = quantities_init[i][:,time]
705                #print i, q, time, q_time[time, :]
706
707           
708            filename = base_nameII + mux_names[i]
709            f = open(filename, 'rb')
710            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
711            #write mux 2 header
712            for latlondep in latlondeps:
713                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
714                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
715                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
716                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
717                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
718                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
719                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
720                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
721                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
722                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
723                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
724                assert abs(az-unpack('f',f.read(4))[0])<epsilon
725                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
726               
727                x = unpack('f', f.read(4))[0]
728                #print time_step
729                #print x
730                assert abs(time_step-x)<epsilon
731                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
732                for j in range(4): # identifier
733                    assert abs(id-unpack('i',f.read(4))[0])<epsilon
734
735            #first_tstep=1
736            #last_tstep=time_step_count
737            for i,latlondep in enumerate(latlondeps):
738                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
739            for i,latlondep in enumerate(latlondeps):
740                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
741
742            # Find when first station starts recording
743            min_tstep = min(first_tstep)
744            # Find when all stations have stopped recording
745            max_tstep = max(last_tstep)
746
747            #for time in  range(time_step_count):
748            for time in range(min_tstep-1,max_tstep):
749                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
750                for point_i in range(points_num):
751                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
752                        x = unpack('f',f.read(4))[0]
753                        #print time, x, q_time[time, point_i]
754                        if q == 'VA': x = -x # South is positive in MUX
755                        assert abs(q_time[time, point_i]-x)<epsilon
756
757            f.close()
758                                               
759        # Create ordering file
760        permutation = ensure_numeric([4,0,2])
761
762       #  _, ordering_filename = tempfile.mkstemp('')
763#         order_fid = open(ordering_filename, 'w') 
764#         order_fid.write('index, longitude, latitude\n')
765#         for index in permutation:
766#             order_fid.write('%d, %f, %f\n' %(index,
767#                                              lat_long_points[index][1],
768#                                              lat_long_points[index][0]))
769#         order_fid.close()
770       
771        # -------------------------------------
772        # Now read files back and check values
773        weights = ensure_numeric([1.0])
774
775        # For each quantity read the associated list of source mux2 file with
776        # extention associated with that quantity
777        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
778        OFFSET = 5
779
780        for j, file in enumerate(filesII):
781            # Read stage, u, v enumerated as j
782            #print 'Reading', j, file
783            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
784
785            #print 'Data received by Python'
786            #print data[1][8]
787            number_of_selected_stations = data.shape[0]
788
789            # Index where data ends and parameters begin
790            parameters_index = data.shape[1]-OFFSET         
791                 
792            quantity=num.zeros((number_of_selected_stations, parameters_index), num.float)
793           
794           
795            for i in range(number_of_selected_stations):
796       
797                #print i, parameters_index
798                #print quantity[i][:]
799                if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :])
800                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
801                if j == 2:
802                    # FIXME (Ole): This is where the output is wrong on Win32
803                   
804                    #print
805                    #print j, i
806                    #print 'Input'
807                    #print 'u', ua1[permutation[i], 8]       
808                    #print 'v', va1[permutation[i], 8]
809               
810                    #print 'Output'
811                    #print 'v ', data[i][:parameters_index][8] 
812
813                    # South is positive in MUX
814                    #print "data[i][:parameters_index]", data[i][:parameters_index]
815                    #print "-va1[permutation[i], :]", -va1[permutation[i], :]
816                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
817       
818        self.delete_mux(filesII)
819           
820    def test_read_mux_platform_problem3(self):
821       
822        # This is to test a situation where read_mux returned
823        # wrong values Win32
824
825       
826        from urs_ext import read_mux2
827       
828        from anuga.config import single_precision as epsilon       
829       
830        verbose = False
831               
832        tide = 1.5
833        time_step_count = 10
834        time_step = 0.02
835
836        '''
837        Win results
838        time_step = 0.2000001
839        This is OK       
840        '''
841       
842        '''
843        Win results
844        time_step = 0.20000001
845
846        ======================================================================
847ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager)
848----------------------------------------------------------------------
849Traceback (most recent call last):
850  File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3
851    ha1[0]=num.sin(times_ref)
852ValueError: matrices are not aligned for copy
853
854        '''
855
856        '''
857        Win results
858        time_step = 0.200000001
859        FAIL
860         assert num.allclose(data[i][:parameters_index],
861         -va1[permutation[i], :])
862        '''
863        times_ref = num.arange(0, time_step_count*time_step, time_step)
864        #print "times_ref", times_ref
865       
866        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
867                           (-21.,115.), (-22., 117.)]
868        stations = len(lat_long_points)
869       
870        # Create different timeseries starting and ending at different times
871        first_tstep=num.ones(stations, num.int)
872        first_tstep[0]+=2   # Point 0 starts at 2
873        first_tstep[1]+=4   # Point 1 starts at 4       
874        first_tstep[2]+=3   # Point 2 starts at 3
875       
876        last_tstep=(time_step_count)*num.ones(stations, num.int)
877        last_tstep[0]-=1    # Point 0 ends 1 step early
878        last_tstep[1]-=2    # Point 1 ends 2 steps early               
879        last_tstep[4]-=3    # Point 4 ends 3 steps early       
880       
881        # Create varying elevation data (positive values for seafloor)
882        gauge_depth=20*num.ones(stations, num.float)
883        for i in range(stations):
884            gauge_depth[i] += i**2
885           
886        # Create data to be written to second mux file       
887        ha1=num.ones((stations,time_step_count), num.float)
888        ha1[0]=num.sin(times_ref)
889        ha1[1]=2*num.sin(times_ref - 3)
890        ha1[2]=5*num.sin(4*times_ref)
891        ha1[3]=num.sin(times_ref)
892        ha1[4]=num.sin(2*times_ref-0.7)
893               
894        ua1=num.zeros((stations,time_step_count),num.float)
895        ua1[0]=3*num.cos(times_ref)       
896        ua1[1]=2*num.sin(times_ref-0.7)   
897        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
898        ua1[4]=2*num.ones(time_step_count)
899       
900        va1=num.zeros((stations,time_step_count),num.float)
901        va1[0]=2*num.cos(times_ref-0.87)       
902        va1[1]=3*num.ones(time_step_count)
903        va1[3]=2*num.sin(times_ref-0.71)       
904        #print "va1[0]", va1[0]  # The 8th element is what will go bad.
905        # Ensure data used to write mux file to be zero when gauges are
906        # not recording
907        for i in range(stations):
908             # For each point
909             for j in range(0, first_tstep[i]-1) + range(last_tstep[i],
910                                                         time_step_count):
911                 # For timesteps before and after recording range
912                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 
913
914
915        #print 'Second station to be written to MUX'
916        #print 'ha', ha1[0,:]
917        #print 'ua', ua1[0,:]
918        #print 'va', va1[0,:]
919       
920        # Write second mux file to be combined by urs2sts
921        base_nameII, filesII = self.write_mux2(lat_long_points,
922                                               time_step_count, time_step,
923                                               first_tstep, last_tstep,
924                                               depth=gauge_depth,
925                                               ha=ha1,
926                                               ua=ua1,
927                                               va=va1)
928        #print "filesII", filesII
929
930
931
932
933        # Read mux file back and verify it's correcness
934
935        ####################################################
936        # FIXME (Ole): This is where the test should
937        # verify that the MUX files are correct.
938
939        #JJ: It appears as though
940        #that certain quantities are not being stored with enough precision
941        #inn muxfile or more likely that they are being cast into a
942        #lower precision when read in using read_mux2 Time step and q_time
943        # are equal but only to approx 1e-7
944        ####################################################
945
946        #define information as it should be stored in mus2 files
947        points_num=len(lat_long_points)
948        depth=gauge_depth
949        ha=ha1
950        ua=ua1
951        va=va1
952       
953        quantities = ['HA','UA','VA']
954        mux_names = [WAVEHEIGHT_MUX2_LABEL,
955                     EAST_VELOCITY_MUX2_LABEL,
956                     NORTH_VELOCITY_MUX2_LABEL]
957        quantities_init = [[],[],[]]
958        latlondeps = []
959        #irrelevant header information
960        ig=ilon=ilat=0
961        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
962        # urs binary is latitude fastest
963        for i,point in enumerate(lat_long_points):
964            lat = point[0]
965            lon = point[1]
966            _ , e, n = redfearn(lat, lon)
967            if depth is None:
968                this_depth = n
969            else:
970                this_depth = depth[i]
971            latlondeps.append([lat, lon, this_depth])
972
973            if ha is None:
974                this_ha = e
975                quantities_init[0].append(num.ones(time_step_count,
976                                                   num.float)*this_ha) # HA
977            else:
978                quantities_init[0].append(ha[i])
979            if ua is None:
980                this_ua = n
981                quantities_init[1].append(num.ones(time_step_count,
982                                                   num.float)*this_ua) # UA
983            else:
984                quantities_init[1].append(ua[i])
985            if va is None:
986                this_va = e
987                quantities_init[2].append(num.ones(time_step_count,
988                                                   num.float)*this_va) #
989            else:
990                quantities_init[2].append(va[i])
991
992        for i, q in enumerate(quantities):
993            #print
994            #print i, q
995           
996            q_time = num.zeros((time_step_count, points_num), num.float)
997            quantities_init[i] = ensure_numeric(quantities_init[i])
998            for time in range(time_step_count):
999                #print i, q, time, quantities_init[i][:,time]
1000                q_time[time,:] = quantities_init[i][:,time]
1001                #print i, q, time, q_time[time, :]
1002
1003           
1004            filename = base_nameII + mux_names[i]
1005            f = open(filename, 'rb')
1006            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
1007            #write mux 2 header
1008            for latlondep in latlondeps:
1009                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
1010                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
1011                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
1012                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
1013                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
1014                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
1015                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
1016                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
1017                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
1018                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
1019                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
1020                assert abs(az-unpack('f',f.read(4))[0])<epsilon
1021                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
1022               
1023                x = unpack('f', f.read(4))[0]
1024                #print time_step
1025                #print x
1026                assert abs(time_step-x)<epsilon
1027                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
1028                for j in range(4): # identifier
1029                    assert abs(id-unpack('i',f.read(4))[0])<epsilon
1030
1031            #first_tstep=1
1032            #last_tstep=time_step_count
1033            for i,latlondep in enumerate(latlondeps):
1034                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
1035            for i,latlondep in enumerate(latlondeps):
1036                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
1037
1038            # Find when first station starts recording
1039            min_tstep = min(first_tstep)
1040            # Find when all stations have stopped recording
1041            max_tstep = max(last_tstep)
1042
1043            #for time in  range(time_step_count):
1044            for time in range(min_tstep-1,max_tstep):
1045                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
1046                for point_i in range(points_num):
1047                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
1048                        x = unpack('f',f.read(4))[0]
1049                        #print time, x, q_time[time, point_i]
1050                        if q == 'VA': x = -x # South is positive in MUX
1051                        #print q+" q_time[%d, %d] = %f" %(time, point_i,
1052                                                      #q_time[time, point_i])
1053                        assert abs(q_time[time, point_i]-x)<epsilon
1054
1055            f.close()
1056                           
1057        permutation = ensure_numeric([4,0,2])
1058                   
1059        # Create ordering file
1060#         _, ordering_filename = tempfile.mkstemp('')
1061#         order_fid = open(ordering_filename, 'w') 
1062#         order_fid.write('index, longitude, latitude\n')
1063#         for index in permutation:
1064#             order_fid.write('%d, %f, %f\n' %(index,
1065#                                              lat_long_points[index][1],
1066#                                              lat_long_points[index][0]))
1067#         order_fid.close()
1068       
1069        # -------------------------------------
1070        # Now read files back and check values
1071        weights = ensure_numeric([1.0])
1072
1073        # For each quantity read the associated list of source mux2 file with
1074        # extention associated with that quantity
1075        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
1076        OFFSET = 5
1077
1078        for j, file in enumerate(filesII):
1079            # Read stage, u, v enumerated as j
1080            #print 'Reading', j, file
1081            #print "file", file
1082            data = read_mux2(1, [file], weights, file_params,
1083                             permutation, verbose)
1084            #print str(j) + "data", data
1085
1086            #print 'Data received by Python'
1087            #print data[1][8]
1088            number_of_selected_stations = data.shape[0]
1089            #print "number_of_selected_stations", number_of_selected_stations
1090            #print "stations", stations
1091
1092            # Index where data ends and parameters begin
1093            parameters_index = data.shape[1]-OFFSET         
1094                 
1095            for i in range(number_of_selected_stations):
1096       
1097                #print i, parameters_index
1098                if j == 0:
1099                    assert num.allclose(data[i][:parameters_index],
1100                                        ha1[permutation[i], :])
1101                   
1102                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
1103                if j == 2:
1104                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
1105       
1106        self.delete_mux(filesII)     
1107
1108
1109         
1110    def test_urs2sts_nonstandard_projection_reverse(self):
1111        """
1112        Test that a point not in the specified zone can occur first
1113        """
1114        tide=0
1115        time_step_count = 3
1116        time_step = 2
1117        lat_long_points =[(-21.,113.5),(-21.,114.5),(-21.,114.), (-21.,115.)]
1118        n=len(lat_long_points)
1119        first_tstep=num.ones(n,num.int)
1120        first_tstep[0]+=1
1121        first_tstep[2]+=1
1122        last_tstep=(time_step_count)*num.ones(n,num.int)
1123        last_tstep[0]-=1
1124
1125        gauge_depth=20*num.ones(n,num.float)
1126        ha=2*num.ones((n,time_step_count),num.float)
1127        ha[0]=num.arange(0,time_step_count)
1128        ha[1]=num.arange(time_step_count,2*time_step_count)
1129        ha[2]=num.arange(2*time_step_count,3*time_step_count)
1130        ha[3]=num.arange(3*time_step_count,4*time_step_count)
1131        ua=5*num.ones((n,time_step_count),num.float)
1132        va=-10*num.ones((n,time_step_count),num.float)
1133
1134        base_name, files = self.write_mux2(lat_long_points,
1135                                      time_step_count, time_step,
1136                                      first_tstep, last_tstep,
1137                                      depth=gauge_depth,
1138                                      ha=ha,
1139                                      ua=ua,
1140                                      va=va)
1141
1142        urs2sts(base_name,
1143                basename_out=base_name, 
1144                zone=50,
1145                mean_stage=tide,verbose=False)
1146
1147        # now I want to check the sts file ...
1148        sts_file = base_name + '.sts'
1149
1150        #Let's interigate the sww file
1151        # Note, the sww info is not gridded.  It is point data.
1152        fid = NetCDFFile(sts_file)
1153
1154        # Make x and y absolute
1155        x = fid.variables['x'][:]
1156        y = fid.variables['y'][:]
1157
1158        geo_reference = Geo_reference(NetCDFObject=fid)
1159        points = geo_reference.get_absolute(map(None, x, y))
1160        points = ensure_numeric(points)
1161
1162        x = points[:,0]
1163        y = points[:,1]
1164
1165        # Check that all coordinate are correctly represented       
1166        # Using the non standard projection (50)
1167        for i in range(4):
1168            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1],
1169                                  zone=50) 
1170            assert num.allclose([x[i],y[i]], [e,n])
1171            assert zone==geo_reference.zone
1172       
1173        self.delete_mux(files)
1174
1175           
1176    def test_urs2stsII(self):
1177        """
1178        Test multiple sources
1179        """
1180        tide=0
1181        time_step_count = 3
1182        time_step = 2
1183        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
1184        n=len(lat_long_points)
1185        first_tstep=num.ones(n,num.int)
1186        first_tstep[0]+=1
1187        first_tstep[2]+=1
1188        last_tstep=(time_step_count)*num.ones(n,num.int)
1189        last_tstep[0]-=1
1190
1191        gauge_depth=20*num.ones(n,num.float)
1192        ha=2*num.ones((n,time_step_count),num.float)
1193        ha[0]=num.arange(0,time_step_count)
1194        ha[1]=num.arange(time_step_count,2*time_step_count)
1195        ha[2]=num.arange(2*time_step_count,3*time_step_count)
1196        ha[3]=num.arange(3*time_step_count,4*time_step_count)
1197        ua=5*num.ones((n,time_step_count),num.float)
1198        va=-10*num.ones((n,time_step_count),num.float)
1199
1200        # Create two identical mux files to be combined by urs2sts
1201        base_nameI, filesI = self.write_mux2(lat_long_points,
1202                                             time_step_count, time_step,
1203                                             first_tstep, last_tstep,
1204                                             depth=gauge_depth,
1205                                             ha=ha,
1206                                             ua=ua,
1207                                             va=va)
1208
1209        base_nameII, filesII = self.write_mux2(lat_long_points,
1210                                               time_step_count, time_step,
1211                                               first_tstep, last_tstep,
1212                                               depth=gauge_depth,
1213                                               ha=ha,
1214                                               ua=ua,
1215                                               va=va)
1216
1217        # Call urs2sts with multiple mux files
1218        urs2sts([base_nameI, base_nameII], 
1219                basename_out=base_nameI, 
1220                weights=[1.0, 1.0],
1221                mean_stage=tide,
1222                verbose=False)
1223
1224        # now I want to check the sts file ...
1225        sts_file = base_nameI + '.sts'
1226
1227        #Let's interrogate the sts file
1228        # Note, the sts info is not gridded.  It is point data.
1229        fid = NetCDFFile(sts_file)
1230
1231        # Make x and y absolute
1232        x = fid.variables['x'][:]
1233        y = fid.variables['y'][:]
1234
1235        geo_reference = Geo_reference(NetCDFObject=fid)
1236        points = geo_reference.get_absolute(map(None, x, y))
1237        points = ensure_numeric(points)
1238
1239        x = points[:,0]
1240        y = points[:,1]
1241
1242        #Check that first coordinate is correctly represented       
1243        #Work out the UTM coordinates for first point
1244        zone, e, n = redfearn(lat_long_points[0][0], lat_long_points[0][1]) 
1245        assert num.allclose([x[0],y[0]], [e,n])
1246
1247        #Check the time vector
1248        times = fid.variables['time'][:]
1249
1250        times_actual = []
1251        for i in range(time_step_count):
1252            times_actual.append(time_step * i)
1253
1254        assert num.allclose(ensure_numeric(times),
1255                            ensure_numeric(times_actual))
1256
1257        #Check first value
1258        stage = fid.variables['stage'][:]
1259        xmomentum = fid.variables['xmomentum'][:]
1260        ymomentum = fid.variables['ymomentum'][:]
1261        elevation = fid.variables['elevation'][:]
1262
1263        # Set original data used to write mux file to be zero when gauges are
1264        # not recdoring
1265       
1266        ha[0][0]=0.0
1267        ha[0][time_step_count-1]=0.0
1268        ha[2][0]=0.0
1269        ua[0][0]=0.0
1270        ua[0][time_step_count-1]=0.0
1271        ua[2][0]=0.0
1272        va[0][0]=0.0
1273        va[0][time_step_count-1]=0.0
1274        va[2][0]=0.0;
1275
1276        # The stage stored in the .sts file should be the sum of the stage
1277        # in the two mux2 files because both have weights = 1. In this case
1278        # the mux2 files are the same so stage == 2.0 * ha
1279        #print 2.0*num.transpose(ha) - stage
1280        assert num.allclose(2.0*num.transpose(ha), stage)  #Meters
1281
1282        #Check the momentums - ua
1283        #momentum = velocity*(stage-elevation)
1284        # elevation = - depth
1285        #momentum = velocity_ua *(stage+depth)
1286
1287        depth=num.zeros((len(lat_long_points),time_step_count),num.float)
1288        for i in range(len(lat_long_points)):
1289            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
1290            #2.0*ha necessary because using two files with weights=1 are used
1291
1292        # The xmomentum stored in the .sts file should be the sum of the ua
1293        # in the two mux2 files multiplied by the depth.
1294        assert num.allclose(2.0*num.transpose(ua*depth), xmomentum) 
1295
1296        #Check the momentums - va
1297        #momentum = velocity*(stage-elevation)
1298        # elevation = - depth
1299        #momentum = velocity_va *(stage+depth)
1300
1301        # The ymomentum stored in the .sts file should be the sum of the va
1302        # in the two mux2 files multiplied by the depth.
1303        assert num.allclose(2.0*num.transpose(va*depth), ymomentum)
1304
1305        # check the elevation values.
1306        # -ve since urs measures depth, sww meshers height,
1307        assert num.allclose(-elevation, gauge_depth)  #Meters
1308
1309        fid.close()
1310        self.delete_mux(filesI)
1311        self.delete_mux(filesII)
1312        os.remove(sts_file)       
1313
1314
1315    def test_urs2sts0(self):
1316        """
1317        Test single source
1318        """
1319        tide=0
1320        time_step_count = 3
1321        time_step = 2
1322        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
1323        n=len(lat_long_points)
1324        first_tstep=num.ones(n,num.int)
1325        first_tstep[0]+=1
1326        first_tstep[2]+=1
1327        last_tstep=(time_step_count)*num.ones(n,num.int)
1328        last_tstep[0]-=1
1329
1330        gauge_depth=20*num.ones(n,num.float)
1331        ha=2*num.ones((n,time_step_count),num.float)
1332        ha[0]=num.arange(0,time_step_count)
1333        ha[1]=num.arange(time_step_count,2*time_step_count)
1334        ha[2]=num.arange(2*time_step_count,3*time_step_count)
1335        ha[3]=num.arange(3*time_step_count,4*time_step_count)
1336        ua=5*num.ones((n,time_step_count),num.float)
1337        va=-10*num.ones((n,time_step_count),num.float)
1338
1339        base_name, files = self.write_mux2(lat_long_points,
1340                                      time_step_count, time_step,
1341                                      first_tstep, last_tstep,
1342                                      depth=gauge_depth,
1343                                      ha=ha,
1344                                      ua=ua,
1345                                      va=va)
1346
1347        urs2sts(base_name,
1348                basename_out=base_name, 
1349                mean_stage=tide,verbose=False)
1350
1351        # now I want to check the sts file ...
1352        sts_file = base_name + '.sts'
1353
1354        #Let's interigate the sww file
1355        # Note, the sww info is not gridded.  It is point data.
1356        fid = NetCDFFile(sts_file)
1357
1358        # Make x and y absolute
1359        x = fid.variables['x'][:]
1360        y = fid.variables['y'][:]
1361
1362        geo_reference = Geo_reference(NetCDFObject=fid)
1363        points = geo_reference.get_absolute(map(None, x, y))
1364        points = ensure_numeric(points)
1365
1366        x = points[:,0]
1367        y = points[:,1]
1368
1369        #Check that first coordinate is correctly represented       
1370        #Work out the UTM coordinates for first point
1371        for i in range(4):
1372            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1]) 
1373            assert num.allclose([x[i],y[i]], [e,n])
1374
1375        #Check the time vector
1376        times = fid.variables['time'][:]
1377
1378        times_actual = []
1379        for i in range(time_step_count):
1380            times_actual.append(time_step * i)
1381
1382        assert num.allclose(ensure_numeric(times),
1383                            ensure_numeric(times_actual))
1384
1385        #Check first value
1386        stage = fid.variables['stage'][:]
1387        xmomentum = fid.variables['xmomentum'][:]
1388        ymomentum = fid.variables['ymomentum'][:]
1389        elevation = fid.variables['elevation'][:]
1390
1391        # Set original data used to write mux file to be zero when gauges are
1392        #not recdoring
1393        ha[0][0]=0.0
1394        ha[0][time_step_count-1]=0.0;
1395        ha[2][0]=0.0;
1396        ua[0][0]=0.0
1397        ua[0][time_step_count-1]=0.0;
1398        ua[2][0]=0.0;
1399        va[0][0]=0.0
1400        va[0][time_step_count-1]=0.0;
1401        va[2][0]=0.0;
1402
1403        assert num.allclose(num.transpose(ha),stage)  #Meters
1404
1405        #Check the momentums - ua
1406        #momentum = velocity*(stage-elevation)
1407        # elevation = - depth
1408        #momentum = velocity_ua *(stage+depth)
1409
1410        depth=num.zeros((len(lat_long_points),time_step_count),num.float)
1411        for i in range(len(lat_long_points)):
1412            depth[i]=gauge_depth[i]+tide+ha[i]
1413        assert num.allclose(num.transpose(ua*depth),xmomentum) 
1414
1415        #Check the momentums - va
1416        #momentum = velocity*(stage-elevation)
1417        # elevation = - depth
1418        #momentum = velocity_va *(stage+depth)
1419
1420        assert num.allclose(num.transpose(va*depth),ymomentum)
1421
1422        # check the elevation values.
1423        # -ve since urs measures depth, sww meshers height,
1424        assert num.allclose(-elevation, gauge_depth)  #Meters
1425
1426        fid.close()
1427        self.delete_mux(files)
1428        os.remove(sts_file)
1429
1430    def test_urs2sts_nonstandard_meridian(self):
1431        """
1432        Test single source using the meridian from zone 50 as a nonstandard meridian
1433        """
1434        tide=0
1435        time_step_count = 3
1436        time_step = 2
1437        lat_long_points =[(-21.,114.5),(-21.,113.5),(-21.,114.), (-21.,115.)]
1438        n=len(lat_long_points)
1439        first_tstep=num.ones(n,num.int)
1440        first_tstep[0]+=1
1441        first_tstep[2]+=1
1442        last_tstep=(time_step_count)*num.ones(n,num.int)
1443        last_tstep[0]-=1
1444
1445        gauge_depth=20*num.ones(n,num.float)
1446        ha=2*num.ones((n,time_step_count),num.float)
1447        ha[0]=num.arange(0,time_step_count)
1448        ha[1]=num.arange(time_step_count,2*time_step_count)
1449        ha[2]=num.arange(2*time_step_count,3*time_step_count)
1450        ha[3]=num.arange(3*time_step_count,4*time_step_count)
1451        ua=5*num.ones((n,time_step_count),num.float)
1452        va=-10*num.ones((n,time_step_count),num.float)
1453
1454        base_name, files = self.write_mux2(lat_long_points,
1455                                           time_step_count, time_step,
1456                                           first_tstep, last_tstep,
1457                                           depth=gauge_depth,
1458                                           ha=ha,
1459                                           ua=ua,
1460                                           va=va)
1461
1462        urs2sts(base_name,
1463                basename_out=base_name, 
1464                central_meridian=123,
1465                mean_stage=tide,
1466                verbose=False)
1467
1468        # now I want to check the sts file ...
1469        sts_file = base_name + '.sts'
1470
1471        #Let's interigate the sww file
1472        # Note, the sww info is not gridded.  It is point data.
1473        fid = NetCDFFile(sts_file)
1474
1475        # Make x and y absolute
1476        x = fid.variables['x'][:]
1477        y = fid.variables['y'][:]
1478
1479        geo_reference = Geo_reference(NetCDFObject=fid)
1480        points = geo_reference.get_absolute(map(None, x, y))
1481        points = ensure_numeric(points)
1482
1483        x = points[:,0]
1484        y = points[:,1]
1485
1486        # Check that all coordinate are correctly represented       
1487        # Using the non standard projection (50)
1488        for i in range(4):
1489            zone, e, n = redfearn(lat_long_points[i][0],
1490                                  lat_long_points[i][1],
1491                                  central_meridian=123)
1492            assert num.allclose([x[i],y[i]], [e,n])
1493            assert zone==-1
1494       
1495        self.delete_mux(files)
1496 
1497    def test_Urs_points(self):
1498        time_step_count = 3
1499        time_step = 2
1500        lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)]
1501        base_name, files = self.write_mux(lat_long_points,
1502                                          time_step_count, time_step)
1503        for file in files:
1504            urs = Read_urs(file)
1505            assert time_step_count == urs.time_step_count
1506            assert time_step == urs.time_step
1507
1508            for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep):
1509                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
1510                    assert num.allclose(n, dep[2])
1511                       
1512            count = 0
1513            for slice in urs:
1514                count += 1
1515                #print slice
1516                for lat_lon, quantity in map(None, lat_long_points, slice):
1517                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
1518                    #print "quantity", quantity
1519                    #print "e", e
1520                    #print "n", n
1521                    if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \
1522                           file[-5:] == NORTH_VELOCITY_LABEL[-5:] :
1523                        assert num.allclose(e, quantity)
1524                    if file[-5:] == EAST_VELOCITY_LABEL[-5:]:
1525                        assert num.allclose(n, quantity)
1526            assert count == time_step_count
1527                     
1528        self.delete_mux(files)       
1529
1530
1531
1532       
1533################################################################################
1534
1535if __name__ == "__main__":
1536    suite = unittest.makeSuite(Test_Mux,'test')
1537    runner = unittest.TextTestRunner() #verbosity=2)
1538    runner.run(suite)
1539       
Note: See TracBrowser for help on using the repository browser.