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

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

File tests passing OK - extra module for .urs files.

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