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

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

Moved mux file module into generic file module folder.

File size: 45.2 KB
Line 
1import unittest
2import tempfile
3import numpy as num
4import os
5from struct import pack, unpack
6
7from anuga.utilities.numerical_tools import ensure_numeric
8from anuga.coordinate_transforms.redfearn import redfearn
9
10from mux import WAVEHEIGHT_MUX2_LABEL, EAST_VELOCITY_MUX2_LABEL, \
11                NORTH_VELOCITY_MUX2_LABEL
12               
13from mux import read_mux2_py
14
15class TestCase(unittest.TestCase):
16    def setUp(self):
17        pass
18
19    def tearDown(self):
20        pass
21
22    def write_mux(self, lat_long_points, time_step_count, time_step,
23                  depth=None, ha=None, ua=None, va=None):
24        """
25        This will write 3 non-gridded mux files, for testing.
26        If no quantities are passed in,
27        na and va quantities will be the Easting values.
28        Depth and ua will be the Northing value.
29       
30        The mux file format has south as positive so
31        this function will swap the sign for va. 
32        """
33
34        #print "lat_long_points", lat_long_points
35        #print "time_step_count",time_step_count
36        #print "time_step",
37
38       
39        points_num = len(lat_long_points)
40        lonlatdeps = []
41        quantities = ['HA','UA','VA']
42       
43        mux_names = [WAVEHEIGHT_MUX_LABEL,
44                     EAST_VELOCITY_LABEL,
45                     NORTH_VELOCITY_LABEL]
46        quantities_init = [[],[],[]]
47        # urs binary is latitude fastest
48        for point in lat_long_points:
49            lat = point[0]
50            lon = point[1]
51            _ , e, n = redfearn(lat, lon)
52            if depth is None:
53                this_depth = n
54            else:
55                this_depth = depth
56            if ha is None:
57                this_ha = e
58            else:
59                this_ha = ha
60            if ua is None:
61                this_ua = n
62            else:
63                this_ua = ua
64            if va is None:
65                this_va = e   
66            else:
67                this_va = va         
68            lonlatdeps.append([lon, lat, this_depth])
69            quantities_init[0].append(this_ha) # HA
70            quantities_init[1].append(this_ua) # UA
71            quantities_init[2].append(this_va) # VA
72               
73        file_handle, base_name = tempfile.mkstemp("")
74        os.close(file_handle)
75        os.remove(base_name)
76
77        files = []       
78        for i, q in enumerate(quantities): 
79            quantities_init[i] = ensure_numeric(quantities_init[i])
80            #print "HA_init", HA_init
81            q_time = num.zeros((time_step_count, points_num), num.float)
82            for time in range(time_step_count):
83                q_time[time,:] = quantities_init[i] #* time * 4
84           
85            #Write C files
86            columns = 3 # long, lat , depth
87            file = base_name + mux_names[i]
88            #print "base_name file",file
89            f = open(file, 'wb')
90            files.append(file)
91            f.write(pack('i',points_num))
92            f.write(pack('i',time_step_count))
93            f.write(pack('f',time_step))
94
95            #write lat/long info
96            for lonlatdep in lonlatdeps:
97                for float in lonlatdep:
98                    f.write(pack('f',float))
99                   
100            # Write quantity info
101            for time in  range(time_step_count):
102                for point_i in range(points_num):
103                    f.write(pack('f',q_time[time,point_i]))
104                    #print " mux_names[i]", mux_names[i]
105                    #print "f.write(pack('f',q_time[time,i]))", q_time[time,point_i]
106            f.close()
107        return base_name, files
108       
109   
110    def delete_mux(self, files):
111        for file in files:
112            os.remove(file)
113       
114    def write_mux2(self, lat_long_points, time_step_count, time_step,
115                   first_tstep, last_tstep,
116                   depth=None, ha=None, ua=None, va=None):
117        """
118        This will write 3 non-gridded mux files, for testing.
119        If no quantities are passed in,
120        na and va quantities will be the Easting values.
121        Depth and ua will be the Northing value.
122        """
123        #print "lat_long_points", lat_long_points
124        #print "time_step_count",time_step_count
125        #print "time_step",
126
127        #irrelevant header information
128        ig=ilon=ilat=0
129        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
130
131        points_num = len(lat_long_points)
132        latlondeps = []
133        quantities = ['HA','UA','VA']
134
135        mux_names = [WAVEHEIGHT_MUX2_LABEL,
136                     EAST_VELOCITY_MUX2_LABEL,
137                     NORTH_VELOCITY_MUX2_LABEL]
138
139        msg='first_tstep and last_step arrays must have same length as number of points'
140        assert len(first_tstep)==points_num,msg
141        assert len(last_tstep)==points_num,msg
142
143        if depth is not None:
144            depth=ensure_numeric(depth)
145            assert len(depth)==points_num
146        if ha is not None:
147            ha=ensure_numeric(ha)
148            assert ha.shape==(points_num,time_step_count)
149        if ua is not None:
150            ua=ensure_numeric(ua)
151            assert ua.shape==(points_num,time_step_count)
152        if va is not None:
153            va=ensure_numeric(va)
154            assert va.shape==(points_num,time_step_count)
155
156        quantities_init = [[],[],[]]
157        # urs binary is latitude fastest
158        for i,point in enumerate(lat_long_points):
159            lat = point[0]
160            lon = point[1]
161            _ , e, n = redfearn(lat, lon)
162            if depth is None:
163                this_depth = n
164            else:
165                this_depth = depth[i]
166            latlondeps.append([lat, lon, this_depth])
167
168            if ha is None:
169                this_ha = e
170                quantities_init[0].append(num.ones(time_step_count,num.float)*this_ha) # HA
171            else:
172                quantities_init[0].append(ha[i])
173            if ua is None:
174                this_ua = n
175                quantities_init[1].append(num.ones(time_step_count,num.float)*this_ua) # UA
176            else:
177                quantities_init[1].append(ua[i])
178            if va is None:
179                this_va = e
180                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
181            else:
182                quantities_init[2].append(-va[i]) # South is negative in MUX
183
184        file_handle, base_name = tempfile.mkstemp("write_mux2")
185        os.close(file_handle)
186        os.remove(base_name)
187
188        files = []       
189        for i, q in enumerate(quantities):
190            q_time = num.zeros((time_step_count, points_num), num.float)
191            quantities_init[i] = ensure_numeric(quantities_init[i])
192            for time in range(time_step_count):
193                #print i, q, time, quantities_init[i][:,time]
194                q_time[time,:] = quantities_init[i][:,time]
195                #print i, q, time, q_time[time, :]               
196
197            #Write C files
198            columns = 3 # long, lat , depth
199            file = base_name + mux_names[i]
200           
201            #print 'base_name file', file
202            f = open(file, 'wb')
203            files.append(file)
204
205            f.write(pack('i',points_num))
206            #write mux 2 header
207            for latlondep in latlondeps:
208                f.write(pack('f',latlondep[0]))
209                f.write(pack('f',latlondep[1]))
210                f.write(pack('f',mcolat))
211                f.write(pack('f',mcolon))
212                f.write(pack('i',ig))
213                f.write(pack('i',ilon))
214                f.write(pack('i',ilat))
215                f.write(pack('f',latlondep[2]))
216                f.write(pack('f',centerlat))
217                f.write(pack('f',centerlon))
218                f.write(pack('f',offset))
219                f.write(pack('f',az))
220                f.write(pack('f',baz))
221                f.write(pack('f',time_step))
222                f.write(pack('i',time_step_count))
223                for j in range(4): # identifier
224                    f.write(pack('f',id))   
225
226            #first_tstep=1
227            #last_tstep=time_step_count
228            for i,latlondep in enumerate(latlondeps):
229                f.write(pack('i',first_tstep[i]))
230            for i,latlondep in enumerate(latlondeps):
231                f.write(pack('i',last_tstep[i]))
232
233            # Find when first station starts recording
234            min_tstep = min(first_tstep)
235            # Find when all stations have stopped recording
236            max_tstep = max(last_tstep)
237
238            #for time in  range(time_step_count):
239            for time in range(min_tstep-1,max_tstep):
240                f.write(pack('f',time*time_step))               
241                for point_i in range(points_num):
242                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
243                        #print 'writing', time, point_i, q_time[time, point_i]
244                        f.write(pack('f', q_time[time, point_i]))
245            f.close()
246
247        return base_name, files
248
249    def test_urs2sts_read_mux2_pyI(self):
250        """test_urs2sts_read_mux2_pyI(self):
251        Constant stage,momentum at each gauge
252        """
253        tide = 1
254        time_step_count = 3
255        time_step = 2
256        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
257        n=len(lat_long_points)
258        first_tstep=num.ones(n,num.int)
259        last_tstep=time_step_count*num.ones(n,num.int)
260        depth=20*num.ones(n,num.float)
261        ha=2*num.ones((n,time_step_count),num.float)
262        ua=5*num.ones((n,time_step_count),num.float)
263        va=-10*num.ones((n,time_step_count),num.float)
264        #-ve added to take into account mux file format where south is positive.
265        base_name, files = self.write_mux2(lat_long_points,
266                                      time_step_count, time_step,
267                                      first_tstep, last_tstep,
268                                      depth=depth,
269                                      ha=ha,
270                                      ua=ua,
271                                      va=va)
272
273        weights=num.ones(1, num.float)
274        #ensure that files are indeed mux2 files
275        times, latitudes, longitudes, elevation, stage, starttime = read_mux2_py([files[0]],weights)
276        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
277        msg='ha and ua have different gauge meta data'
278        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
279        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity, starttime_va=read_mux2_py([files[2]],weights)
280        msg='ha and va have different gauge meta data'
281        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
282
283        self.delete_mux(files)
284
285        msg='time array has incorrect length'
286        assert times.shape[0]==time_step_count,msg
287       
288        msg = 'time array is incorrect'
289        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
290        assert num.allclose(times,time_step*num.arange(time_step_count)), msg
291       
292        msg='Incorrect gauge positions returned'
293        for i,point in enumerate(lat_long_points):
294            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
295
296        msg='Incorrect gauge depths returned'
297        assert num.allclose(elevation,-depth),msg
298        msg='incorrect gauge height time series returned'
299        assert num.allclose(stage,ha)
300        msg='incorrect gauge ua time series returned'
301        assert num.allclose(xvelocity,ua)
302        msg='incorrect gauge va time series returned'
303        assert num.allclose(yvelocity, -va)
304
305    def test_urs2sts_read_mux2_pyII(self):
306        """Spatially varing stage
307        """
308        tide = 1
309        time_step_count = 3
310        time_step = 2
311        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
312        n=len(lat_long_points)
313        first_tstep=num.ones(n,num.int)
314        last_tstep=(time_step_count)*num.ones(n,num.int)
315        depth=20*num.ones(n,num.float)
316        ha=2*num.ones((n,time_step_count),num.float)
317        ha[0]=num.arange(0,time_step_count)+1
318        ha[1]=time_step_count-num.arange(1,time_step_count+1)
319        ha[1]=num.arange(time_step_count,2*time_step_count)
320        ha[2]=num.arange(2*time_step_count,3*time_step_count)
321        ha[3]=num.arange(3*time_step_count,4*time_step_count)
322        ua=5*num.ones((n,time_step_count),num.float)
323        va=-10*num.ones((n,time_step_count),num.float)
324        #-ve added to take into account mux file format where south is positive.
325        base_name, files = self.write_mux2(lat_long_points,
326                                      time_step_count, time_step,
327                                      first_tstep, last_tstep,
328                                      depth=depth,
329                                      ha=ha,
330                                      ua=ua,
331                                      va=va)
332
333        weights=num.ones(1, num.float)
334        #ensure that files are indeed mux2 files
335        times, latitudes, longitudes, elevation, stage,starttime=read_mux2_py([files[0]],weights)
336        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity,starttime_ua=read_mux2_py([files[1]],weights)
337        msg='ha and ua have different gauge meta data'
338        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
339        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
340        msg='ha and va have different gauge meta data'
341        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
342
343
344        self.delete_mux(files)
345
346        msg='time array has incorrect length'
347        #assert times.shape[0]==time_step_count,msg
348        msg = 'time array is incorrect'
349        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
350        msg='Incorrect gauge positions returned'
351        for i,point in enumerate(lat_long_points):
352            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
353
354        msg='Incorrect gauge depths returned'
355        assert num.allclose(elevation, -depth),msg
356        msg='incorrect gauge height time series returned'
357        assert num.allclose(stage, ha)
358        msg='incorrect gauge ua time series returned'
359        assert num.allclose(xvelocity, ua)
360        msg='incorrect gauge va time series returned'
361        assert num.allclose(yvelocity, -va) # South is positive in MUX
362
363    def test_urs2sts_read_mux2_pyIII(self):
364        """Varying start and finish times
365        """
366        tide = 1
367        time_step_count = 3
368        time_step = 2
369        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
370        n=len(lat_long_points)
371        first_tstep=num.ones(n,num.int)
372        first_tstep[0]+=1
373        first_tstep[2]+=1
374        last_tstep=(time_step_count)*num.ones(n,num.int)
375        last_tstep[0]-=1
376
377        depth=20*num.ones(n,num.float)
378        ha=2*num.ones((n,time_step_count),num.float)
379        ha[0]=num.arange(0,time_step_count)
380        ha[1]=num.arange(time_step_count,2*time_step_count)
381        ha[2]=num.arange(2*time_step_count,3*time_step_count)
382        ha[3]=num.arange(3*time_step_count,4*time_step_count)
383        ua=5*num.ones((n,time_step_count),num.float)
384        va=-10*num.ones((n,time_step_count),num.float)
385        #-ve added to take into account mux file format where south is positive.
386        base_name, files = self.write_mux2(lat_long_points,
387                                      time_step_count, time_step,
388                                      first_tstep, last_tstep,
389                                      depth=depth,
390                                      ha=ha,
391                                      ua=ua,
392                                      va=va)
393
394        weights=num.ones(1, num.float)
395        #ensure that files are indeed mux2 files
396        times, latitudes, longitudes, elevation, stage, starttime=read_mux2_py([files[0]],weights)
397        ua_times, ua_latitudes, ua_longitudes, ua_elevation, xvelocity, starttime_ua=read_mux2_py([files[1]],weights)
398        msg='ha and ua have different gauge meta data'
399        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
400        va_times, va_latitudes, va_longitudes, va_elevation, yvelocity,starttime_va=read_mux2_py([files[2]],weights)
401        msg='ha and va have different gauge meta data'
402        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
403
404        self.delete_mux(files)
405
406        msg='time array has incorrect length'
407        #assert times.shape[0]==time_step_count,msg
408        msg = 'time array is incorrect'
409        #assert allclose(times,time_step*num.arange(1,time_step_count+1)),msg
410        msg='Incorrect gauge positions returned'
411        for i,point in enumerate(lat_long_points):
412            assert num.allclose(latitudes[i],point[0]) and num.allclose(longitudes[i],point[1]),msg
413
414
415        # Set original data used to write mux file to be zero when gauges are
416        #not recdoring
417        ha[0][0]=0.0
418        ha[0][time_step_count-1]=0.0;
419        ha[2][0]=0.0;
420        ua[0][0]=0.0
421        ua[0][time_step_count-1]=0.0;
422        ua[2][0]=0.0;
423        va[0][0]=0.0
424        va[0][time_step_count-1]=0.0;
425        va[2][0]=0.0;
426        msg='Incorrect gauge depths returned'
427        assert num.allclose(elevation,-depth),msg
428        msg='incorrect gauge height time series returned'
429        assert num.allclose(stage,ha)
430        msg='incorrect gauge ua time series returned'
431        assert num.allclose(xvelocity,ua)
432        msg='incorrect gauge va time series returned'
433        assert num.allclose(yvelocity, -va) # South is positive in mux
434       
435
436       
437    def test_read_mux_platform_problem1(self):
438        """test_read_mux_platform_problem1
439       
440        This is to test a situation where read_mux returned
441        wrong values Win32
442
443        This test passes on Windows but test_read_mux_platform_problem2
444        does not
445        """
446       
447        from urs_ext import read_mux2
448       
449        verbose = False
450               
451        tide = 1.5
452        time_step_count = 10
453        time_step = 0.2
454        times_ref = num.arange(0, time_step_count*time_step, time_step)
455
456        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
457        n = len(lat_long_points)
458       
459        # Create different timeseries starting and ending at different times
460        first_tstep=num.ones(n, num.int)
461        first_tstep[0]+=2   # Point 0 starts at 2
462        first_tstep[1]+=4   # Point 1 starts at 4       
463        first_tstep[2]+=3   # Point 2 starts at 3
464       
465        last_tstep=(time_step_count)*num.ones(n,num.int)
466        last_tstep[0]-=1    # Point 0 ends 1 step early
467        last_tstep[1]-=2    # Point 1 ends 2 steps early               
468        last_tstep[4]-=3    # Point 4 ends 3 steps early       
469       
470        # Create varying elevation data (positive values for seafloor)
471        gauge_depth=20*num.ones(n,num.float)
472        for i in range(n):
473            gauge_depth[i] += i**2
474           
475        # Create data to be written to first mux file       
476        ha0=2*num.ones((n,time_step_count),num.float)
477        ha0[0]=num.arange(0,time_step_count)
478        ha0[1]=num.arange(time_step_count,2*time_step_count)
479        ha0[2]=num.arange(2*time_step_count,3*time_step_count)
480        ha0[3]=num.arange(3*time_step_count,4*time_step_count)
481        ua0=5*num.ones((n,time_step_count),num.float)
482        va0=-10*num.ones((n,time_step_count),num.float)
483
484        # Ensure data used to write mux file to be zero when gauges are
485        # not recording
486        for i in range(n):
487             # For each point
488             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
489                 # For timesteps before and after recording range
490                 ha0[i][j] = ua0[i][j] = va0[i][j] = 0.0                                 
491       
492        # Write first mux file to be combined by urs2sts
493        base_nameI, filesI = self.write_mux2(lat_long_points,
494                                             time_step_count, time_step,
495                                             first_tstep, last_tstep,
496                                             depth=gauge_depth,
497                                             ha=ha0,
498                                             ua=ua0,
499                                             va=va0)
500
501        # Create ordering file
502        permutation = ensure_numeric([4,0,2])
503
504        _, ordering_filename = tempfile.mkstemp('')
505        order_fid = open(ordering_filename, 'w') 
506        order_fid.write('index, longitude, latitude\n')
507        for index in permutation:
508            order_fid.write('%d, %f, %f\n' %(index, 
509                                             lat_long_points[index][1], 
510                                             lat_long_points[index][0]))
511        order_fid.close()
512       
513       
514
515        # -------------------------------------
516        # Now read files back and check values
517        weights = ensure_numeric([1.0])
518
519        # For each quantity read the associated list of source mux2 file with
520        # extention associated with that quantity
521        file_params=-1*num.ones(3,num.float) #[nsta,dt,nt]
522        OFFSET = 5
523
524        for j, file in enumerate(filesI):
525            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
526
527            number_of_selected_stations = data.shape[0]
528
529            # Index where data ends and parameters begin
530            parameters_index = data.shape[1]-OFFSET         
531         
532            for i in range(number_of_selected_stations):
533                if j == 0: assert num.allclose(data[i][:parameters_index], ha0[permutation[i], :])
534                if j == 1: assert num.allclose(data[i][:parameters_index], ua0[permutation[i], :])
535                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
536       
537        self.delete_mux(filesI)
538
539
540       
541       
542    def test_read_mux_platform_problem2(self):
543        """test_read_mux_platform_problem2
544       
545        This is to test a situation where read_mux returned
546        wrong values Win32
547
548        This test does not pass on Windows but test_read_mux_platform_problem1
549        does
550        """
551       
552        from urs_ext import read_mux2
553       
554        from anuga.config import single_precision as epsilon       
555       
556        verbose = False
557               
558        tide = 1.5
559        time_step_count = 10
560        time_step = 0.2
561       
562        times_ref = num.arange(0, time_step_count*time_step, time_step)
563       
564        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
565                           (-21.,115.), (-22., 117.)]
566        n = len(lat_long_points)
567       
568        # Create different timeseries starting and ending at different times
569        first_tstep=num.ones(n,num.int)
570        first_tstep[0]+=2   # Point 0 starts at 2
571        first_tstep[1]+=4   # Point 1 starts at 4       
572        first_tstep[2]+=3   # Point 2 starts at 3
573       
574        last_tstep=(time_step_count)*num.ones(n,num.int)
575        last_tstep[0]-=1    # Point 0 ends 1 step early
576        last_tstep[1]-=2    # Point 1 ends 2 steps early               
577        last_tstep[4]-=3    # Point 4 ends 3 steps early       
578       
579        # Create varying elevation data (positive values for seafloor)
580        gauge_depth=20*num.ones(n,num.float)
581        for i in range(n):
582            gauge_depth[i] += i**2
583           
584        # Create data to be written to second mux file       
585        ha1=num.ones((n,time_step_count),num.float)
586        ha1[0]=num.sin(times_ref)
587        ha1[1]=2*num.sin(times_ref - 3)
588        ha1[2]=5*num.sin(4*times_ref)
589        ha1[3]=num.sin(times_ref)
590        ha1[4]=num.sin(2*times_ref-0.7)
591               
592        ua1=num.zeros((n,time_step_count),num.float)
593        ua1[0]=3*num.cos(times_ref)       
594        ua1[1]=2*num.sin(times_ref-0.7)   
595        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
596        ua1[4]=2*num.ones(time_step_count)
597       
598        va1=num.zeros((n,time_step_count),num.float)
599        va1[0]=2*num.cos(times_ref-0.87)       
600        va1[1]=3*num.ones(time_step_count)
601        va1[3]=2*num.sin(times_ref-0.71)       
602       
603        # Ensure data used to write mux file to be zero when gauges are
604        # not recording
605        for i in range(n):
606             # For each point
607             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
608                 # For timesteps before and after recording range
609                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 
610
611
612        #print 'Second station to be written to MUX'
613        #print 'ha', ha1[0,:]
614        #print 'ua', ua1[0,:]
615        #print 'va', va1[0,:]
616       
617        # Write second mux file to be combined by urs2sts
618        base_nameII, filesII = self.write_mux2(lat_long_points,
619                                               time_step_count, time_step,
620                                               first_tstep, last_tstep,
621                                               depth=gauge_depth,
622                                               ha=ha1,
623                                               ua=ua1,
624                                               va=va1)
625
626
627
628
629        # Read mux file back and verify it's correcness
630
631        ####################################################
632        # FIXME (Ole): This is where the test should
633        # verify that the MUX files are correct.
634
635        #JJ: It appears as though
636        #that certain quantities are not being stored with enough precision
637        #inn muxfile or more likely that they are being cast into a
638        #lower precision when read in using read_mux2 Time step and q_time
639        # are equal but only to approx 1e-7
640        ####################################################
641
642        #define information as it should be stored in mus2 files
643        points_num=len(lat_long_points)
644        depth=gauge_depth
645        ha=ha1
646        ua=ua1
647        va=va1
648       
649        quantities = ['HA','UA','VA']
650        mux_names = [WAVEHEIGHT_MUX2_LABEL,
651                     EAST_VELOCITY_MUX2_LABEL,
652                     NORTH_VELOCITY_MUX2_LABEL]
653        quantities_init = [[],[],[]]
654        latlondeps = []
655        #irrelevant header information
656        ig=ilon=ilat=0
657        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
658        # urs binary is latitude fastest
659        for i,point in enumerate(lat_long_points):
660            lat = point[0]
661            lon = point[1]
662            _ , e, n = redfearn(lat, lon)
663            if depth is None:
664                this_depth = n
665            else:
666                this_depth = depth[i]
667            latlondeps.append([lat, lon, this_depth])
668
669            if ha is None:
670                this_ha = e
671                quantities_init[0].append(num.ones(time_step_count,num.float)*this_ha) # HA
672            else:
673                quantities_init[0].append(ha[i])
674            if ua is None:
675                this_ua = n
676                quantities_init[1].append(num.ones(time_step_count,num.float)*this_ua) # UA
677            else:
678                quantities_init[1].append(ua[i])
679            if va is None:
680                this_va = e
681                quantities_init[2].append(num.ones(time_step_count,num.float)*this_va) #
682            else:
683                quantities_init[2].append(va[i])
684
685        for i, q in enumerate(quantities):
686            #print
687            #print i, q
688           
689            q_time = num.zeros((time_step_count, points_num), num.float)
690            quantities_init[i] = ensure_numeric(quantities_init[i])
691            for time in range(time_step_count):
692                #print i, q, time, quantities_init[i][:,time]
693                q_time[time,:] = quantities_init[i][:,time]
694                #print i, q, time, q_time[time, :]
695
696           
697            filename = base_nameII + mux_names[i]
698            f = open(filename, 'rb')
699            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
700            #write mux 2 header
701            for latlondep in latlondeps:
702                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
703                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
704                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
705                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
706                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
707                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
708                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
709                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
710                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
711                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
712                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
713                assert abs(az-unpack('f',f.read(4))[0])<epsilon
714                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
715               
716                x = unpack('f', f.read(4))[0]
717                #print time_step
718                #print x
719                assert abs(time_step-x)<epsilon
720                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
721                for j in range(4): # identifier
722                    assert abs(id-unpack('i',f.read(4))[0])<epsilon
723
724            #first_tstep=1
725            #last_tstep=time_step_count
726            for i,latlondep in enumerate(latlondeps):
727                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
728            for i,latlondep in enumerate(latlondeps):
729                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
730
731            # Find when first station starts recording
732            min_tstep = min(first_tstep)
733            # Find when all stations have stopped recording
734            max_tstep = max(last_tstep)
735
736            #for time in  range(time_step_count):
737            for time in range(min_tstep-1,max_tstep):
738                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
739                for point_i in range(points_num):
740                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
741                        x = unpack('f',f.read(4))[0]
742                        #print time, x, q_time[time, point_i]
743                        if q == 'VA': x = -x # South is positive in MUX
744                        assert abs(q_time[time, point_i]-x)<epsilon
745
746            f.close()
747                                               
748        # Create ordering file
749        permutation = ensure_numeric([4,0,2])
750
751       #  _, ordering_filename = tempfile.mkstemp('')
752#         order_fid = open(ordering_filename, 'w') 
753#         order_fid.write('index, longitude, latitude\n')
754#         for index in permutation:
755#             order_fid.write('%d, %f, %f\n' %(index,
756#                                              lat_long_points[index][1],
757#                                              lat_long_points[index][0]))
758#         order_fid.close()
759       
760        # -------------------------------------
761        # Now read files back and check values
762        weights = ensure_numeric([1.0])
763
764        # For each quantity read the associated list of source mux2 file with
765        # extention associated with that quantity
766        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
767        OFFSET = 5
768
769        for j, file in enumerate(filesII):
770            # Read stage, u, v enumerated as j
771            #print 'Reading', j, file
772            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
773
774            #print 'Data received by Python'
775            #print data[1][8]
776            number_of_selected_stations = data.shape[0]
777
778            # Index where data ends and parameters begin
779            parameters_index = data.shape[1]-OFFSET         
780                 
781            quantity=num.zeros((number_of_selected_stations, parameters_index), num.float)
782           
783           
784            for i in range(number_of_selected_stations):
785       
786                #print i, parameters_index
787                #print quantity[i][:]
788                if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :])
789                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
790                if j == 2:
791                    # FIXME (Ole): This is where the output is wrong on Win32
792                   
793                    #print
794                    #print j, i
795                    #print 'Input'
796                    #print 'u', ua1[permutation[i], 8]       
797                    #print 'v', va1[permutation[i], 8]
798               
799                    #print 'Output'
800                    #print 'v ', data[i][:parameters_index][8] 
801
802                    # South is positive in MUX
803                    #print "data[i][:parameters_index]", data[i][:parameters_index]
804                    #print "-va1[permutation[i], :]", -va1[permutation[i], :]
805                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
806       
807        self.delete_mux(filesII)
808           
809    def test_read_mux_platform_problem3(self):
810       
811        # This is to test a situation where read_mux returned
812        # wrong values Win32
813
814       
815        from urs_ext import read_mux2
816       
817        from anuga.config import single_precision as epsilon       
818       
819        verbose = False
820               
821        tide = 1.5
822        time_step_count = 10
823        time_step = 0.02
824
825        '''
826        Win results
827        time_step = 0.2000001
828        This is OK       
829        '''
830       
831        '''
832        Win results
833        time_step = 0.20000001
834
835        ======================================================================
836ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager)
837----------------------------------------------------------------------
838Traceback (most recent call last):
839  File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3
840    ha1[0]=num.sin(times_ref)
841ValueError: matrices are not aligned for copy
842
843        '''
844
845        '''
846        Win results
847        time_step = 0.200000001
848        FAIL
849         assert num.allclose(data[i][:parameters_index],
850         -va1[permutation[i], :])
851        '''
852        times_ref = num.arange(0, time_step_count*time_step, time_step)
853        #print "times_ref", times_ref
854       
855        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
856                           (-21.,115.), (-22., 117.)]
857        stations = len(lat_long_points)
858       
859        # Create different timeseries starting and ending at different times
860        first_tstep=num.ones(stations, num.int)
861        first_tstep[0]+=2   # Point 0 starts at 2
862        first_tstep[1]+=4   # Point 1 starts at 4       
863        first_tstep[2]+=3   # Point 2 starts at 3
864       
865        last_tstep=(time_step_count)*num.ones(stations, num.int)
866        last_tstep[0]-=1    # Point 0 ends 1 step early
867        last_tstep[1]-=2    # Point 1 ends 2 steps early               
868        last_tstep[4]-=3    # Point 4 ends 3 steps early       
869       
870        # Create varying elevation data (positive values for seafloor)
871        gauge_depth=20*num.ones(stations, num.float)
872        for i in range(stations):
873            gauge_depth[i] += i**2
874           
875        # Create data to be written to second mux file       
876        ha1=num.ones((stations,time_step_count), num.float)
877        ha1[0]=num.sin(times_ref)
878        ha1[1]=2*num.sin(times_ref - 3)
879        ha1[2]=5*num.sin(4*times_ref)
880        ha1[3]=num.sin(times_ref)
881        ha1[4]=num.sin(2*times_ref-0.7)
882               
883        ua1=num.zeros((stations,time_step_count),num.float)
884        ua1[0]=3*num.cos(times_ref)       
885        ua1[1]=2*num.sin(times_ref-0.7)   
886        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
887        ua1[4]=2*num.ones(time_step_count)
888       
889        va1=num.zeros((stations,time_step_count),num.float)
890        va1[0]=2*num.cos(times_ref-0.87)       
891        va1[1]=3*num.ones(time_step_count)
892        va1[3]=2*num.sin(times_ref-0.71)       
893        #print "va1[0]", va1[0]  # The 8th element is what will go bad.
894        # Ensure data used to write mux file to be zero when gauges are
895        # not recording
896        for i in range(stations):
897             # For each point
898             for j in range(0, first_tstep[i]-1) + range(last_tstep[i],
899                                                         time_step_count):
900                 # For timesteps before and after recording range
901                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0 
902
903
904        #print 'Second station to be written to MUX'
905        #print 'ha', ha1[0,:]
906        #print 'ua', ua1[0,:]
907        #print 'va', va1[0,:]
908       
909        # Write second mux file to be combined by urs2sts
910        base_nameII, filesII = self.write_mux2(lat_long_points,
911                                               time_step_count, time_step,
912                                               first_tstep, last_tstep,
913                                               depth=gauge_depth,
914                                               ha=ha1,
915                                               ua=ua1,
916                                               va=va1)
917        #print "filesII", filesII
918
919
920
921
922        # Read mux file back and verify it's correcness
923
924        ####################################################
925        # FIXME (Ole): This is where the test should
926        # verify that the MUX files are correct.
927
928        #JJ: It appears as though
929        #that certain quantities are not being stored with enough precision
930        #inn muxfile or more likely that they are being cast into a
931        #lower precision when read in using read_mux2 Time step and q_time
932        # are equal but only to approx 1e-7
933        ####################################################
934
935        #define information as it should be stored in mus2 files
936        points_num=len(lat_long_points)
937        depth=gauge_depth
938        ha=ha1
939        ua=ua1
940        va=va1
941       
942        quantities = ['HA','UA','VA']
943        mux_names = [WAVEHEIGHT_MUX2_LABEL,
944                     EAST_VELOCITY_MUX2_LABEL,
945                     NORTH_VELOCITY_MUX2_LABEL]
946        quantities_init = [[],[],[]]
947        latlondeps = []
948        #irrelevant header information
949        ig=ilon=ilat=0
950        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
951        # urs binary is latitude fastest
952        for i,point in enumerate(lat_long_points):
953            lat = point[0]
954            lon = point[1]
955            _ , e, n = redfearn(lat, lon)
956            if depth is None:
957                this_depth = n
958            else:
959                this_depth = depth[i]
960            latlondeps.append([lat, lon, this_depth])
961
962            if ha is None:
963                this_ha = e
964                quantities_init[0].append(num.ones(time_step_count,
965                                                   num.float)*this_ha) # HA
966            else:
967                quantities_init[0].append(ha[i])
968            if ua is None:
969                this_ua = n
970                quantities_init[1].append(num.ones(time_step_count,
971                                                   num.float)*this_ua) # UA
972            else:
973                quantities_init[1].append(ua[i])
974            if va is None:
975                this_va = e
976                quantities_init[2].append(num.ones(time_step_count,
977                                                   num.float)*this_va) #
978            else:
979                quantities_init[2].append(va[i])
980
981        for i, q in enumerate(quantities):
982            #print
983            #print i, q
984           
985            q_time = num.zeros((time_step_count, points_num), num.float)
986            quantities_init[i] = ensure_numeric(quantities_init[i])
987            for time in range(time_step_count):
988                #print i, q, time, quantities_init[i][:,time]
989                q_time[time,:] = quantities_init[i][:,time]
990                #print i, q, time, q_time[time, :]
991
992           
993            filename = base_nameII + mux_names[i]
994            f = open(filename, 'rb')
995            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
996            #write mux 2 header
997            for latlondep in latlondeps:
998                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
999                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
1000                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
1001                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
1002                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
1003                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
1004                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
1005                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
1006                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
1007                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
1008                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
1009                assert abs(az-unpack('f',f.read(4))[0])<epsilon
1010                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
1011               
1012                x = unpack('f', f.read(4))[0]
1013                #print time_step
1014                #print x
1015                assert abs(time_step-x)<epsilon
1016                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
1017                for j in range(4): # identifier
1018                    assert abs(id-unpack('i',f.read(4))[0])<epsilon
1019
1020            #first_tstep=1
1021            #last_tstep=time_step_count
1022            for i,latlondep in enumerate(latlondeps):
1023                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
1024            for i,latlondep in enumerate(latlondeps):
1025                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
1026
1027            # Find when first station starts recording
1028            min_tstep = min(first_tstep)
1029            # Find when all stations have stopped recording
1030            max_tstep = max(last_tstep)
1031
1032            #for time in  range(time_step_count):
1033            for time in range(min_tstep-1,max_tstep):
1034                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
1035                for point_i in range(points_num):
1036                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
1037                        x = unpack('f',f.read(4))[0]
1038                        #print time, x, q_time[time, point_i]
1039                        if q == 'VA': x = -x # South is positive in MUX
1040                        #print q+" q_time[%d, %d] = %f" %(time, point_i,
1041                                                      #q_time[time, point_i])
1042                        assert abs(q_time[time, point_i]-x)<epsilon
1043
1044            f.close()
1045                           
1046        permutation = ensure_numeric([4,0,2])
1047                   
1048        # Create ordering file
1049#         _, ordering_filename = tempfile.mkstemp('')
1050#         order_fid = open(ordering_filename, 'w') 
1051#         order_fid.write('index, longitude, latitude\n')
1052#         for index in permutation:
1053#             order_fid.write('%d, %f, %f\n' %(index,
1054#                                              lat_long_points[index][1],
1055#                                              lat_long_points[index][0]))
1056#         order_fid.close()
1057       
1058        # -------------------------------------
1059        # Now read files back and check values
1060        weights = ensure_numeric([1.0])
1061
1062        # For each quantity read the associated list of source mux2 file with
1063        # extention associated with that quantity
1064        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
1065        OFFSET = 5
1066
1067        for j, file in enumerate(filesII):
1068            # Read stage, u, v enumerated as j
1069            #print 'Reading', j, file
1070            #print "file", file
1071            data = read_mux2(1, [file], weights, file_params,
1072                             permutation, verbose)
1073            #print str(j) + "data", data
1074
1075            #print 'Data received by Python'
1076            #print data[1][8]
1077            number_of_selected_stations = data.shape[0]
1078            #print "number_of_selected_stations", number_of_selected_stations
1079            #print "stations", stations
1080
1081            # Index where data ends and parameters begin
1082            parameters_index = data.shape[1]-OFFSET         
1083                 
1084            for i in range(number_of_selected_stations):
1085       
1086                #print i, parameters_index
1087                if j == 0:
1088                    assert num.allclose(data[i][:parameters_index],
1089                                        ha1[permutation[i], :])
1090                   
1091                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
1092                if j == 2:
1093                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
1094       
1095        self.delete_mux(filesII)     
1096       
1097
1098################################################################################
1099
1100if __name__ == "__main__":
1101    suite = unittest.makeSuite(TestCase,'test')
1102    runner = unittest.TextTestRunner() #verbosity=2)
1103    runner.run(suite)
1104       
Note: See TracBrowser for help on using the repository browser.