source: inundation/pyvolution/test_data_manager.py @ 2639

Last change on this file since 2639 was 2639, checked in by ole, 18 years ago

Added new test for ferret generated NetCDF. This test breaks deliberately, since the code has yet to be written.

File size: 136.7 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import copy
6from Numeric import zeros, array, allclose, Float
7from utilities.numerical_tools import mean
8import tempfile
9import os
10from Scientific.IO.NetCDF import NetCDFFile
11
12from data_manager import *
13from shallow_water import *
14from config import epsilon
15import data_manager
16
17from coordinate_transforms.geo_reference import Geo_reference
18
19class Test_Data_Manager(unittest.TestCase):
20    def setUp(self):
21        import time
22        from mesh_factory import rectangular
23
24        #Create basic mesh
25        points, vertices, boundary = rectangular(2, 2)
26
27        #Create shallow water domain
28        domain = Domain(points, vertices, boundary)
29        domain.default_order=2
30
31
32        #Set some field values
33        domain.set_quantity('elevation', lambda x,y: -x)
34        domain.set_quantity('friction', 0.03)
35
36
37        ######################
38        # Boundary conditions
39        B = Transmissive_boundary(domain)
40        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
41
42
43        ######################
44        #Initial condition - with jumps
45
46
47        bed = domain.quantities['elevation'].vertex_values
48        stage = zeros(bed.shape, Float)
49
50        h = 0.3
51        for i in range(stage.shape[0]):
52            if i % 2 == 0:
53                stage[i,:] = bed[i,:] + h
54            else:
55                stage[i,:] = bed[i,:]
56
57        domain.set_quantity('stage', stage)
58        self.initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
59
60        domain.distribute_to_vertices_and_edges()
61
62
63        self.domain = domain
64
65        C = domain.get_vertex_coordinates()
66        self.X = C[:,0:6:2].copy()
67        self.Y = C[:,1:6:2].copy()
68
69        self.F = bed
70
71
72
73
74        #Write A testfile (not realistic. Values aren't realistic)
75
76        self.test_MOST_file = 'most_small'
77
78        longitudes = [150.66667, 150.83334, 151., 151.16667]
79        latitudes = [-34.5, -34.33333, -34.16667, -34]
80
81        long_name = 'LON'
82        lat_name = 'LAT'
83
84        nx = 4
85        ny = 4
86        six = 6
87
88
89        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
90            fid = NetCDFFile(self.test_MOST_file + ext, 'w')
91       
92            fid.createDimension(long_name,nx)
93            fid.createVariable(long_name,'d',(long_name,))
94            fid.variables[long_name].point_spacing='uneven'
95            fid.variables[long_name].units='degrees_east'
96            fid.variables[long_name].assignValue(longitudes)
97
98            fid.createDimension(lat_name,ny)
99            fid.createVariable(lat_name,'d',(lat_name,))
100            fid.variables[lat_name].point_spacing='uneven'
101            fid.variables[lat_name].units='degrees_north'
102            fid.variables[lat_name].assignValue(latitudes)
103
104            fid.createDimension('TIME',six)
105            fid.createVariable('TIME','d',('TIME',))
106            fid.variables['TIME'].point_spacing='uneven'
107            fid.variables['TIME'].units='seconds'
108            fid.variables['TIME'].assignValue([0.0, 0.1, 0.6, 1.1, 1.6, 2.1])
109           
110
111            name = ext[1:3].upper()
112            if name == 'E.': name = 'ELEVATION'
113            fid.createVariable(name,'d',('TIME', lat_name, long_name))
114            fid.variables[name].units='CENTIMETERS'
115            fid.variables[name].missing_value=-1.e+034       
116           
117            fid.variables[name].assignValue([[[0.3400644, 0, -46.63519, -6.50198],
118                                              [-0.1214216, 0, 0, 0],
119                                              [0, 0, 0, 0],
120                                              [0, 0, 0, 0]],
121                                             [[0.3400644, 2.291054e-005, -23.33335, -6.50198],
122                                              [-0.1213987, 4.581959e-005, -1.594838e-007, 1.421085e-012],
123                                              [2.291054e-005, 4.582107e-005, 4.581715e-005, 1.854517e-009],
124                                              [0, 2.291054e-005, 2.291054e-005, 0]],
125                                             [[0.3400644, 0.0001374632, -23.31503, -6.50198],
126                                              [-0.1212842, 0.0002756907, 0.006325484, 1.380492e-006],
127                                              [0.0001374632, 0.0002749264, 0.0002742863, 6.665601e-008],
128                                              [0, 0.0001374632, 0.0001374632, 0]],
129                                             [[0.3400644, 0.0002520159, -23.29672, -6.50198],
130                                              [-0.1211696, 0.0005075303, 0.01264618, 6.208276e-006],
131                                              [0.0002520159, 0.0005040318, 0.0005027961, 2.23865e-007],
132                                              [0, 0.0002520159, 0.0002520159, 0]],
133                                             [[0.3400644, 0.0003665686, -23.27842, -6.50198],
134                                              [-0.1210551, 0.0007413362, 0.01896192, 1.447638e-005],
135                                              [0.0003665686, 0.0007331371, 0.0007313463, 4.734126e-007],
136                                              [0, 0.0003665686, 0.0003665686, 0]],
137                                             [[0.3400644, 0.0004811212, -23.26012, -6.50198],
138                                              [-0.1209405, 0.0009771062, 0.02527271, 2.617787e-005],
139                                              [0.0004811212, 0.0009622425, 0.0009599366, 8.152277e-007],
140                                              [0, 0.0004811212, 0.0004811212, 0]]])
141           
142       
143            fid.close()
144
145       
146
147
148    def tearDown(self):
149        import os
150        for ext in ['_ha.nc', '_ua.nc', '_va.nc', '_e.nc']:
151            #print 'Trying to remove', self.test_MOST_file + ext
152            os.remove(self.test_MOST_file + ext)       
153
154    def test_sww_constant(self):
155        """Test that constant sww information can be written correctly
156        (non smooth)
157        """
158
159        import time, os
160        from Numeric import array, zeros, allclose, Float, concatenate
161        from Scientific.IO.NetCDF import NetCDFFile
162
163        self.domain.filename = 'datatest' + str(id(self))
164        self.domain.format = 'sww'
165        self.domain.smooth = False
166
167        sww = get_dataobject(self.domain)
168        sww.store_connectivity()
169
170        #Check contents
171        #Get NetCDF
172        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
173
174        # Get the variables
175        x = fid.variables['x']
176        y = fid.variables['y']
177        z = fid.variables['elevation']
178
179        volumes = fid.variables['volumes']
180
181
182        assert allclose (x[:], self.X.flat)
183        assert allclose (y[:], self.Y.flat)
184        assert allclose (z[:], self.F.flat)
185
186        V = volumes
187
188        P = len(self.domain)
189        for k in range(P):
190            assert V[k, 0] == 3*k
191            assert V[k, 1] == 3*k+1
192            assert V[k, 2] == 3*k+2
193
194
195        fid.close()
196
197        #Cleanup
198        os.remove(sww.filename)
199
200
201    def test_sww_constant_smooth(self):
202        """Test that constant sww information can be written correctly
203        (non smooth)
204        """
205
206        import time, os
207        from Numeric import array, zeros, allclose, Float, concatenate
208        from Scientific.IO.NetCDF import NetCDFFile
209
210        self.domain.filename = 'datatest' + str(id(self))
211        self.domain.format = 'sww'
212        self.domain.smooth = True
213
214        sww = get_dataobject(self.domain)
215        sww.store_connectivity()
216
217        #Check contents
218        #Get NetCDF
219        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
220
221        # Get the variables
222        x = fid.variables['x']
223        y = fid.variables['y']
224        z = fid.variables['elevation']
225
226        volumes = fid.variables['volumes']
227
228        X = x[:]
229        Y = y[:]
230
231        assert allclose([X[0], Y[0]], array([0.0, 0.0]))
232        assert allclose([X[1], Y[1]], array([0.0, 0.5]))
233        assert allclose([X[2], Y[2]], array([0.0, 1.0]))
234
235        assert allclose([X[4], Y[4]], array([0.5, 0.5]))
236
237        assert allclose([X[7], Y[7]], array([1.0, 0.5]))
238
239        Z = z[:]
240        assert Z[4] == -0.5
241
242        V = volumes
243        assert V[2,0] == 4
244        assert V[2,1] == 5
245        assert V[2,2] == 1
246
247        assert V[4,0] == 6
248        assert V[4,1] == 7
249        assert V[4,2] == 3
250
251
252        fid.close()
253
254        #Cleanup
255        os.remove(sww.filename)
256
257
258
259    def test_sww_variable(self):
260        """Test that sww information can be written correctly
261        """
262
263        import time, os
264        from Numeric import array, zeros, allclose, Float, concatenate
265        from Scientific.IO.NetCDF import NetCDFFile
266
267        self.domain.filename = 'datatest' + str(id(self))
268        self.domain.format = 'sww'
269        self.domain.smooth = True
270        self.domain.reduction = mean
271
272        sww = get_dataobject(self.domain)
273        sww.store_connectivity()
274        sww.store_timestep('stage')
275
276        #Check contents
277        #Get NetCDF
278        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
279
280
281        # Get the variables
282        x = fid.variables['x']
283        y = fid.variables['y']
284        z = fid.variables['elevation']
285        time = fid.variables['time']
286        stage = fid.variables['stage']
287
288
289        Q = self.domain.quantities['stage']
290        Q0 = Q.vertex_values[:,0]
291        Q1 = Q.vertex_values[:,1]
292        Q2 = Q.vertex_values[:,2]
293
294        A = stage[0,:]
295        #print A[0], (Q2[0,0] + Q1[1,0])/2
296        assert allclose(A[0], (Q2[0] + Q1[1])/2)
297        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
298        assert allclose(A[2], Q0[3])
299        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
300
301        #Center point
302        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
303                                 Q0[5] + Q2[6] + Q1[7])/6)
304
305
306
307        fid.close()
308
309        #Cleanup
310        os.remove(sww.filename)
311
312
313    def test_sww_variable2(self):
314        """Test that sww information can be written correctly
315        multiple timesteps. Use average as reduction operator
316        """
317
318        import time, os
319        from Numeric import array, zeros, allclose, Float, concatenate
320        from Scientific.IO.NetCDF import NetCDFFile
321
322        self.domain.filename = 'datatest' + str(id(self))
323        self.domain.format = 'sww'
324        self.domain.smooth = True
325
326        self.domain.reduction = mean
327
328        sww = get_dataobject(self.domain)
329        sww.store_connectivity()
330        sww.store_timestep('stage')
331        self.domain.evolve_to_end(finaltime = 0.01)
332        sww.store_timestep('stage')
333
334
335        #Check contents
336        #Get NetCDF
337        fid = NetCDFFile(sww.filename, 'r')  #Open existing file for append
338
339        # Get the variables
340        x = fid.variables['x']
341        y = fid.variables['y']
342        z = fid.variables['elevation']
343        time = fid.variables['time']
344        stage = fid.variables['stage']
345
346        #Check values
347        Q = self.domain.quantities['stage']
348        Q0 = Q.vertex_values[:,0]
349        Q1 = Q.vertex_values[:,1]
350        Q2 = Q.vertex_values[:,2]
351
352        A = stage[1,:]
353        assert allclose(A[0], (Q2[0] + Q1[1])/2)
354        assert allclose(A[1], (Q0[1] + Q1[3] + Q2[2])/3)
355        assert allclose(A[2], Q0[3])
356        assert allclose(A[3], (Q0[0] + Q1[5] + Q2[4])/3)
357
358        #Center point
359        assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
360                                 Q0[5] + Q2[6] + Q1[7])/6)
361
362
363        fid.close()
364
365        #Cleanup
366        os.remove(sww.filename)
367
368    def test_sww_variable3(self):
369        """Test that sww information can be written correctly
370        multiple timesteps using a different reduction operator (min)
371        """
372
373        import time, os
374        from Numeric import array, zeros, allclose, Float, concatenate
375        from Scientific.IO.NetCDF import NetCDFFile
376
377        self.domain.filename = 'datatest' + str(id(self))
378        self.domain.format = 'sww'
379        self.domain.smooth = True
380        self.domain.reduction = min
381
382        sww = get_dataobject(self.domain)
383        sww.store_connectivity()
384        sww.store_timestep('stage')
385
386        self.domain.evolve_to_end(finaltime = 0.01)
387        sww.store_timestep('stage')
388
389
390        #Check contents
391        #Get NetCDF
392        fid = NetCDFFile(sww.filename, 'r')
393
394
395        # Get the variables
396        x = fid.variables['x']
397        y = fid.variables['y']
398        z = fid.variables['elevation']
399        time = fid.variables['time']
400        stage = fid.variables['stage']
401
402        #Check values
403        Q = self.domain.quantities['stage']
404        Q0 = Q.vertex_values[:,0]
405        Q1 = Q.vertex_values[:,1]
406        Q2 = Q.vertex_values[:,2]
407
408        A = stage[1,:]
409        assert allclose(A[0], min(Q2[0], Q1[1]))
410        assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
411        assert allclose(A[2], Q0[3])
412        assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
413
414        #Center point
415        assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
416                                  Q0[5], Q2[6], Q1[7]))
417
418
419        fid.close()
420
421        #Cleanup
422        os.remove(sww.filename)
423
424
425    def test_sync(self):
426        """Test info stored at each timestep is as expected (incl initial condition)
427        """
428
429        import time, os, config
430        from Numeric import array, zeros, allclose, Float, concatenate
431        from Scientific.IO.NetCDF import NetCDFFile
432
433        self.domain.filename = 'synctest'
434        self.domain.format = 'sww'
435        self.domain.smooth = False
436        self.domain.store = True
437        self.domain.beta_h = 0
438
439        #Evolution
440        for t in self.domain.evolve(yieldstep = 1.0, finaltime = 4.0):
441            stage = self.domain.quantities['stage'].vertex_values
442
443            #Get NetCDF
444            fid = NetCDFFile(self.domain.writer.filename, 'r')
445            stage_file = fid.variables['stage']
446
447            if t == 0.0:
448                assert allclose(stage, self.initial_stage)
449                assert allclose(stage_file[:], stage.flat)
450            else:
451                assert not allclose(stage, self.initial_stage)
452                assert not allclose(stage_file[:], stage.flat)
453
454            fid.close()
455
456        os.remove(self.domain.writer.filename)
457
458
459
460    def test_sww_DSG(self):
461        """Not a test, rather a look at the sww format
462        """
463
464        import time, os
465        from Numeric import array, zeros, allclose, Float, concatenate
466        from Scientific.IO.NetCDF import NetCDFFile
467
468        self.domain.filename = 'datatest' + str(id(self))
469        self.domain.format = 'sww'
470        self.domain.smooth = True
471        self.domain.reduction = mean
472
473        sww = get_dataobject(self.domain)
474        sww.store_connectivity()
475        sww.store_timestep('stage')
476
477        #Check contents
478        #Get NetCDF
479        fid = NetCDFFile(sww.filename, 'r')
480
481        # Get the variables
482        x = fid.variables['x']
483        y = fid.variables['y']
484        z = fid.variables['elevation']
485
486        volumes = fid.variables['volumes']
487        time = fid.variables['time']
488
489        # 2D
490        stage = fid.variables['stage']
491
492        X = x[:]
493        Y = y[:]
494        Z = z[:]
495        V = volumes[:]
496        T = time[:]
497        S = stage[:,:]
498
499#         print "****************************"
500#         print "X ",X
501#         print "****************************"
502#         print "Y ",Y
503#         print "****************************"
504#         print "Z ",Z
505#         print "****************************"
506#         print "V ",V
507#         print "****************************"
508#         print "Time ",T
509#         print "****************************"
510#         print "Stage ",S
511#         print "****************************"
512
513
514        fid.close()
515
516        #Cleanup
517        os.remove(sww.filename)
518
519
520    #def test_write_pts(self):
521    #    #Obsolete
522    #   
523    #    #Get (enough) datapoints
524    #
525    #    from Numeric import array
526    #    points = array([[ 0.66666667, 0.66666667],
527    #                    [ 1.33333333, 1.33333333],
528    #                    [ 2.66666667, 0.66666667],
529    #                    [ 0.66666667, 2.66666667],
530    #                    [ 0.0, 1.0],
531    #                    [ 0.0, 3.0],
532    #                    [ 1.0, 0.0],
533    #                    [ 1.0, 1.0],
534    #                    [ 1.0, 2.0],
535    #                    [ 1.0, 3.0],
536    #                    [ 2.0, 1.0],
537    #                    [ 3.0, 0.0],
538    #                    [ 3.0, 1.0]])
539    #
540    #    z = points[:,0] + 2*points[:,1]
541    #
542    #    ptsfile = 'testptsfile.pts'
543    #    write_ptsfile(ptsfile, points, z,
544    #                  attribute_name = 'linear_combination')
545    #
546    #    #Check contents
547    #    #Get NetCDF
548    #    from Scientific.IO.NetCDF import NetCDFFile       
549    #    fid = NetCDFFile(ptsfile, 'r')
550    #
551    #    # Get the variables
552    #    #print fid.variables.keys()
553    #    points1 = fid.variables['points']
554    #    z1 = fid.variables['linear_combination']
555    #
556    #    #Check values#
557    #
558    #    #print points[:]
559    #    #print ref_points
560    #    assert allclose(points, points1) 
561    #
562    #    #print attributes[:]
563    #    #print ref_elevation
564    #    assert allclose(z, z1)
565    #
566    #    #Cleanup
567    #    fid.close()
568    #
569    #    import os
570    #    os.remove(ptsfile)
571       
572
573    def test_dem2pts_bounding_box_v2(self):
574        """Test conversion from dem in ascii format to native NetCDF xya format
575        """
576
577        import time, os
578        from Numeric import array, zeros, allclose, Float, concatenate, ones
579        from Scientific.IO.NetCDF import NetCDFFile
580
581        #Write test asc file
582        root = 'demtest'
583
584        filename = root+'.asc'
585        fid = open(filename, 'w')
586        fid.write("""ncols         10
587nrows         10
588xllcorner     2000
589yllcorner     3000
590cellsize      1
591NODATA_value  -9999
592""")
593        #Create linear function       
594        ref_points = []
595        ref_elevation = []
596        x0 = 2000
597        y = 3010
598        yvec = range(10)
599        xvec = range(10)
600        z = -1
601        for i in range(10):
602            y = y - 1
603            for j in range(10):
604                x = x0 + xvec[j]
605                z += 1
606                ref_points.append ([x,y])
607                ref_elevation.append(z)
608                fid.write('%f ' %z)
609            fid.write('\n')
610
611        fid.close()
612           
613        #print 'sending pts', ref_points
614        #print 'sending elev', ref_elevation
615
616        #Write prj file with metadata
617        metafilename = root+'.prj'
618        fid = open(metafilename, 'w')
619
620
621        fid.write("""Projection UTM
622Zone 56
623Datum WGS84
624Zunits NO
625Units METERS
626Spheroid WGS84
627Xshift 0.0000000000
628Yshift 10000000.0000000000
629Parameters
630""")
631        fid.close()
632
633        #Convert to NetCDF pts
634        convert_dem_from_ascii2netcdf(root)
635        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
636                northing_min=3003.0, northing_max=3006.0)
637
638        #Check contents
639        #Get NetCDF
640        fid = NetCDFFile(root+'.pts', 'r')
641
642        # Get the variables
643        #print fid.variables.keys()
644        points = fid.variables['points']
645        elevation = fid.variables['elevation']
646
647        #Check values
648        assert fid.xllcorner[0] == 2002.0
649        assert fid.yllcorner[0] == 3003.0
650
651        #create new reference points
652        newz = []       
653        newz[0:5] = ref_elevation[32:38]
654        newz[6:11] = ref_elevation[42:48]
655        newz[12:17] = ref_elevation[52:58]
656        newz[18:23] = ref_elevation[62:68]
657        ref_elevation = []
658        ref_elevation = newz
659        ref_points = []
660        x0 = 2002
661        y = 3007
662        yvec = range(4)
663        xvec = range(6)
664        for i in range(4):
665            y = y - 1
666            ynew = y - 3003.0
667            for j in range(6):
668                x = x0 + xvec[j]
669                xnew = x - 2002.0
670                ref_points.append ([xnew,ynew]) #Relative point values       
671
672        assert allclose(points, ref_points)
673
674        assert allclose(elevation, ref_elevation)
675
676        #Cleanup
677        fid.close()
678
679
680        os.remove(root + '.pts')
681        os.remove(root + '.dem')
682        os.remove(root + '.asc')
683        os.remove(root + '.prj')
684
685
686    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
687        """Test conversion from dem in ascii format to native NetCDF xya format
688        """
689
690        import time, os
691        from Numeric import array, zeros, allclose, Float, concatenate, ones
692        from Scientific.IO.NetCDF import NetCDFFile
693
694        #Write test asc file
695        root = 'demtest'
696
697        filename = root+'.asc'
698        fid = open(filename, 'w')
699        fid.write("""ncols         10
700nrows         10
701xllcorner     2000
702yllcorner     3000
703cellsize      1
704NODATA_value  -9999
705""")
706        #Create linear function       
707        ref_points = []
708        ref_elevation = []
709        x0 = 2000
710        y = 3010
711        yvec = range(10)
712        xvec = range(10)
713        #z = range(100)
714        z = zeros(100)
715        NODATA_value = -9999
716        count = -1
717        for i in range(10):
718            y = y - 1
719            for j in range(10):
720                x = x0 + xvec[j]
721                ref_points.append ([x,y])
722                count += 1
723                z[count] = (4*i - 3*j)%13 
724                if j == 4: z[count] = NODATA_value #column inside clipping region
725                if j == 8: z[count] = NODATA_value #column outside clipping region
726                if i == 9: z[count] = NODATA_value #row outside clipping region
727                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region               
728                ref_elevation.append( z[count] )
729                fid.write('%f ' %z[count])
730            fid.write('\n')
731
732        fid.close()             
733       
734        #print 'sending elev', ref_elevation
735
736        #Write prj file with metadata
737        metafilename = root+'.prj'
738        fid = open(metafilename, 'w')
739
740
741        fid.write("""Projection UTM
742Zone 56
743Datum WGS84
744Zunits NO
745Units METERS
746Spheroid WGS84
747Xshift 0.0000000000
748Yshift 10000000.0000000000
749Parameters
750""")
751        fid.close()
752
753        #Convert to NetCDF pts
754        convert_dem_from_ascii2netcdf(root)
755        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
756                northing_min=3003.0, northing_max=3006.0)
757
758        #Check contents
759        #Get NetCDF
760        fid = NetCDFFile(root+'.pts', 'r')
761
762        # Get the variables
763        #print fid.variables.keys()
764        points = fid.variables['points']
765        elevation = fid.variables['elevation']
766
767        #Check values
768        assert fid.xllcorner[0] == 2002.0
769        assert fid.yllcorner[0] == 3003.0
770
771        #create new reference points
772        newz = zeros(19)
773        newz[0:2] = ref_elevation[32:34]
774        newz[2:5] = ref_elevation[35:38]
775        newz[5:7] = ref_elevation[42:44]
776        newz[7] = ref_elevation[45]
777        newz[8] = ref_elevation[47]
778        newz[9:11] = ref_elevation[52:54]
779        newz[11:14] = ref_elevation[55:58]
780        newz[14:16] = ref_elevation[62:64]
781        newz[16:19] = ref_elevation[65:68]
782
783       
784        ref_elevation = newz
785        ref_points = []
786        new_ref_points = []
787        x0 = 2002
788        y = 3007
789        yvec = range(4)
790        xvec = range(6)
791        for i in range(4):
792            y = y - 1
793            ynew = y - 3003.0
794            for j in range(6):
795                x = x0 + xvec[j]
796                xnew = x - 2002.0
797                if j <> 2 and (i<>1 or j<>4):   
798                    ref_points.append([x,y])
799                    new_ref_points.append ([xnew,ynew])
800
801               
802        assert allclose(points, new_ref_points)
803        assert allclose(elevation, ref_elevation)
804
805        #Cleanup
806        fid.close()
807
808
809        os.remove(root + '.pts')
810        os.remove(root + '.dem')
811        os.remove(root + '.asc')
812        os.remove(root + '.prj')
813
814
815    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
816        """Test conversion from dem in ascii format to native NetCDF xya format
817        Check missing values on clipping boundary
818        """
819
820        import time, os
821        from Numeric import array, zeros, allclose, Float, concatenate, ones
822        from Scientific.IO.NetCDF import NetCDFFile
823
824        #Write test asc file
825        root = 'demtest'
826
827        filename = root+'.asc'
828        fid = open(filename, 'w')
829        fid.write("""ncols         10
830nrows         10
831xllcorner     2000
832yllcorner     3000
833cellsize      1
834NODATA_value  -9999
835""")
836        #Create linear function       
837        ref_points = []
838        ref_elevation = []
839        x0 = 2000
840        y = 3010
841        yvec = range(10)
842        xvec = range(10)
843        #z = range(100)
844        z = zeros(100)
845        NODATA_value = -9999
846        count = -1
847        for i in range(10):
848            y = y - 1
849            for j in range(10):
850                x = x0 + xvec[j]
851                ref_points.append ([x,y])
852                count += 1
853                z[count] = (4*i - 3*j)%13 
854                if j == 4: z[count] = NODATA_value #column inside clipping region
855                if j == 8: z[count] = NODATA_value #column outside clipping region
856                if i == 6: z[count] = NODATA_value #row on clipping boundary
857                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region               
858                ref_elevation.append( z[count] )
859                fid.write('%f ' %z[count])
860            fid.write('\n')
861
862        fid.close()             
863       
864        #print 'sending elev', ref_elevation
865
866        #Write prj file with metadata
867        metafilename = root+'.prj'
868        fid = open(metafilename, 'w')
869
870
871        fid.write("""Projection UTM
872Zone 56
873Datum WGS84
874Zunits NO
875Units METERS
876Spheroid WGS84
877Xshift 0.0000000000
878Yshift 10000000.0000000000
879Parameters
880""")
881        fid.close()
882
883        #Convert to NetCDF pts
884        convert_dem_from_ascii2netcdf(root)
885        dem2pts(root, easting_min=2002.0, easting_max=2007.0,
886                northing_min=3003.0, northing_max=3006.0)
887
888        #Check contents
889        #Get NetCDF
890        fid = NetCDFFile(root+'.pts', 'r')
891
892        # Get the variables
893        #print fid.variables.keys()
894        points = fid.variables['points']
895        elevation = fid.variables['elevation']
896
897        #Check values
898        assert fid.xllcorner[0] == 2002.0
899        assert fid.yllcorner[0] == 3003.0
900
901        #create new reference points
902        newz = zeros(14)
903        newz[0:2] = ref_elevation[32:34]
904        newz[2:5] = ref_elevation[35:38]
905        newz[5:7] = ref_elevation[42:44]
906        newz[7] = ref_elevation[45]
907        newz[8] = ref_elevation[47]
908        newz[9:11] = ref_elevation[52:54]
909        newz[11:14] = ref_elevation[55:58]
910
911
912       
913        ref_elevation = newz
914        ref_points = []
915        new_ref_points = []
916        x0 = 2002
917        y = 3007
918        yvec = range(4)
919        xvec = range(6)
920        for i in range(4):
921            y = y - 1
922            ynew = y - 3003.0
923            for j in range(6):
924                x = x0 + xvec[j]
925                xnew = x - 2002.0
926                if j <> 2 and (i<>1 or j<>4) and i<>3:   
927                    ref_points.append([x,y])
928                    new_ref_points.append ([xnew,ynew])
929
930
931        #print points[:],points[:].shape
932        #print new_ref_points, len(new_ref_points)
933
934        assert allclose(elevation, ref_elevation)                   
935        assert allclose(points, new_ref_points)
936
937
938        #Cleanup
939        fid.close()
940
941
942        os.remove(root + '.pts')
943        os.remove(root + '.dem')
944        os.remove(root + '.asc')
945        os.remove(root + '.prj')
946
947
948    def test_hecras_cross_sections2pts(self):
949        """Test conversion from HECRAS cross sections in ascii format
950        to native NetCDF pts format
951        """
952
953        import time, os
954        from Numeric import array, zeros, allclose, Float, concatenate
955        from Scientific.IO.NetCDF import NetCDFFile
956
957        #Write test asc file
958        root = 'hecrastest'
959
960        filename = root+'.sdf'
961        fid = open(filename, 'w')
962        fid.write("""
963# RAS export file created on Mon 15Aug2005 11:42
964# by HEC-RAS Version 3.1.1
965
966BEGIN HEADER:
967  UNITS: METRIC
968  DTM TYPE: TIN
969  DTM: v:\1\cit\perth_topo\river_tin
970  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
971  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
972  MAP PROJECTION: UTM
973  PROJECTION ZONE: 50
974  DATUM: AGD66
975  VERTICAL DATUM:
976  NUMBER OF REACHES:  19
977  NUMBER OF CROSS-SECTIONS:  2
978END HEADER:
979
980
981BEGIN CROSS-SECTIONS:
982
983  CROSS-SECTION:
984    STREAM ID:Southern-Wungong
985    REACH ID:Southern-Wungong
986    STATION:21410   
987    CUT LINE:
988      407546.08 , 6437277.542
989      407329.32 , 6437489.482
990      407283.11 , 6437541.232
991    SURFACE LINE:
992     407546.08,   6437277.54,   52.14
993     407538.88,   6437284.58,   51.07
994     407531.68,   6437291.62,   50.56
995     407524.48,   6437298.66,   49.58
996     407517.28,   6437305.70,   49.09
997     407510.08,   6437312.74,   48.76
998  END:
999
1000  CROSS-SECTION:
1001    STREAM ID:Swan River     
1002    REACH ID:Swan Mouth     
1003    STATION:840.*   
1004    CUT LINE:
1005      381178.0855 , 6452559.0685
1006      380485.4755 , 6453169.272
1007    SURFACE LINE:
1008     381178.09,   6452559.07,   4.17
1009     381169.49,   6452566.64,   4.26
1010     381157.78,   6452576.96,   4.34
1011     381155.97,   6452578.56,   4.35
1012     381143.72,   6452589.35,   4.43
1013     381136.69,   6452595.54,   4.58
1014     381114.74,   6452614.88,   4.41
1015     381075.53,   6452649.43,   4.17
1016     381071.47,   6452653.00,   3.99
1017     381063.46,   6452660.06,   3.67
1018     381054.41,   6452668.03,   3.67
1019  END: 
1020END CROSS-SECTIONS:
1021""")
1022
1023        fid.close()
1024
1025
1026        #Convert to NetCDF pts
1027        hecras_cross_sections2pts(root)
1028
1029        #Check contents
1030        #Get NetCDF
1031        fid = NetCDFFile(root+'.pts', 'r')
1032
1033        # Get the variables
1034        #print fid.variables.keys()
1035        points = fid.variables['points']
1036        elevation = fid.variables['elevation']
1037
1038        #Check values
1039        ref_points = [[407546.08, 6437277.54],
1040                      [407538.88, 6437284.58],
1041                      [407531.68, 6437291.62],
1042                      [407524.48, 6437298.66],
1043                      [407517.28, 6437305.70],
1044                      [407510.08, 6437312.74]]
1045
1046        ref_points += [[381178.09, 6452559.07], 
1047                       [381169.49, 6452566.64], 
1048                       [381157.78, 6452576.96], 
1049                       [381155.97, 6452578.56], 
1050                       [381143.72, 6452589.35], 
1051                       [381136.69, 6452595.54], 
1052                       [381114.74, 6452614.88], 
1053                       [381075.53, 6452649.43], 
1054                       [381071.47, 6452653.00], 
1055                       [381063.46, 6452660.06], 
1056                       [381054.41, 6452668.03]] 
1057       
1058                     
1059        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
1060        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
1061
1062        #print points[:]
1063        #print ref_points
1064        assert allclose(points, ref_points)
1065
1066        #print attributes[:]
1067        #print ref_elevation
1068        assert allclose(elevation, ref_elevation)
1069
1070        #Cleanup
1071        fid.close()
1072
1073
1074        os.remove(root + '.sdf')
1075        os.remove(root + '.pts')
1076
1077
1078
1079    def test_sww2dem_asc_elevation(self):
1080        """Test that sww information can be converted correctly to asc/prj
1081        format readable by e.g. ArcView
1082        """
1083
1084        import time, os
1085        from Numeric import array, zeros, allclose, Float, concatenate
1086        from Scientific.IO.NetCDF import NetCDFFile
1087
1088        #Setup
1089        self.domain.filename = 'datatest'
1090
1091        prjfile = self.domain.filename + '_elevation.prj'
1092        ascfile = self.domain.filename + '_elevation.asc'
1093        swwfile = self.domain.filename + '.sww'
1094
1095        self.domain.set_datadir('.')
1096        self.domain.format = 'sww'
1097        self.domain.smooth = True
1098        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1099
1100        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1101
1102        sww = get_dataobject(self.domain)
1103        sww.store_connectivity()
1104        sww.store_timestep('stage')
1105
1106        self.domain.evolve_to_end(finaltime = 0.01)
1107        sww.store_timestep('stage')
1108
1109        cellsize = 0.25
1110        #Check contents
1111        #Get NetCDF
1112
1113        fid = NetCDFFile(sww.filename, 'r')
1114
1115        # Get the variables
1116        x = fid.variables['x'][:]
1117        y = fid.variables['y'][:]
1118        z = fid.variables['elevation'][:]
1119        time = fid.variables['time'][:]
1120        stage = fid.variables['stage'][:]
1121
1122
1123        #Export to ascii/prj files
1124        sww2dem(self.domain.filename,
1125                quantity = 'elevation',
1126                cellsize = cellsize,
1127                verbose = False,
1128                format = 'asc')
1129               
1130        #Check prj (meta data)
1131        prjid = open(prjfile)
1132        lines = prjid.readlines()
1133        prjid.close()
1134
1135        L = lines[0].strip().split()
1136        assert L[0].strip().lower() == 'projection'
1137        assert L[1].strip().lower() == 'utm'
1138
1139        L = lines[1].strip().split()
1140        assert L[0].strip().lower() == 'zone'
1141        assert L[1].strip().lower() == '56'
1142
1143        L = lines[2].strip().split()
1144        assert L[0].strip().lower() == 'datum'
1145        assert L[1].strip().lower() == 'wgs84'
1146
1147        L = lines[3].strip().split()
1148        assert L[0].strip().lower() == 'zunits'
1149        assert L[1].strip().lower() == 'no'
1150
1151        L = lines[4].strip().split()
1152        assert L[0].strip().lower() == 'units'
1153        assert L[1].strip().lower() == 'meters'
1154
1155        L = lines[5].strip().split()
1156        assert L[0].strip().lower() == 'spheroid'
1157        assert L[1].strip().lower() == 'wgs84'
1158
1159        L = lines[6].strip().split()
1160        assert L[0].strip().lower() == 'xshift'
1161        assert L[1].strip().lower() == '500000'
1162
1163        L = lines[7].strip().split()
1164        assert L[0].strip().lower() == 'yshift'
1165        assert L[1].strip().lower() == '10000000'
1166
1167        L = lines[8].strip().split()
1168        assert L[0].strip().lower() == 'parameters'
1169
1170
1171        #Check asc file
1172        ascid = open(ascfile)
1173        lines = ascid.readlines()
1174        ascid.close()
1175
1176        L = lines[0].strip().split()
1177        assert L[0].strip().lower() == 'ncols'
1178        assert L[1].strip().lower() == '5'
1179
1180        L = lines[1].strip().split()
1181        assert L[0].strip().lower() == 'nrows'
1182        assert L[1].strip().lower() == '5'
1183
1184        L = lines[2].strip().split()
1185        assert L[0].strip().lower() == 'xllcorner'
1186        assert allclose(float(L[1].strip().lower()), 308500)
1187
1188        L = lines[3].strip().split()
1189        assert L[0].strip().lower() == 'yllcorner'
1190        assert allclose(float(L[1].strip().lower()), 6189000)
1191
1192        L = lines[4].strip().split()
1193        assert L[0].strip().lower() == 'cellsize'
1194        assert allclose(float(L[1].strip().lower()), cellsize)
1195
1196        L = lines[5].strip().split()
1197        assert L[0].strip() == 'NODATA_value'
1198        assert L[1].strip().lower() == '-9999'
1199
1200        #Check grid values
1201        for j in range(5):
1202            L = lines[6+j].strip().split()
1203            y = (4-j) * cellsize
1204            for i in range(5):
1205                assert allclose(float(L[i]), -i*cellsize - y)
1206
1207
1208        fid.close()
1209
1210        #Cleanup
1211        os.remove(prjfile)
1212        os.remove(ascfile)
1213        os.remove(swwfile)
1214
1215
1216
1217    def test_sww2dem_larger(self):
1218        """Test that sww information can be converted correctly to asc/prj
1219        format readable by e.g. ArcView. Here:
1220
1221        ncols         11
1222        nrows         11
1223        xllcorner     308500
1224        yllcorner     6189000
1225        cellsize      10.000000
1226        NODATA_value  -9999
1227        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1228         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1229         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1230         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
1231         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
1232         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
1233         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
1234         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
1235         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
1236         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
1237           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
1238       
1239        """
1240
1241        import time, os
1242        from Numeric import array, zeros, allclose, Float, concatenate
1243        from Scientific.IO.NetCDF import NetCDFFile
1244
1245        #Setup
1246
1247        from mesh_factory import rectangular
1248
1249        #Create basic mesh (100m x 100m)
1250        points, vertices, boundary = rectangular(2, 2, 100, 100)
1251
1252        #Create shallow water domain
1253        domain = Domain(points, vertices, boundary)
1254        domain.default_order = 2
1255
1256        domain.filename = 'datatest'
1257
1258        prjfile = domain.filename + '_elevation.prj'
1259        ascfile = domain.filename + '_elevation.asc'
1260        swwfile = domain.filename + '.sww'
1261
1262        domain.set_datadir('.')
1263        domain.format = 'sww'
1264        domain.smooth = True
1265        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1266
1267        #
1268        domain.set_quantity('elevation', lambda x,y: -x-y)
1269        domain.set_quantity('stage', 0)       
1270
1271        B = Transmissive_boundary(domain)
1272        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1273
1274
1275        #
1276        sww = get_dataobject(domain)
1277        sww.store_connectivity()
1278        sww.store_timestep('stage')
1279
1280        domain.evolve_to_end(finaltime = 0.01)
1281        sww.store_timestep('stage')
1282
1283        cellsize = 10  #10m grid
1284
1285       
1286        #Check contents
1287        #Get NetCDF
1288
1289        fid = NetCDFFile(sww.filename, 'r')
1290
1291        # Get the variables
1292        x = fid.variables['x'][:]
1293        y = fid.variables['y'][:]
1294        z = fid.variables['elevation'][:]
1295        time = fid.variables['time'][:]
1296        stage = fid.variables['stage'][:]
1297
1298
1299        #Export to ascii/prj files
1300        sww2dem(domain.filename,
1301                quantity = 'elevation',
1302                cellsize = cellsize,
1303                verbose = False,
1304                format = 'asc')
1305       
1306               
1307        #Check prj (meta data)
1308        prjid = open(prjfile)
1309        lines = prjid.readlines()
1310        prjid.close()
1311
1312        L = lines[0].strip().split()
1313        assert L[0].strip().lower() == 'projection'
1314        assert L[1].strip().lower() == 'utm'
1315
1316        L = lines[1].strip().split()
1317        assert L[0].strip().lower() == 'zone'
1318        assert L[1].strip().lower() == '56'
1319
1320        L = lines[2].strip().split()
1321        assert L[0].strip().lower() == 'datum'
1322        assert L[1].strip().lower() == 'wgs84'
1323
1324        L = lines[3].strip().split()
1325        assert L[0].strip().lower() == 'zunits'
1326        assert L[1].strip().lower() == 'no'
1327
1328        L = lines[4].strip().split()
1329        assert L[0].strip().lower() == 'units'
1330        assert L[1].strip().lower() == 'meters'
1331
1332        L = lines[5].strip().split()
1333        assert L[0].strip().lower() == 'spheroid'
1334        assert L[1].strip().lower() == 'wgs84'
1335
1336        L = lines[6].strip().split()
1337        assert L[0].strip().lower() == 'xshift'
1338        assert L[1].strip().lower() == '500000'
1339
1340        L = lines[7].strip().split()
1341        assert L[0].strip().lower() == 'yshift'
1342        assert L[1].strip().lower() == '10000000'
1343
1344        L = lines[8].strip().split()
1345        assert L[0].strip().lower() == 'parameters'
1346
1347
1348        #Check asc file
1349        ascid = open(ascfile)
1350        lines = ascid.readlines()
1351        ascid.close()
1352
1353        L = lines[0].strip().split()
1354        assert L[0].strip().lower() == 'ncols'
1355        assert L[1].strip().lower() == '11'
1356
1357        L = lines[1].strip().split()
1358        assert L[0].strip().lower() == 'nrows'
1359        assert L[1].strip().lower() == '11'
1360
1361        L = lines[2].strip().split()
1362        assert L[0].strip().lower() == 'xllcorner'
1363        assert allclose(float(L[1].strip().lower()), 308500)
1364
1365        L = lines[3].strip().split()
1366        assert L[0].strip().lower() == 'yllcorner'
1367        assert allclose(float(L[1].strip().lower()), 6189000)
1368
1369        L = lines[4].strip().split()
1370        assert L[0].strip().lower() == 'cellsize'
1371        assert allclose(float(L[1].strip().lower()), cellsize)
1372
1373        L = lines[5].strip().split()
1374        assert L[0].strip() == 'NODATA_value'
1375        assert L[1].strip().lower() == '-9999'
1376
1377        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1378        for i, line in enumerate(lines[6:]):
1379            for j, value in enumerate( line.split() ):
1380                #assert allclose(float(value), -(10-i+j)*cellsize)
1381                assert float(value) == -(10-i+j)*cellsize
1382       
1383
1384        fid.close()
1385
1386        #Cleanup
1387        os.remove(prjfile)
1388        os.remove(ascfile)
1389        os.remove(swwfile)
1390
1391
1392
1393    def test_sww2dem_boundingbox(self):
1394        """Test that sww information can be converted correctly to asc/prj
1395        format readable by e.g. ArcView.
1396        This will test that mesh can be restricted by bounding box
1397
1398        Original extent is 100m x 100m:
1399
1400        Eastings:   308500 -  308600
1401        Northings: 6189000 - 6189100
1402
1403        Bounding box changes this to the 50m x 50m square defined by
1404
1405        Eastings:   308530 -  308570
1406        Northings: 6189050 - 6189100
1407
1408        The cropped values should be
1409
1410         -130 -140 -150 -160 -170 
1411         -120 -130 -140 -150 -160 
1412         -110 -120 -130 -140 -150 
1413         -100 -110 -120 -130 -140 
1414          -90 -100 -110 -120 -130 
1415          -80  -90 -100 -110 -120
1416
1417        and the new lower reference point should be   
1418        Eastings:   308530
1419        Northings: 6189050
1420       
1421        Original dataset is the same as in test_sww2dem_larger()
1422       
1423        """
1424
1425        import time, os
1426        from Numeric import array, zeros, allclose, Float, concatenate
1427        from Scientific.IO.NetCDF import NetCDFFile
1428
1429        #Setup
1430
1431        from mesh_factory import rectangular
1432
1433        #Create basic mesh (100m x 100m)
1434        points, vertices, boundary = rectangular(2, 2, 100, 100)
1435
1436        #Create shallow water domain
1437        domain = Domain(points, vertices, boundary)
1438        domain.default_order = 2
1439
1440        domain.filename = 'datatest'
1441
1442        prjfile = domain.filename + '_elevation.prj'
1443        ascfile = domain.filename + '_elevation.asc'
1444        swwfile = domain.filename + '.sww'
1445
1446        domain.set_datadir('.')
1447        domain.format = 'sww'
1448        domain.smooth = True
1449        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1450
1451        #
1452        domain.set_quantity('elevation', lambda x,y: -x-y)
1453        domain.set_quantity('stage', 0)       
1454
1455        B = Transmissive_boundary(domain)
1456        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1457
1458
1459        #
1460        sww = get_dataobject(domain)
1461        sww.store_connectivity()
1462        sww.store_timestep('stage')
1463
1464        domain.evolve_to_end(finaltime = 0.01)
1465        sww.store_timestep('stage')
1466
1467        cellsize = 10  #10m grid
1468
1469       
1470        #Check contents
1471        #Get NetCDF
1472
1473        fid = NetCDFFile(sww.filename, 'r')
1474
1475        # Get the variables
1476        x = fid.variables['x'][:]
1477        y = fid.variables['y'][:]
1478        z = fid.variables['elevation'][:]
1479        time = fid.variables['time'][:]
1480        stage = fid.variables['stage'][:]
1481
1482
1483        #Export to ascii/prj files
1484        sww2dem(domain.filename,
1485                quantity = 'elevation',
1486                cellsize = cellsize,
1487                easting_min = 308530,
1488                easting_max = 308570,
1489                northing_min = 6189050, 
1490                northing_max = 6189100,                               
1491                verbose = False,
1492                format = 'asc')
1493       
1494        fid.close()
1495
1496       
1497        #Check prj (meta data)
1498        prjid = open(prjfile)
1499        lines = prjid.readlines()
1500        prjid.close()
1501
1502        L = lines[0].strip().split()
1503        assert L[0].strip().lower() == 'projection'
1504        assert L[1].strip().lower() == 'utm'
1505
1506        L = lines[1].strip().split()
1507        assert L[0].strip().lower() == 'zone'
1508        assert L[1].strip().lower() == '56'
1509
1510        L = lines[2].strip().split()
1511        assert L[0].strip().lower() == 'datum'
1512        assert L[1].strip().lower() == 'wgs84'
1513
1514        L = lines[3].strip().split()
1515        assert L[0].strip().lower() == 'zunits'
1516        assert L[1].strip().lower() == 'no'
1517
1518        L = lines[4].strip().split()
1519        assert L[0].strip().lower() == 'units'
1520        assert L[1].strip().lower() == 'meters'
1521
1522        L = lines[5].strip().split()
1523        assert L[0].strip().lower() == 'spheroid'
1524        assert L[1].strip().lower() == 'wgs84'
1525
1526        L = lines[6].strip().split()
1527        assert L[0].strip().lower() == 'xshift'
1528        assert L[1].strip().lower() == '500000'
1529
1530        L = lines[7].strip().split()
1531        assert L[0].strip().lower() == 'yshift'
1532        assert L[1].strip().lower() == '10000000'
1533
1534        L = lines[8].strip().split()
1535        assert L[0].strip().lower() == 'parameters'
1536
1537
1538        #Check asc file
1539        ascid = open(ascfile)
1540        lines = ascid.readlines()
1541        ascid.close()
1542
1543        L = lines[0].strip().split()
1544        assert L[0].strip().lower() == 'ncols'
1545        assert L[1].strip().lower() == '5'
1546
1547        L = lines[1].strip().split()
1548        assert L[0].strip().lower() == 'nrows'
1549        assert L[1].strip().lower() == '6'
1550
1551        L = lines[2].strip().split()
1552        assert L[0].strip().lower() == 'xllcorner'
1553        assert allclose(float(L[1].strip().lower()), 308530)
1554
1555        L = lines[3].strip().split()
1556        assert L[0].strip().lower() == 'yllcorner'
1557        assert allclose(float(L[1].strip().lower()), 6189050)
1558
1559        L = lines[4].strip().split()
1560        assert L[0].strip().lower() == 'cellsize'
1561        assert allclose(float(L[1].strip().lower()), cellsize)
1562
1563        L = lines[5].strip().split()
1564        assert L[0].strip() == 'NODATA_value'
1565        assert L[1].strip().lower() == '-9999'
1566
1567        #Check grid values
1568        for i, line in enumerate(lines[6:]):
1569            for j, value in enumerate( line.split() ):
1570                #assert float(value) == -(10-i+j)*cellsize               
1571                assert float(value) == -(10-i+j+3)*cellsize
1572       
1573
1574
1575        #Cleanup
1576        os.remove(prjfile)
1577        os.remove(ascfile)
1578        os.remove(swwfile)
1579
1580
1581
1582    def test_sww2dem_asc_stage_reduction(self):
1583        """Test that sww information can be converted correctly to asc/prj
1584        format readable by e.g. ArcView
1585
1586        This tests the reduction of quantity stage using min
1587        """
1588
1589        import time, os
1590        from Numeric import array, zeros, allclose, Float, concatenate
1591        from Scientific.IO.NetCDF import NetCDFFile
1592
1593        #Setup
1594        self.domain.filename = 'datatest'
1595
1596        prjfile = self.domain.filename + '_stage.prj'
1597        ascfile = self.domain.filename + '_stage.asc'
1598        swwfile = self.domain.filename + '.sww'
1599
1600        self.domain.set_datadir('.')
1601        self.domain.format = 'sww'
1602        self.domain.smooth = True
1603        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1604
1605        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1606
1607
1608        sww = get_dataobject(self.domain)
1609        sww.store_connectivity()
1610        sww.store_timestep('stage')
1611
1612        self.domain.evolve_to_end(finaltime = 0.01)
1613        sww.store_timestep('stage')
1614
1615        cellsize = 0.25
1616        #Check contents
1617        #Get NetCDF
1618
1619        fid = NetCDFFile(sww.filename, 'r')
1620
1621        # Get the variables
1622        x = fid.variables['x'][:]
1623        y = fid.variables['y'][:]
1624        z = fid.variables['elevation'][:]
1625        time = fid.variables['time'][:]
1626        stage = fid.variables['stage'][:]
1627
1628
1629        #Export to ascii/prj files
1630        sww2dem(self.domain.filename,
1631                quantity = 'stage',
1632                cellsize = cellsize,
1633                reduction = min,
1634                format = 'asc')
1635
1636
1637        #Check asc file
1638        ascid = open(ascfile)
1639        lines = ascid.readlines()
1640        ascid.close()
1641
1642        L = lines[0].strip().split()
1643        assert L[0].strip().lower() == 'ncols'
1644        assert L[1].strip().lower() == '5'
1645
1646        L = lines[1].strip().split()
1647        assert L[0].strip().lower() == 'nrows'
1648        assert L[1].strip().lower() == '5'
1649
1650        L = lines[2].strip().split()
1651        assert L[0].strip().lower() == 'xllcorner'
1652        assert allclose(float(L[1].strip().lower()), 308500)
1653
1654        L = lines[3].strip().split()
1655        assert L[0].strip().lower() == 'yllcorner'
1656        assert allclose(float(L[1].strip().lower()), 6189000)
1657
1658        L = lines[4].strip().split()
1659        assert L[0].strip().lower() == 'cellsize'
1660        assert allclose(float(L[1].strip().lower()), cellsize)
1661
1662        L = lines[5].strip().split()
1663        assert L[0].strip() == 'NODATA_value'
1664        assert L[1].strip().lower() == '-9999'
1665
1666
1667        #Check grid values (where applicable)
1668        for j in range(5):
1669            if j%2 == 0:
1670                L = lines[6+j].strip().split()
1671                jj = 4-j
1672                for i in range(5):
1673                    if i%2 == 0:
1674                        index = jj/2 + i/2*3
1675                        val0 = stage[0,index]
1676                        val1 = stage[1,index]
1677
1678                        #print i, j, index, ':', L[i], val0, val1
1679                        assert allclose(float(L[i]), min(val0, val1))
1680
1681
1682        fid.close()
1683
1684        #Cleanup
1685        os.remove(prjfile)
1686        os.remove(ascfile)
1687        #os.remove(swwfile)
1688
1689
1690
1691    def test_sww2dem_asc_derived_quantity(self):
1692        """Test that sww information can be converted correctly to asc/prj
1693        format readable by e.g. ArcView
1694
1695        This tests the use of derived quantities
1696        """
1697
1698        import time, os
1699        from Numeric import array, zeros, allclose, Float, concatenate
1700        from Scientific.IO.NetCDF import NetCDFFile
1701
1702        #Setup
1703        self.domain.filename = 'datatest'
1704
1705        prjfile = self.domain.filename + '_depth.prj'
1706        ascfile = self.domain.filename + '_depth.asc'
1707        swwfile = self.domain.filename + '.sww'
1708
1709        self.domain.set_datadir('.')
1710        self.domain.format = 'sww'
1711        self.domain.smooth = True
1712        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1713        self.domain.set_quantity('stage', 0.0)
1714
1715        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1716
1717
1718        sww = get_dataobject(self.domain)
1719        sww.store_connectivity()
1720        sww.store_timestep('stage')
1721
1722        self.domain.evolve_to_end(finaltime = 0.01)
1723        sww.store_timestep('stage')
1724
1725        cellsize = 0.25
1726        #Check contents
1727        #Get NetCDF
1728
1729        fid = NetCDFFile(sww.filename, 'r')
1730
1731        # Get the variables
1732        x = fid.variables['x'][:]
1733        y = fid.variables['y'][:]
1734        z = fid.variables['elevation'][:]
1735        time = fid.variables['time'][:]
1736        stage = fid.variables['stage'][:]
1737
1738
1739        #Export to ascii/prj files
1740        sww2dem(self.domain.filename,
1741                basename_out = 'datatest_depth',
1742                quantity = 'stage - elevation',
1743                cellsize = cellsize,
1744                reduction = min,
1745                format = 'asc',
1746                verbose = False)
1747
1748
1749        #Check asc file
1750        ascid = open(ascfile)
1751        lines = ascid.readlines()
1752        ascid.close()
1753
1754        L = lines[0].strip().split()
1755        assert L[0].strip().lower() == 'ncols'
1756        assert L[1].strip().lower() == '5'
1757
1758        L = lines[1].strip().split()
1759        assert L[0].strip().lower() == 'nrows'
1760        assert L[1].strip().lower() == '5'
1761
1762        L = lines[2].strip().split()
1763        assert L[0].strip().lower() == 'xllcorner'
1764        assert allclose(float(L[1].strip().lower()), 308500)
1765
1766        L = lines[3].strip().split()
1767        assert L[0].strip().lower() == 'yllcorner'
1768        assert allclose(float(L[1].strip().lower()), 6189000)
1769
1770        L = lines[4].strip().split()
1771        assert L[0].strip().lower() == 'cellsize'
1772        assert allclose(float(L[1].strip().lower()), cellsize)
1773
1774        L = lines[5].strip().split()
1775        assert L[0].strip() == 'NODATA_value'
1776        assert L[1].strip().lower() == '-9999'
1777
1778
1779        #Check grid values (where applicable)
1780        for j in range(5):
1781            if j%2 == 0:
1782                L = lines[6+j].strip().split()
1783                jj = 4-j
1784                for i in range(5):
1785                    if i%2 == 0:
1786                        index = jj/2 + i/2*3
1787                        val0 = stage[0,index] - z[index]
1788                        val1 = stage[1,index] - z[index]
1789
1790                        #print i, j, index, ':', L[i], val0, val1
1791                        assert allclose(float(L[i]), min(val0, val1))
1792
1793
1794        fid.close()
1795
1796        #Cleanup
1797        os.remove(prjfile)
1798        os.remove(ascfile)
1799        #os.remove(swwfile)
1800
1801
1802
1803
1804
1805    def test_sww2dem_asc_missing_points(self):
1806        """Test that sww information can be converted correctly to asc/prj
1807        format readable by e.g. ArcView
1808
1809        This test includes the writing of missing values
1810        """
1811
1812        import time, os
1813        from Numeric import array, zeros, allclose, Float, concatenate
1814        from Scientific.IO.NetCDF import NetCDFFile
1815
1816        #Setup mesh not coinciding with rectangle.
1817        #This will cause missing values to occur in gridded data
1818
1819
1820        points = [                        [1.0, 1.0],
1821                              [0.5, 0.5], [1.0, 0.5],
1822                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1823
1824        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1825
1826        #Create shallow water domain
1827        domain = Domain(points, vertices)
1828        domain.default_order=2
1829
1830
1831        #Set some field values
1832        domain.set_quantity('elevation', lambda x,y: -x-y)
1833        domain.set_quantity('friction', 0.03)
1834
1835
1836        ######################
1837        # Boundary conditions
1838        B = Transmissive_boundary(domain)
1839        domain.set_boundary( {'exterior': B} )
1840
1841
1842        ######################
1843        #Initial condition - with jumps
1844
1845        bed = domain.quantities['elevation'].vertex_values
1846        stage = zeros(bed.shape, Float)
1847
1848        h = 0.3
1849        for i in range(stage.shape[0]):
1850            if i % 2 == 0:
1851                stage[i,:] = bed[i,:] + h
1852            else:
1853                stage[i,:] = bed[i,:]
1854
1855        domain.set_quantity('stage', stage)
1856        domain.distribute_to_vertices_and_edges()
1857
1858        domain.filename = 'datatest'
1859
1860        prjfile = domain.filename + '_elevation.prj'
1861        ascfile = domain.filename + '_elevation.asc'
1862        swwfile = domain.filename + '.sww'
1863
1864        domain.set_datadir('.')
1865        domain.format = 'sww'
1866        domain.smooth = True
1867
1868        domain.geo_reference = Geo_reference(56,308500,6189000)
1869
1870        sww = get_dataobject(domain)
1871        sww.store_connectivity()
1872        sww.store_timestep('stage')
1873
1874        cellsize = 0.25
1875        #Check contents
1876        #Get NetCDF
1877
1878        fid = NetCDFFile(swwfile, 'r')
1879
1880        # Get the variables
1881        x = fid.variables['x'][:]
1882        y = fid.variables['y'][:]
1883        z = fid.variables['elevation'][:]
1884        time = fid.variables['time'][:]
1885
1886        try:
1887            geo_reference = Geo_reference(NetCDFObject=fid)
1888        except AttributeError, e:
1889            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1890
1891        #Export to ascii/prj files
1892        sww2dem(domain.filename,
1893                quantity = 'elevation',
1894                cellsize = cellsize,
1895                verbose = False,
1896                format = 'asc')
1897
1898
1899        #Check asc file
1900        ascid = open(ascfile)
1901        lines = ascid.readlines()
1902        ascid.close()
1903
1904        L = lines[0].strip().split()
1905        assert L[0].strip().lower() == 'ncols'
1906        assert L[1].strip().lower() == '5'
1907
1908        L = lines[1].strip().split()
1909        assert L[0].strip().lower() == 'nrows'
1910        assert L[1].strip().lower() == '5'
1911
1912        L = lines[2].strip().split()
1913        assert L[0].strip().lower() == 'xllcorner'
1914        assert allclose(float(L[1].strip().lower()), 308500)
1915
1916        L = lines[3].strip().split()
1917        assert L[0].strip().lower() == 'yllcorner'
1918        assert allclose(float(L[1].strip().lower()), 6189000)
1919
1920        L = lines[4].strip().split()
1921        assert L[0].strip().lower() == 'cellsize'
1922        assert allclose(float(L[1].strip().lower()), cellsize)
1923
1924        L = lines[5].strip().split()
1925        assert L[0].strip() == 'NODATA_value'
1926        assert L[1].strip().lower() == '-9999'
1927
1928        #Check grid values
1929        for j in range(5):
1930            L = lines[6+j].strip().split()
1931            assert len(L) == 5
1932            y = (4-j) * cellsize
1933
1934            for i in range(5):
1935                #print i
1936                if i+j >= 4:
1937                    assert allclose(float(L[i]), -i*cellsize - y)
1938                else:
1939                    #Missing values
1940                    assert allclose(float(L[i]), -9999)
1941
1942
1943
1944        fid.close()
1945
1946        #Cleanup
1947        os.remove(prjfile)
1948        os.remove(ascfile)
1949        os.remove(swwfile)
1950
1951    def test_sww2ers_simple(self):
1952        """Test that sww information can be converted correctly to asc/prj
1953        format readable by e.g. ArcView
1954        """
1955
1956        import time, os
1957        from Numeric import array, zeros, allclose, Float, concatenate
1958        from Scientific.IO.NetCDF import NetCDFFile
1959
1960
1961        NODATA_value = 1758323
1962
1963        #Setup
1964        self.domain.filename = 'datatest'
1965
1966        headerfile = self.domain.filename + '.ers'
1967        swwfile = self.domain.filename + '.sww'
1968
1969        self.domain.set_datadir('.')
1970        self.domain.format = 'sww'
1971        self.domain.smooth = True
1972        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1973
1974        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1975
1976        sww = get_dataobject(self.domain)
1977        sww.store_connectivity()
1978        sww.store_timestep('stage')
1979
1980        self.domain.evolve_to_end(finaltime = 0.01)
1981        sww.store_timestep('stage')
1982
1983        cellsize = 0.25
1984        #Check contents
1985        #Get NetCDF
1986
1987        fid = NetCDFFile(sww.filename, 'r')
1988
1989        # Get the variables
1990        x = fid.variables['x'][:]
1991        y = fid.variables['y'][:]
1992        z = fid.variables['elevation'][:]
1993        time = fid.variables['time'][:]
1994        stage = fid.variables['stage'][:]
1995
1996
1997        #Export to ers files
1998        #sww2ers(self.domain.filename,
1999        #        quantity = 'elevation',
2000        #        cellsize = cellsize,
2001        #        verbose = False)
2002               
2003        sww2dem(self.domain.filename,
2004                quantity = 'elevation',
2005                cellsize = cellsize,
2006                NODATA_value = NODATA_value,
2007                verbose = False,
2008                format = 'ers')
2009
2010        #Check header data
2011        from ermapper_grids import read_ermapper_header, read_ermapper_data
2012       
2013        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
2014        #print header
2015        assert header['projection'].lower() == '"utm-56"'
2016        assert header['datum'].lower() == '"wgs84"'
2017        assert header['units'].lower() == '"meters"'   
2018        assert header['value'].lower() == '"elevation"'         
2019        assert header['xdimension'] == '0.25'
2020        assert header['ydimension'] == '0.25'   
2021        assert float(header['eastings']) == 308500.0   #xllcorner
2022        assert float(header['northings']) == 6189000.0 #yllcorner       
2023        assert int(header['nroflines']) == 5
2024        assert int(header['nrofcellsperline']) == 5     
2025        assert int(header['nullcellvalue']) == NODATA_value
2026        #FIXME - there is more in the header                   
2027
2028               
2029        #Check grid data               
2030        grid = read_ermapper_data(self.domain.filename + '_elevation') 
2031       
2032        #FIXME (Ole): Why is this the desired reference grid for -x-y?
2033        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
2034                    -1,    -1.25, -1.5,  -1.75, -2.0,
2035                    -0.75, -1.0,  -1.25, -1.5,  -1.75,             
2036                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
2037                    -0.25, -0.5,  -0.75, -1.0,  -1.25]             
2038                                         
2039
2040        #print grid
2041        assert allclose(grid, ref_grid)
2042
2043        fid.close()
2044       
2045        #Cleanup
2046        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
2047        #Done (Ole) - it was because sww2ers didn't close it's sww file
2048        os.remove(sww.filename)
2049        os.remove(self.domain.filename + '_elevation')
2050        os.remove(self.domain.filename + '_elevation.ers')                 
2051
2052
2053
2054    def test_ferret2sww1(self):
2055        """Test that georeferencing etc works when converting from
2056        ferret format (lat/lon) to sww format (UTM)
2057        """
2058        from Scientific.IO.NetCDF import NetCDFFile
2059        import os, sys
2060
2061        #The test file has
2062        # LON = 150.66667, 150.83334, 151, 151.16667
2063        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2064        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2065        #
2066        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2067        # Fourth value (index==3) is -6.50198 cm
2068
2069       
2070
2071        #Read
2072        from coordinate_transforms.redfearn import redfearn
2073        #fid = NetCDFFile(self.test_MOST_file)
2074        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')       
2075        first_value = fid.variables['HA'][:][0,0,0]
2076        fourth_value = fid.variables['HA'][:][0,0,3]
2077        fid.close()
2078
2079
2080        #Call conversion (with zero origin)
2081        #ferret2sww('small', verbose=False,
2082        #           origin = (56, 0, 0))
2083        ferret2sww(self.test_MOST_file, verbose=False,
2084                   origin = (56, 0, 0))
2085
2086        #Work out the UTM coordinates for first point
2087        zone, e, n = redfearn(-34.5, 150.66667)
2088        #print zone, e, n
2089
2090        #Read output file 'small.sww'
2091        #fid = NetCDFFile('small.sww')
2092        fid = NetCDFFile(self.test_MOST_file + '.sww')       
2093
2094        x = fid.variables['x'][:]
2095        y = fid.variables['y'][:]
2096
2097        #Check that first coordinate is correctly represented
2098        assert allclose(x[0], e)
2099        assert allclose(y[0], n)
2100
2101        #Check first value
2102        stage = fid.variables['stage'][:]
2103        xmomentum = fid.variables['xmomentum'][:]
2104        ymomentum = fid.variables['ymomentum'][:]
2105
2106        #print ymomentum
2107
2108        assert allclose(stage[0,0], first_value/100)  #Meters
2109
2110        #Check fourth value
2111        assert allclose(stage[0,3], fourth_value/100)  #Meters
2112
2113        fid.close()
2114
2115        #Cleanup
2116        import os
2117        os.remove(self.test_MOST_file + '.sww')
2118
2119
2120    def test_ferret2sww_2(self):
2121        """Test that georeferencing etc works when converting from
2122        ferret format (lat/lon) to sww format (UTM)
2123        """
2124        from Scientific.IO.NetCDF import NetCDFFile
2125
2126        #The test file has
2127        # LON = 150.66667, 150.83334, 151, 151.16667
2128        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2129        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2130        #
2131        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2132        # Fourth value (index==3) is -6.50198 cm
2133
2134
2135        from coordinate_transforms.redfearn import redfearn
2136
2137        #fid = NetCDFFile('small_ha.nc')
2138        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2139
2140        #Pick a coordinate and a value
2141
2142        time_index = 1
2143        lat_index = 0
2144        lon_index = 2
2145
2146        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
2147        test_time = fid.variables['TIME'][:][time_index]
2148        test_lat = fid.variables['LAT'][:][lat_index]
2149        test_lon = fid.variables['LON'][:][lon_index]
2150
2151        linear_point_index = lat_index*4 + lon_index
2152        fid.close()
2153
2154        #Call conversion (with zero origin)
2155        ferret2sww(self.test_MOST_file, verbose=False,
2156                   origin = (56, 0, 0))
2157
2158
2159        #Work out the UTM coordinates for test point
2160        zone, e, n = redfearn(test_lat, test_lon)
2161
2162        #Read output file 'small.sww'
2163        fid = NetCDFFile(self.test_MOST_file + '.sww')
2164
2165        x = fid.variables['x'][:]
2166        y = fid.variables['y'][:]
2167
2168        #Check that test coordinate is correctly represented
2169        assert allclose(x[linear_point_index], e)
2170        assert allclose(y[linear_point_index], n)
2171
2172        #Check test value
2173        stage = fid.variables['stage'][:]
2174
2175        assert allclose(stage[time_index, linear_point_index], test_value/100)
2176
2177        fid.close()
2178
2179        #Cleanup
2180        import os
2181        os.remove(self.test_MOST_file + '.sww')
2182
2183
2184
2185    def test_ferret2sww3(self):
2186        """Elevation included
2187        """
2188        from Scientific.IO.NetCDF import NetCDFFile
2189
2190        #The test file has
2191        # LON = 150.66667, 150.83334, 151, 151.16667
2192        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2193        # ELEVATION = [-1 -2 -3 -4
2194        #             -5 -6 -7 -8
2195        #              ...
2196        #              ...      -16]
2197        # where the top left corner is -1m,
2198        # and the ll corner is -13.0m
2199        #
2200        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2201        # Fourth value (index==3) is -6.50198 cm
2202
2203        from coordinate_transforms.redfearn import redfearn
2204        import os
2205        fid1 = NetCDFFile('test_ha.nc','w')
2206        fid2 = NetCDFFile('test_ua.nc','w')
2207        fid3 = NetCDFFile('test_va.nc','w')
2208        fid4 = NetCDFFile('test_e.nc','w')
2209
2210        h1_list = [150.66667,150.83334,151.]
2211        h2_list = [-34.5,-34.33333]
2212
2213        long_name = 'LON'
2214        lat_name = 'LAT'
2215        time_name = 'TIME'       
2216
2217        nx = 3
2218        ny = 2
2219
2220        for fid in [fid1,fid2,fid3]:
2221            fid.createDimension(long_name,nx)
2222            fid.createVariable(long_name,'d',(long_name,))
2223            fid.variables[long_name].point_spacing='uneven'
2224            fid.variables[long_name].units='degrees_east'
2225            fid.variables[long_name].assignValue(h1_list)
2226
2227            fid.createDimension(lat_name,ny)
2228            fid.createVariable(lat_name,'d',(lat_name,))
2229            fid.variables[lat_name].point_spacing='uneven'
2230            fid.variables[lat_name].units='degrees_north'
2231            fid.variables[lat_name].assignValue(h2_list)
2232
2233            fid.createDimension(time_name,2)
2234            fid.createVariable(time_name,'d',(time_name,))
2235            fid.variables[time_name].point_spacing='uneven'
2236            fid.variables[time_name].units='seconds'
2237            fid.variables[time_name].assignValue([0.,1.])
2238            if fid == fid3: break
2239
2240
2241        for fid in [fid4]:
2242            fid.createDimension(long_name,nx)
2243            fid.createVariable(long_name,'d',(long_name,))
2244            fid.variables[long_name].point_spacing='uneven'
2245            fid.variables[long_name].units='degrees_east'
2246            fid.variables[long_name].assignValue(h1_list)
2247
2248            fid.createDimension(lat_name,ny)
2249            fid.createVariable(lat_name,'d',(lat_name,))
2250            fid.variables[lat_name].point_spacing='uneven'
2251            fid.variables[lat_name].units='degrees_north'
2252            fid.variables[lat_name].assignValue(h2_list)
2253
2254        name = {}
2255        name[fid1]='HA'
2256        name[fid2]='UA'
2257        name[fid3]='VA'
2258        name[fid4]='ELEVATION'
2259
2260        units = {}
2261        units[fid1]='cm'
2262        units[fid2]='cm/s'
2263        units[fid3]='cm/s'
2264        units[fid4]='m'
2265
2266        values = {}
2267        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2268        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2269        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2270        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2271
2272        for fid in [fid1,fid2,fid3]:
2273          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
2274          fid.variables[name[fid]].point_spacing='uneven'
2275          fid.variables[name[fid]].units=units[fid]
2276          fid.variables[name[fid]].assignValue(values[fid])
2277          fid.variables[name[fid]].missing_value = -99999999.
2278          if fid == fid3: break
2279
2280        for fid in [fid4]:
2281            fid.createVariable(name[fid],'d',(lat_name,long_name))
2282            fid.variables[name[fid]].point_spacing='uneven'
2283            fid.variables[name[fid]].units=units[fid]
2284            fid.variables[name[fid]].assignValue(values[fid])
2285            fid.variables[name[fid]].missing_value = -99999999.
2286
2287
2288        fid1.sync(); fid1.close()
2289        fid2.sync(); fid2.close()
2290        fid3.sync(); fid3.close()
2291        fid4.sync(); fid4.close()
2292
2293        fid1 = NetCDFFile('test_ha.nc','r')
2294        fid2 = NetCDFFile('test_e.nc','r')
2295        fid3 = NetCDFFile('test_va.nc','r')
2296
2297
2298        first_amp = fid1.variables['HA'][:][0,0,0]
2299        third_amp = fid1.variables['HA'][:][0,0,2]
2300        first_elevation = fid2.variables['ELEVATION'][0,0]
2301        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2302        first_speed = fid3.variables['VA'][0,0,0]
2303        third_speed = fid3.variables['VA'][:][0,0,2]
2304
2305        fid1.close()
2306        fid2.close()
2307        fid3.close()
2308
2309        #Call conversion (with zero origin)
2310        ferret2sww('test', verbose=False,
2311                   origin = (56, 0, 0))
2312
2313        os.remove('test_va.nc')
2314        os.remove('test_ua.nc')
2315        os.remove('test_ha.nc')
2316        os.remove('test_e.nc')
2317
2318        #Read output file 'test.sww'
2319        fid = NetCDFFile('test.sww')
2320
2321
2322        #Check first value
2323        elevation = fid.variables['elevation'][:]
2324        stage = fid.variables['stage'][:]
2325        xmomentum = fid.variables['xmomentum'][:]
2326        ymomentum = fid.variables['ymomentum'][:]
2327
2328        #print ymomentum
2329        first_height = first_amp/100 - first_elevation
2330        third_height = third_amp/100 - third_elevation
2331        first_momentum=first_speed*first_height/100
2332        third_momentum=third_speed*third_height/100
2333
2334        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2335        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2336
2337        fid.close()
2338
2339        #Cleanup
2340        os.remove('test.sww')
2341
2342
2343
2344    def test_ferret2sww4(self):
2345        """Like previous but with augmented variable names as
2346        in files produced by ferret as opposed to MOST
2347        """
2348        from Scientific.IO.NetCDF import NetCDFFile
2349
2350        #The test file has
2351        # LON = 150.66667, 150.83334, 151, 151.16667
2352        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2353        # ELEVATION = [-1 -2 -3 -4
2354        #             -5 -6 -7 -8
2355        #              ...
2356        #              ...      -16]
2357        # where the top left corner is -1m,
2358        # and the ll corner is -13.0m
2359        #
2360        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2361        # Fourth value (index==3) is -6.50198 cm
2362
2363        from coordinate_transforms.redfearn import redfearn
2364        import os
2365        fid1 = NetCDFFile('test_ha.nc','w')
2366        fid2 = NetCDFFile('test_ua.nc','w')
2367        fid3 = NetCDFFile('test_va.nc','w')
2368        fid4 = NetCDFFile('test_e.nc','w')
2369
2370        h1_list = [150.66667,150.83334,151.]
2371        h2_list = [-34.5,-34.33333]
2372
2373        long_name = 'LON961_1261'
2374        lat_name = 'LAT481_841'
2375        time_name = 'TIME1'
2376
2377        nx = 3
2378        ny = 2
2379
2380        for fid in [fid1,fid2,fid3]:
2381            fid.createDimension(long_name,nx)
2382            fid.createVariable(long_name,'d',(long_name,))
2383            fid.variables[long_name].point_spacing='uneven'
2384            fid.variables[long_name].units='degrees_east'
2385            fid.variables[long_name].assignValue(h1_list)
2386
2387            fid.createDimension(lat_name,ny)
2388            fid.createVariable(lat_name,'d',(lat_name,))
2389            fid.variables[lat_name].point_spacing='uneven'
2390            fid.variables[lat_name].units='degrees_north'
2391            fid.variables[lat_name].assignValue(h2_list)
2392
2393            fid.createDimension(time_name,2)
2394            fid.createVariable(time_name,'d',(time_name,))
2395            fid.variables[time_name].point_spacing='uneven'
2396            fid.variables[time_name].units='seconds'
2397            fid.variables[time_name].assignValue([0.,1.])
2398            if fid == fid3: break
2399
2400
2401        for fid in [fid4]:
2402            fid.createDimension(long_name,nx)
2403            fid.createVariable(long_name,'d',(long_name,))
2404            fid.variables[long_name].point_spacing='uneven'
2405            fid.variables[long_name].units='degrees_east'
2406            fid.variables[long_name].assignValue(h1_list)
2407
2408            fid.createDimension(lat_name,ny)
2409            fid.createVariable(lat_name,'d',(lat_name,))
2410            fid.variables[lat_name].point_spacing='uneven'
2411            fid.variables[lat_name].units='degrees_north'
2412            fid.variables[lat_name].assignValue(h2_list)
2413
2414        name = {}
2415        name[fid1]='HA'
2416        name[fid2]='UA'
2417        name[fid3]='VA'
2418        name[fid4]='ELEVATION'
2419
2420        units = {}
2421        units[fid1]='cm'
2422        units[fid2]='cm/s'
2423        units[fid3]='cm/s'
2424        units[fid4]='m'
2425
2426        values = {}
2427        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2428        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2429        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2430        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2431
2432        for fid in [fid1,fid2,fid3]:
2433          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
2434          fid.variables[name[fid]].point_spacing='uneven'
2435          fid.variables[name[fid]].units=units[fid]
2436          fid.variables[name[fid]].assignValue(values[fid])
2437          fid.variables[name[fid]].missing_value = -99999999.
2438          if fid == fid3: break
2439
2440        for fid in [fid4]:
2441            fid.createVariable(name[fid],'d',(lat_name,long_name))
2442            fid.variables[name[fid]].point_spacing='uneven'
2443            fid.variables[name[fid]].units=units[fid]
2444            fid.variables[name[fid]].assignValue(values[fid])
2445            fid.variables[name[fid]].missing_value = -99999999.
2446
2447
2448        fid1.sync(); fid1.close()
2449        fid2.sync(); fid2.close()
2450        fid3.sync(); fid3.close()
2451        fid4.sync(); fid4.close()
2452
2453        fid1 = NetCDFFile('test_ha.nc','r')
2454        fid2 = NetCDFFile('test_e.nc','r')
2455        fid3 = NetCDFFile('test_va.nc','r')
2456
2457
2458        first_amp = fid1.variables['HA'][:][0,0,0]
2459        third_amp = fid1.variables['HA'][:][0,0,2]
2460        first_elevation = fid2.variables['ELEVATION'][0,0]
2461        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2462        first_speed = fid3.variables['VA'][0,0,0]
2463        third_speed = fid3.variables['VA'][:][0,0,2]
2464
2465        fid1.close()
2466        fid2.close()
2467        fid3.close()
2468
2469        #Call conversion (with zero origin)
2470        ferret2sww('test', verbose=False,
2471                   origin = (56, 0, 0))
2472
2473        os.remove('test_va.nc')
2474        os.remove('test_ua.nc')
2475        os.remove('test_ha.nc')
2476        os.remove('test_e.nc')
2477
2478        #Read output file 'test.sww'
2479        fid = NetCDFFile('test.sww')
2480
2481
2482        #Check first value
2483        elevation = fid.variables['elevation'][:]
2484        stage = fid.variables['stage'][:]
2485        xmomentum = fid.variables['xmomentum'][:]
2486        ymomentum = fid.variables['ymomentum'][:]
2487
2488        #print ymomentum
2489        first_height = first_amp/100 - first_elevation
2490        third_height = third_amp/100 - third_elevation
2491        first_momentum=first_speed*first_height/100
2492        third_momentum=third_speed*third_height/100
2493
2494        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2495        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2496
2497        fid.close()
2498
2499        #Cleanup
2500        os.remove('test.sww')
2501
2502
2503
2504
2505    def test_ferret2sww_nz_origin(self):
2506        from Scientific.IO.NetCDF import NetCDFFile
2507        from coordinate_transforms.redfearn import redfearn
2508
2509        #Call conversion (with nonzero origin)
2510        ferret2sww(self.test_MOST_file, verbose=False,
2511                   origin = (56, 100000, 200000))
2512
2513
2514        #Work out the UTM coordinates for first point
2515        zone, e, n = redfearn(-34.5, 150.66667)
2516
2517        #Read output file 'small.sww'
2518        #fid = NetCDFFile('small.sww', 'r')
2519        fid = NetCDFFile(self.test_MOST_file + '.sww')
2520       
2521        x = fid.variables['x'][:]
2522        y = fid.variables['y'][:]
2523
2524        #Check that first coordinate is correctly represented
2525        assert allclose(x[0], e-100000)
2526        assert allclose(y[0], n-200000)
2527
2528        fid.close()
2529
2530        #Cleanup
2531        os.remove(self.test_MOST_file + '.sww')
2532
2533
2534
2535    def test_sww_extent(self):
2536        """Not a test, rather a look at the sww format
2537        """
2538
2539        import time, os
2540        from Numeric import array, zeros, allclose, Float, concatenate
2541        from Scientific.IO.NetCDF import NetCDFFile
2542
2543        self.domain.filename = 'datatest' + str(id(self))
2544        self.domain.format = 'sww'
2545        self.domain.smooth = True
2546        self.domain.reduction = mean
2547        self.domain.set_datadir('.')
2548
2549
2550        sww = get_dataobject(self.domain)
2551        sww.store_connectivity()
2552        sww.store_timestep('stage')
2553        self.domain.time = 2.
2554
2555        #Modify stage at second timestep
2556        stage = self.domain.quantities['stage'].vertex_values
2557        self.domain.set_quantity('stage', stage/2)
2558
2559        sww.store_timestep('stage')
2560
2561        file_and_extension_name = self.domain.filename + ".sww"
2562        #print "file_and_extension_name",file_and_extension_name
2563        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2564               extent_sww(file_and_extension_name )
2565
2566        assert allclose(xmin, 0.0)
2567        assert allclose(xmax, 1.0)
2568        assert allclose(ymin, 0.0)
2569        assert allclose(ymax, 1.0)
2570        assert allclose(stagemin, -0.85)
2571        assert allclose(stagemax, 0.15)
2572
2573
2574        #Cleanup
2575        os.remove(sww.filename)
2576
2577
2578
2579    def test_sww2domain1(self):
2580        ################################################
2581        #Create a test domain, and evolve and save it.
2582        ################################################
2583        from mesh_factory import rectangular
2584        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2585             Constant_height, Time_boundary, Transmissive_boundary
2586        from Numeric import array
2587
2588        #Create basic mesh
2589
2590        yiel=0.01
2591        points, vertices, boundary = rectangular(10,10)
2592
2593        #Create shallow water domain
2594        domain = Domain(points, vertices, boundary)
2595        domain.geo_reference = Geo_reference(56,11,11)
2596        domain.smooth = False
2597        domain.visualise = False
2598        domain.store = True
2599        domain.filename = 'bedslope'
2600        domain.default_order=2
2601        #Bed-slope and friction
2602        domain.set_quantity('elevation', lambda x,y: -x/3)
2603        domain.set_quantity('friction', 0.1)
2604        # Boundary conditions
2605        from math import sin, pi
2606        Br = Reflective_boundary(domain)
2607        Bt = Transmissive_boundary(domain)
2608        Bd = Dirichlet_boundary([0.2,0.,0.])
2609        Bw = Time_boundary(domain=domain,
2610                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2611
2612        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2613        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2614
2615        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2616        #Initial condition
2617        h = 0.05
2618        elevation = domain.quantities['elevation'].vertex_values
2619        domain.set_quantity('stage', elevation + h)
2620
2621        domain.check_integrity()
2622        #Evolution
2623        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2624        #    domain.write_time()
2625            pass
2626
2627
2628        ##########################################
2629        #Import the example's file as a new domain
2630        ##########################################
2631        from data_manager import sww2domain
2632        from Numeric import allclose
2633        import os
2634
2635        filename = domain.datadir+os.sep+domain.filename+'.sww'
2636        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2637        #points, vertices, boundary = rectangular(15,15)
2638        #domain2.boundary = boundary
2639        ###################
2640        ##NOW TEST IT!!!
2641        ###################
2642
2643        #os.remove(domain.filename + '.sww')
2644        os.remove(filename)       
2645
2646        bits = ['vertex_coordinates']
2647        for quantity in ['elevation']+domain.quantities_to_be_stored:
2648            bits.append('get_quantity("%s").get_integral()' %quantity)
2649            bits.append('get_quantity("%s").get_values()' %quantity)
2650
2651        for bit in bits:
2652            #print 'testing that domain.'+bit+' has been restored'
2653            #print bit
2654        #print 'done'
2655            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2656
2657        ######################################
2658        #Now evolve them both, just to be sure
2659        ######################################x
2660        visualise = False
2661        #visualise = True
2662        domain.visualise = visualise
2663        domain.time = 0.
2664        from time import sleep
2665
2666        final = .1
2667        domain.set_quantity('friction', 0.1)
2668        domain.store = False
2669        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2670
2671
2672        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2673            if visualise: sleep(1.)
2674            #domain.write_time()
2675            pass
2676
2677        final = final - (domain2.starttime-domain.starttime)
2678        #BUT since domain1 gets time hacked back to 0:
2679        final = final + (domain2.starttime-domain.starttime)
2680
2681        domain2.smooth = False
2682        domain2.visualise = visualise
2683        domain2.store = False
2684        domain2.default_order=2
2685        domain2.set_quantity('friction', 0.1)
2686        #Bed-slope and friction
2687        # Boundary conditions
2688        Bd2=Dirichlet_boundary([0.2,0.,0.])
2689        domain2.boundary = domain.boundary
2690        #print 'domain2.boundary'
2691        #print domain2.boundary
2692        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2693        #domain2.set_boundary({'exterior': Bd})
2694
2695        domain2.check_integrity()
2696
2697        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2698            if visualise: sleep(1.)
2699            #domain2.write_time()
2700            pass
2701
2702        ###################
2703        ##NOW TEST IT!!!
2704        ##################
2705
2706        bits = [ 'vertex_coordinates']
2707
2708        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2709            bits.append('get_quantity("%s").get_integral()' %quantity)
2710            bits.append('get_quantity("%s").get_values()' %quantity)
2711
2712        for bit in bits:
2713            #print bit
2714            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2715
2716
2717    def test_sww2domain2(self):
2718        ##################################################################
2719        #Same as previous test, but this checks how NaNs are handled.
2720        ##################################################################
2721
2722
2723        from mesh_factory import rectangular
2724        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2725             Constant_height, Time_boundary, Transmissive_boundary
2726        from Numeric import array
2727
2728        #Create basic mesh
2729        points, vertices, boundary = rectangular(2,2)
2730
2731        #Create shallow water domain
2732        domain = Domain(points, vertices, boundary)
2733        domain.smooth = False
2734        domain.visualise = False
2735        domain.store = True
2736        domain.set_name('test_file')
2737        domain.set_datadir('.')
2738        domain.default_order=2
2739        domain.quantities_to_be_stored=['stage']
2740
2741        domain.set_quantity('elevation', lambda x,y: -x/3)
2742        domain.set_quantity('friction', 0.1)
2743
2744        from math import sin, pi
2745        Br = Reflective_boundary(domain)
2746        Bt = Transmissive_boundary(domain)
2747        Bd = Dirichlet_boundary([0.2,0.,0.])
2748        Bw = Time_boundary(domain=domain,
2749                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2750
2751        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2752
2753        h = 0.05
2754        elevation = domain.quantities['elevation'].vertex_values
2755        domain.set_quantity('stage', elevation + h)
2756
2757        domain.check_integrity()
2758
2759        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2760            pass
2761            #domain.write_time()
2762
2763
2764
2765        ##################################
2766        #Import the file as a new domain
2767        ##################################
2768        from data_manager import sww2domain
2769        from Numeric import allclose
2770        import os
2771
2772        filename = domain.datadir+os.sep+domain.filename+'.sww'
2773
2774        #Fail because NaNs are present
2775        try:
2776            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2777        except:
2778            #Now import it, filling NaNs to be 0
2779            filler = 0
2780            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2781
2782        #Clean up
2783        os.remove(filename)
2784           
2785
2786        bits = [ 'geo_reference.get_xllcorner()',
2787                'geo_reference.get_yllcorner()',
2788                'vertex_coordinates']
2789
2790        for quantity in ['elevation']+domain.quantities_to_be_stored:
2791            bits.append('get_quantity("%s").get_integral()' %quantity)
2792            bits.append('get_quantity("%s").get_values()' %quantity)
2793
2794        for bit in bits:
2795        #    print 'testing that domain.'+bit+' has been restored'
2796            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2797
2798        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
2799        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
2800        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
2801        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
2802
2803
2804
2805    #def test_weed(self):
2806        from data_manager import weed
2807
2808        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2809        volumes1 = [[0,1,2],[3,4,5]]
2810        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2811        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2812
2813        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2814
2815        assert len(points2)==len(coordinates2)
2816        for i in range(len(coordinates2)):
2817            coordinate = tuple(coordinates2[i])
2818            assert points2.has_key(coordinate)
2819            points2[coordinate]=i
2820
2821        for triangle in volumes1:
2822            for coordinate in triangle:
2823                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2824                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2825
2826
2827    #FIXME This fails - smooth makes the comparism too hard for allclose
2828    def ztest_sww2domain3(self):
2829        ################################################
2830        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2831        ################################################
2832        from mesh_factory import rectangular
2833        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2834             Constant_height, Time_boundary, Transmissive_boundary
2835        from Numeric import array
2836        #Create basic mesh
2837
2838        yiel=0.01
2839        points, vertices, boundary = rectangular(10,10)
2840
2841        #Create shallow water domain
2842        domain = Domain(points, vertices, boundary)
2843        domain.geo_reference = Geo_reference(56,11,11)
2844        domain.smooth = True
2845        domain.visualise = False
2846        domain.store = True
2847        domain.filename = 'bedslope'
2848        domain.default_order=2
2849        #Bed-slope and friction
2850        domain.set_quantity('elevation', lambda x,y: -x/3)
2851        domain.set_quantity('friction', 0.1)
2852        # Boundary conditions
2853        from math import sin, pi
2854        Br = Reflective_boundary(domain)
2855        Bt = Transmissive_boundary(domain)
2856        Bd = Dirichlet_boundary([0.2,0.,0.])
2857        Bw = Time_boundary(domain=domain,
2858                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2859
2860        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2861
2862        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2863        #Initial condition
2864        h = 0.05
2865        elevation = domain.quantities['elevation'].vertex_values
2866        domain.set_quantity('stage', elevation + h)
2867
2868
2869        domain.check_integrity()
2870        #Evolution
2871        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2872        #    domain.write_time()
2873            pass
2874
2875
2876        ##########################################
2877        #Import the example's file as a new domain
2878        ##########################################
2879        from data_manager import sww2domain
2880        from Numeric import allclose
2881        import os
2882
2883        filename = domain.datadir+os.sep+domain.filename+'.sww'
2884        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2885        #points, vertices, boundary = rectangular(15,15)
2886        #domain2.boundary = boundary
2887        ###################
2888        ##NOW TEST IT!!!
2889        ###################
2890
2891        os.remove(domain.filename + '.sww')
2892
2893        #FIXME smooth domain so that they can be compared
2894
2895
2896        bits = []#'vertex_coordinates']
2897        for quantity in ['elevation']+domain.quantities_to_be_stored:
2898            bits.append('quantities["%s"].get_integral()'%quantity)
2899
2900
2901        for bit in bits:
2902            #print 'testing that domain.'+bit+' has been restored'
2903            #print bit
2904            #print 'done'
2905            #print ('domain.'+bit), eval('domain.'+bit)
2906            #print ('domain2.'+bit), eval('domain2.'+bit)
2907            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2908            pass
2909
2910        ######################################
2911        #Now evolve them both, just to be sure
2912        ######################################x
2913        visualise = False
2914        visualise = True
2915        domain.visualise = visualise
2916        domain.time = 0.
2917        from time import sleep
2918
2919        final = .5
2920        domain.set_quantity('friction', 0.1)
2921        domain.store = False
2922        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2923
2924        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2925            if visualise: sleep(.03)
2926            #domain.write_time()
2927            pass
2928
2929        domain2.smooth = True
2930        domain2.visualise = visualise
2931        domain2.store = False
2932        domain2.default_order=2
2933        domain2.set_quantity('friction', 0.1)
2934        #Bed-slope and friction
2935        # Boundary conditions
2936        Bd2=Dirichlet_boundary([0.2,0.,0.])
2937        Br2 = Reflective_boundary(domain2)
2938        domain2.boundary = domain.boundary
2939        #print 'domain2.boundary'
2940        #print domain2.boundary
2941        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2942        #domain2.boundary = domain.boundary
2943        #domain2.set_boundary({'exterior': Bd})
2944
2945        domain2.check_integrity()
2946
2947        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2948            if visualise: sleep(.03)
2949            #domain2.write_time()
2950            pass
2951
2952        ###################
2953        ##NOW TEST IT!!!
2954        ##################
2955
2956        print '><><><><>>'
2957        bits = [ 'vertex_coordinates']
2958
2959        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2960            #bits.append('quantities["%s"].get_integral()'%quantity)
2961            bits.append('get_quantity("%s").get_values()' %quantity)
2962
2963        for bit in bits:
2964            print bit
2965            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2966
2967
2968    def test_decimate_dem(self):
2969        """Test decimation of dem file
2970        """
2971
2972        import os
2973        from Numeric import ones, allclose, Float, arange
2974        from Scientific.IO.NetCDF import NetCDFFile
2975
2976        #Write test dem file
2977        root = 'decdemtest'
2978
2979        filename = root + '.dem'
2980        fid = NetCDFFile(filename, 'w')
2981
2982        fid.institution = 'Geoscience Australia'
2983        fid.description = 'NetCDF DEM format for compact and portable ' +\
2984                          'storage of spatial point data'
2985
2986        nrows = 15
2987        ncols = 18
2988
2989        fid.ncols = ncols
2990        fid.nrows = nrows
2991        fid.xllcorner = 2000.5
2992        fid.yllcorner = 3000.5
2993        fid.cellsize = 25
2994        fid.NODATA_value = -9999
2995
2996        fid.zone = 56
2997        fid.false_easting = 0.0
2998        fid.false_northing = 0.0
2999        fid.projection = 'UTM'
3000        fid.datum = 'WGS84'
3001        fid.units = 'METERS'
3002
3003        fid.createDimension('number_of_points', nrows*ncols)
3004
3005        fid.createVariable('elevation', Float, ('number_of_points',))
3006
3007        elevation = fid.variables['elevation']
3008
3009        elevation[:] = (arange(nrows*ncols))
3010
3011        fid.close()
3012
3013        #generate the elevation values expected in the decimated file
3014        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
3015                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3016                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3017                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3018                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3019                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
3020                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3021                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3022                         (144+145+146+162+163+164+180+181+182) / 9.0,
3023                         (148+149+150+166+167+168+184+185+186) / 9.0,
3024                         (152+153+154+170+171+172+188+189+190) / 9.0,
3025                         (156+157+158+174+175+176+192+193+194) / 9.0,
3026                         (216+217+218+234+235+236+252+253+254) / 9.0,
3027                         (220+221+222+238+239+240+256+257+258) / 9.0,
3028                         (224+225+226+242+243+244+260+261+262) / 9.0,
3029                         (228+229+230+246+247+248+264+265+266) / 9.0]
3030
3031        #generate a stencil for computing the decimated values
3032        stencil = ones((3,3), Float) / 9.0
3033
3034        decimate_dem(root, stencil=stencil, cellsize_new=100)
3035
3036        #Open decimated NetCDF file
3037        fid = NetCDFFile(root + '_100.dem', 'r')
3038
3039        # Get decimated elevation
3040        elevation = fid.variables['elevation']
3041
3042        #Check values
3043        assert allclose(elevation, ref_elevation)
3044
3045        #Cleanup
3046        fid.close()
3047
3048        os.remove(root + '.dem')
3049        os.remove(root + '_100.dem')
3050
3051    def test_decimate_dem_NODATA(self):
3052        """Test decimation of dem file that includes NODATA values
3053        """
3054
3055        import os
3056        from Numeric import ones, allclose, Float, arange, reshape
3057        from Scientific.IO.NetCDF import NetCDFFile
3058
3059        #Write test dem file
3060        root = 'decdemtest'
3061
3062        filename = root + '.dem'
3063        fid = NetCDFFile(filename, 'w')
3064
3065        fid.institution = 'Geoscience Australia'
3066        fid.description = 'NetCDF DEM format for compact and portable ' +\
3067                          'storage of spatial point data'
3068
3069        nrows = 15
3070        ncols = 18
3071        NODATA_value = -9999
3072
3073        fid.ncols = ncols
3074        fid.nrows = nrows
3075        fid.xllcorner = 2000.5
3076        fid.yllcorner = 3000.5
3077        fid.cellsize = 25
3078        fid.NODATA_value = NODATA_value
3079
3080        fid.zone = 56
3081        fid.false_easting = 0.0
3082        fid.false_northing = 0.0
3083        fid.projection = 'UTM'
3084        fid.datum = 'WGS84'
3085        fid.units = 'METERS'
3086
3087        fid.createDimension('number_of_points', nrows*ncols)
3088
3089        fid.createVariable('elevation', Float, ('number_of_points',))
3090
3091        elevation = fid.variables['elevation']
3092
3093        #generate initial elevation values
3094        elevation_tmp = (arange(nrows*ncols))
3095        #add some NODATA values
3096        elevation_tmp[0]   = NODATA_value
3097        elevation_tmp[95]  = NODATA_value
3098        elevation_tmp[188] = NODATA_value
3099        elevation_tmp[189] = NODATA_value
3100        elevation_tmp[190] = NODATA_value
3101        elevation_tmp[209] = NODATA_value
3102        elevation_tmp[252] = NODATA_value
3103
3104        elevation[:] = elevation_tmp
3105
3106        fid.close()
3107
3108        #generate the elevation values expected in the decimated file
3109        ref_elevation = [NODATA_value,
3110                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3111                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3112                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3113                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3114                         NODATA_value,
3115                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3116                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3117                         (144+145+146+162+163+164+180+181+182) / 9.0,
3118                         (148+149+150+166+167+168+184+185+186) / 9.0,
3119                         NODATA_value,
3120                         (156+157+158+174+175+176+192+193+194) / 9.0,
3121                         NODATA_value,
3122                         (220+221+222+238+239+240+256+257+258) / 9.0,
3123                         (224+225+226+242+243+244+260+261+262) / 9.0,
3124                         (228+229+230+246+247+248+264+265+266) / 9.0]
3125
3126        #generate a stencil for computing the decimated values
3127        stencil = ones((3,3), Float) / 9.0
3128
3129        decimate_dem(root, stencil=stencil, cellsize_new=100)
3130
3131        #Open decimated NetCDF file
3132        fid = NetCDFFile(root + '_100.dem', 'r')
3133
3134        # Get decimated elevation
3135        elevation = fid.variables['elevation']
3136
3137        #Check values
3138        assert allclose(elevation, ref_elevation)
3139
3140        #Cleanup
3141        fid.close()
3142
3143        os.remove(root + '.dem')
3144        os.remove(root + '_100.dem')
3145
3146    def xxxtestz_sww2ers_real(self):
3147        """Test that sww information can be converted correctly to asc/prj
3148        format readable by e.g. ArcView
3149        """
3150
3151        import time, os
3152        from Numeric import array, zeros, allclose, Float, concatenate
3153        from Scientific.IO.NetCDF import NetCDFFile
3154
3155        # the memory optimised least squares
3156        #  cellsize = 20,   # this one seems to hang
3157        #  cellsize = 200000, # Ran 1 test in 269.703s
3158                                #Ran 1 test in 267.344s
3159        #  cellsize = 20000,  # Ran 1 test in 460.922s
3160        #  cellsize = 2000   #Ran 1 test in 5340.250s
3161        #  cellsize = 200   #this one seems to hang, building matirx A
3162
3163        # not optimised
3164        # seems to hang
3165        #  cellsize = 2000   # Ran 1 test in 5334.563s
3166        #Export to ascii/prj files
3167        sww2dem('karratha_100m',
3168                quantity = 'depth',
3169                cellsize = 200000,
3170                verbose = True)
3171
3172    def test_read_asc(self):
3173        """Test conversion from dem in ascii format to native NetCDF xya format
3174        """
3175
3176        import time, os
3177        from Numeric import array, zeros, allclose, Float, concatenate
3178        from Scientific.IO.NetCDF import NetCDFFile
3179
3180        import data_manager
3181        #Write test asc file
3182        filename = tempfile.mktemp(".000")
3183        fid = open(filename, 'w')
3184        fid.write("""ncols         7
3185nrows         4
3186xllcorner     2000.5
3187yllcorner     3000.5
3188cellsize      25
3189NODATA_value  -9999
3190    97.921    99.285   125.588   180.830   258.645   342.872   415.836
3191   473.157   514.391   553.893   607.120   678.125   777.283   883.038
3192   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
3193   502.645   516.230   504.739   450.604   388.500   338.097   514.980
3194""")
3195        fid.close()
3196        bath_metadata, grid = \
3197                   data_manager._read_asc(filename, verbose=False)
3198        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
3199        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
3200        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
3201        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
3202        self.failUnless(grid[0][0]  == 97.921,  'Failed')
3203        self.failUnless(grid[3][6]  == 514.980,  'Failed')
3204
3205        os.remove(filename)
3206       
3207    def test_asc_csiro2sww(self):
3208        import tempfile
3209
3210        bath_dir = tempfile.mkdtemp()
3211        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3212        #bath_dir = 'bath_data_manager_test'
3213        #print "os.getcwd( )",os.getcwd( )
3214        elevation_dir =  tempfile.mkdtemp()
3215        #elevation_dir = 'elev_expanded'
3216        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3217        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3218       
3219        fid = open(bath_dir_filename, 'w')
3220        fid.write(""" ncols             3
3221 nrows             2
3222 xllcorner    148.00000
3223 yllcorner    -38.00000
3224 cellsize       0.25
3225 nodata_value   -9999.0
3226    9000.000    -1000.000    3000.0
3227   -1000.000    9000.000  -1000.000
3228""")
3229        fid.close()
3230       
3231        fid = open(elevation_dir_filename1, 'w')
3232        fid.write(""" ncols             3
3233 nrows             2
3234 xllcorner    148.00000
3235 yllcorner    -38.00000
3236 cellsize       0.25
3237 nodata_value   -9999.0
3238    9000.000    0.000    3000.0
3239     0.000     9000.000     0.000
3240""")
3241        fid.close()
3242
3243        fid = open(elevation_dir_filename2, 'w')
3244        fid.write(""" ncols             3
3245 nrows             2
3246 xllcorner    148.00000
3247 yllcorner    -38.00000
3248 cellsize       0.25
3249 nodata_value   -9999.0
3250    9000.000    4000.000    4000.0
3251    4000.000    9000.000    4000.000
3252""")
3253        fid.close()
3254         
3255        ucur_dir =  tempfile.mkdtemp()
3256        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3257        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3258       
3259        fid = open(ucur_dir_filename1, 'w')
3260        fid.write(""" ncols             3
3261 nrows             2
3262 xllcorner    148.00000
3263 yllcorner    -38.00000
3264 cellsize       0.25
3265 nodata_value   -9999.0
3266    90.000    60.000    30.0
3267    10.000    10.000    10.000
3268""")
3269        fid.close()
3270        fid = open(ucur_dir_filename2, 'w')
3271        fid.write(""" ncols             3
3272 nrows             2
3273 xllcorner    148.00000
3274 yllcorner    -38.00000
3275 cellsize       0.25
3276 nodata_value   -9999.0
3277    90.000    60.000    30.0
3278    10.000    10.000    10.000
3279""")
3280        fid.close()
3281       
3282        vcur_dir =  tempfile.mkdtemp()
3283        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3284        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3285       
3286        fid = open(vcur_dir_filename1, 'w')
3287        fid.write(""" ncols             3
3288 nrows             2
3289 xllcorner    148.00000
3290 yllcorner    -38.00000
3291 cellsize       0.25
3292 nodata_value   -9999.0
3293    90.000    60.000    30.0
3294    10.000    10.000    10.000
3295""")
3296        fid.close()
3297        fid = open(vcur_dir_filename2, 'w')
3298        fid.write(""" ncols             3
3299 nrows             2
3300 xllcorner    148.00000
3301 yllcorner    -38.00000
3302 cellsize       0.25
3303 nodata_value   -9999.0
3304    90.000    60.000    30.0
3305    10.000    10.000    10.000
3306""")
3307        fid.close()
3308       
3309        sww_file = 'a_test.sww'
3310        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
3311
3312        # check the sww file
3313       
3314        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3315        x = fid.variables['x'][:]
3316        y = fid.variables['y'][:]
3317        z = fid.variables['z'][:]
3318        stage = fid.variables['stage'][:]
3319        xmomentum = fid.variables['xmomentum'][:]
3320        geo_ref = Geo_reference(NetCDFObject=fid)
3321        #print "geo_ref",geo_ref
3322        x_ref = geo_ref.get_xllcorner()
3323        y_ref = geo_ref.get_yllcorner()
3324        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3325        assert allclose(x_ref, 587798.418) # (-38, 148)
3326        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3327       
3328        #Zone:   55   
3329        #Easting:  588095.674  Northing: 5821451.722
3330        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3331        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3332
3333        #Zone:   55   
3334        #Easting:  632145.632  Northing: 5820863.269
3335        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3336        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3337
3338        #Zone:   55   
3339        #Easting:  609748.788  Northing: 5793447.860
3340        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3341        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3342
3343        assert allclose(z[0],9000.0 )
3344        assert allclose(stage[0][1],0.0 )
3345
3346        #(4000+1000)*60
3347        assert allclose(xmomentum[1][1],300000.0 )
3348       
3349       
3350        fid.close()
3351   
3352        #tidy up
3353        os.remove(bath_dir_filename)
3354        os.rmdir(bath_dir)
3355
3356        os.remove(elevation_dir_filename1)
3357        os.remove(elevation_dir_filename2)
3358        os.rmdir(elevation_dir)
3359       
3360        os.remove(ucur_dir_filename1)
3361        os.remove(ucur_dir_filename2)
3362        os.rmdir(ucur_dir)
3363       
3364        os.remove(vcur_dir_filename1)
3365        os.remove(vcur_dir_filename2)
3366        os.rmdir(vcur_dir)
3367
3368
3369        # remove sww file
3370        os.remove(sww_file)
3371       
3372    def test_asc_csiro2sww2(self):
3373        import tempfile
3374
3375        bath_dir = tempfile.mkdtemp()
3376        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3377        #bath_dir = 'bath_data_manager_test'
3378        #print "os.getcwd( )",os.getcwd( )
3379        elevation_dir =  tempfile.mkdtemp()
3380        #elevation_dir = 'elev_expanded'
3381        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3382        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3383       
3384        fid = open(bath_dir_filename, 'w')
3385        fid.write(""" ncols             3
3386 nrows             2
3387 xllcorner    148.00000
3388 yllcorner    -38.00000
3389 cellsize       0.25
3390 nodata_value   -9999.0
3391    9000.000    -1000.000    3000.0
3392   -1000.000    9000.000  -1000.000
3393""")
3394        fid.close()
3395       
3396        fid = open(elevation_dir_filename1, 'w')
3397        fid.write(""" ncols             3
3398 nrows             2
3399 xllcorner    148.00000
3400 yllcorner    -38.00000
3401 cellsize       0.25
3402 nodata_value   -9999.0
3403    9000.000    0.000    3000.0
3404     0.000     -9999.000     -9999.000
3405""")
3406        fid.close()
3407
3408        fid = open(elevation_dir_filename2, 'w')
3409        fid.write(""" ncols             3
3410 nrows             2
3411 xllcorner    148.00000
3412 yllcorner    -38.00000
3413 cellsize       0.25
3414 nodata_value   -9999.0
3415    9000.000    4000.000    4000.0
3416    4000.000    9000.000    4000.000
3417""")
3418        fid.close()
3419         
3420        ucur_dir =  tempfile.mkdtemp()
3421        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3422        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3423       
3424        fid = open(ucur_dir_filename1, 'w')
3425        fid.write(""" ncols             3
3426 nrows             2
3427 xllcorner    148.00000
3428 yllcorner    -38.00000
3429 cellsize       0.25
3430 nodata_value   -9999.0
3431    90.000    60.000    30.0
3432    10.000    10.000    10.000
3433""")
3434        fid.close()
3435        fid = open(ucur_dir_filename2, 'w')
3436        fid.write(""" ncols             3
3437 nrows             2
3438 xllcorner    148.00000
3439 yllcorner    -38.00000
3440 cellsize       0.25
3441 nodata_value   -9999.0
3442    90.000    60.000    30.0
3443    10.000    10.000    10.000
3444""")
3445        fid.close()
3446       
3447        vcur_dir =  tempfile.mkdtemp()
3448        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3449        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3450       
3451        fid = open(vcur_dir_filename1, 'w')
3452        fid.write(""" ncols             3
3453 nrows             2
3454 xllcorner    148.00000
3455 yllcorner    -38.00000
3456 cellsize       0.25
3457 nodata_value   -9999.0
3458    90.000    60.000    30.0
3459    10.000    10.000    10.000
3460""")
3461        fid.close()
3462        fid = open(vcur_dir_filename2, 'w')
3463        fid.write(""" ncols             3
3464 nrows             2
3465 xllcorner    148.00000
3466 yllcorner    -38.00000
3467 cellsize       0.25
3468 nodata_value   -9999.0
3469    90.000    60.000    30.0
3470    10.000    10.000    10.000
3471""")
3472        fid.close()
3473       
3474        try:
3475            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3476                          vcur_dir, sww_file)
3477        except:
3478            #tidy up
3479            os.remove(bath_dir_filename)
3480            os.rmdir(bath_dir)
3481
3482            os.remove(elevation_dir_filename1)
3483            os.remove(elevation_dir_filename2)
3484            os.rmdir(elevation_dir)
3485       
3486            os.remove(ucur_dir_filename1)
3487            os.remove(ucur_dir_filename2)
3488            os.rmdir(ucur_dir)
3489       
3490            os.remove(vcur_dir_filename1)
3491            os.remove(vcur_dir_filename2)
3492            os.rmdir(vcur_dir)
3493        else:
3494            #tidy up
3495            os.remove(bath_dir_filename)
3496            os.rmdir(bath_dir)
3497
3498            os.remove(elevation_dir_filename1)
3499            os.remove(elevation_dir_filename2)
3500            os.rmdir(elevation_dir)
3501            raise 'Should raise exception'
3502       
3503            os.remove(ucur_dir_filename1)
3504            os.remove(ucur_dir_filename2)
3505            os.rmdir(ucur_dir)
3506       
3507            os.remove(vcur_dir_filename1)
3508            os.remove(vcur_dir_filename2)
3509            os.rmdir(vcur_dir)
3510
3511       
3512     
3513    def test_asc_csiro2sww3(self):
3514        import tempfile
3515
3516        bath_dir = tempfile.mkdtemp()
3517        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3518        #bath_dir = 'bath_data_manager_test'
3519        #print "os.getcwd( )",os.getcwd( )
3520        elevation_dir =  tempfile.mkdtemp()
3521        #elevation_dir = 'elev_expanded'
3522        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3523        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3524       
3525        fid = open(bath_dir_filename, 'w')
3526        fid.write(""" ncols             3
3527 nrows             2
3528 xllcorner    148.00000
3529 yllcorner    -38.00000
3530 cellsize       0.25
3531 nodata_value   -9999.0
3532    9000.000    -1000.000    3000.0
3533   -1000.000    9000.000  -1000.000
3534""")
3535        fid.close()
3536       
3537        fid = open(elevation_dir_filename1, 'w')
3538        fid.write(""" ncols             3
3539 nrows             2
3540 xllcorner    148.00000
3541 yllcorner    -38.00000
3542 cellsize       0.25
3543 nodata_value   -9999.0
3544    9000.000    0.000    3000.0
3545     0.000     -9999.000     -9999.000
3546""")
3547        fid.close()
3548
3549        fid = open(elevation_dir_filename2, 'w')
3550        fid.write(""" ncols             3
3551 nrows             2
3552 xllcorner    148.00000
3553 yllcorner    -38.00000
3554 cellsize       0.25
3555 nodata_value   -9999.0
3556    9000.000    4000.000    4000.0
3557    4000.000    9000.000    4000.000
3558""")
3559        fid.close()
3560         
3561        ucur_dir =  tempfile.mkdtemp()
3562        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3563        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3564       
3565        fid = open(ucur_dir_filename1, 'w')
3566        fid.write(""" ncols             3
3567 nrows             2
3568 xllcorner    148.00000
3569 yllcorner    -38.00000
3570 cellsize       0.25
3571 nodata_value   -9999.0
3572    90.000    60.000    30.0
3573    10.000    10.000    10.000
3574""")
3575        fid.close()
3576        fid = open(ucur_dir_filename2, 'w')
3577        fid.write(""" ncols             3
3578 nrows             2
3579 xllcorner    148.00000
3580 yllcorner    -38.00000
3581 cellsize       0.25
3582 nodata_value   -9999.0
3583    90.000    60.000    30.0
3584    10.000    10.000    10.000
3585""")
3586        fid.close()
3587       
3588        vcur_dir =  tempfile.mkdtemp()
3589        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3590        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3591       
3592        fid = open(vcur_dir_filename1, 'w')
3593        fid.write(""" ncols             3
3594 nrows             2
3595 xllcorner    148.00000
3596 yllcorner    -38.00000
3597 cellsize       0.25
3598 nodata_value   -9999.0
3599    90.000    60.000    30.0
3600    10.000    10.000    10.000
3601""")
3602        fid.close()
3603        fid = open(vcur_dir_filename2, 'w')
3604        fid.write(""" ncols             3
3605 nrows             2
3606 xllcorner    148.00000
3607 yllcorner    -38.00000
3608 cellsize       0.25
3609 nodata_value   -9999.0
3610    90.000    60.000    30.0
3611    10.000    10.000    10.000
3612""")
3613        fid.close()
3614       
3615        sww_file = 'a_test.sww'
3616        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3617                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3618                      mean_stage = 100)
3619
3620        # check the sww file
3621       
3622        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3623        x = fid.variables['x'][:]
3624        y = fid.variables['y'][:]
3625        z = fid.variables['z'][:]
3626        stage = fid.variables['stage'][:]
3627        xmomentum = fid.variables['xmomentum'][:]
3628        geo_ref = Geo_reference(NetCDFObject=fid)
3629        #print "geo_ref",geo_ref
3630        x_ref = geo_ref.get_xllcorner()
3631        y_ref = geo_ref.get_yllcorner()
3632        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3633        assert allclose(x_ref, 587798.418) # (-38, 148)
3634        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3635       
3636        #Zone:   55   
3637        #Easting:  588095.674  Northing: 5821451.722
3638        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3639        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3640
3641        #Zone:   55   
3642        #Easting:  632145.632  Northing: 5820863.269
3643        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3644        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3645
3646        #Zone:   55   
3647        #Easting:  609748.788  Northing: 5793447.860
3648        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3649        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3650
3651        assert allclose(z[0],9000.0 )
3652        assert allclose(stage[0][4],100.0 )
3653        assert allclose(stage[0][5],100.0 )
3654
3655        #(100.0 - 9000)*10
3656        assert allclose(xmomentum[0][4], -89000.0 )
3657       
3658        #(100.0 - -1000.000)*10
3659        assert allclose(xmomentum[0][5], 11000.0 )
3660       
3661        fid.close()
3662   
3663        #tidy up
3664        os.remove(bath_dir_filename)
3665        os.rmdir(bath_dir)
3666
3667        os.remove(elevation_dir_filename1)
3668        os.remove(elevation_dir_filename2)
3669        os.rmdir(elevation_dir)
3670       
3671        os.remove(ucur_dir_filename1)
3672        os.remove(ucur_dir_filename2)
3673        os.rmdir(ucur_dir)
3674       
3675        os.remove(vcur_dir_filename1)
3676        os.remove(vcur_dir_filename2)
3677        os.rmdir(vcur_dir)
3678
3679        # remove sww file
3680        os.remove(sww_file)
3681   
3682     
3683    def test_asc_csiro2sww4(self):
3684        """
3685        Test specifying the extent
3686        """
3687       
3688        import tempfile
3689
3690        bath_dir = tempfile.mkdtemp()
3691        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3692        #bath_dir = 'bath_data_manager_test'
3693        #print "os.getcwd( )",os.getcwd( )
3694        elevation_dir =  tempfile.mkdtemp()
3695        #elevation_dir = 'elev_expanded'
3696        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3697        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3698       
3699        fid = open(bath_dir_filename, 'w')
3700        fid.write(""" ncols             4
3701 nrows             4
3702 xllcorner    148.00000
3703 yllcorner    -38.00000
3704 cellsize       0.25
3705 nodata_value   -9999.0
3706   -9000.000    -1000.000   -3000.0 -2000.000
3707   -1000.000    9000.000  -1000.000 -3000.000
3708   -4000.000    6000.000   2000.000 -5000.000
3709   -9000.000    -1000.000   -3000.0 -2000.000
3710""")
3711        fid.close()
3712       
3713        fid = open(elevation_dir_filename1, 'w')
3714        fid.write(""" ncols             4
3715 nrows             4
3716 xllcorner    148.00000
3717 yllcorner    -38.00000
3718 cellsize       0.25
3719 nodata_value   -9999.0
3720   -900.000    -100.000   -300.0 -200.000
3721   -100.000    900.000  -100.000 -300.000
3722   -400.000    600.000   200.000 -500.000
3723   -900.000    -100.000   -300.0 -200.000
3724""")
3725        fid.close()
3726
3727        fid = open(elevation_dir_filename2, 'w')
3728        fid.write(""" ncols             4
3729 nrows             4
3730 xllcorner    148.00000
3731 yllcorner    -38.00000
3732 cellsize       0.25
3733 nodata_value   -9999.0
3734   -990.000    -110.000   -330.0 -220.000
3735   -110.000    990.000  -110.000 -330.000
3736   -440.000    660.000   220.000 -550.000
3737   -990.000    -110.000   -330.0 -220.000
3738""")
3739        fid.close()
3740         
3741        ucur_dir =  tempfile.mkdtemp()
3742        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3743        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3744       
3745        fid = open(ucur_dir_filename1, 'w')
3746        fid.write(""" ncols             4
3747 nrows             4
3748 xllcorner    148.00000
3749 yllcorner    -38.00000
3750 cellsize       0.25
3751 nodata_value   -9999.0
3752   -90.000    -10.000   -30.0 -20.000
3753   -10.000    90.000  -10.000 -30.000
3754   -40.000    60.000   20.000 -50.000
3755   -90.000    -10.000   -30.0 -20.000
3756""")
3757        fid.close()
3758        fid = open(ucur_dir_filename2, 'w')
3759        fid.write(""" ncols             4
3760 nrows             4
3761 xllcorner    148.00000
3762 yllcorner    -38.00000
3763 cellsize       0.25
3764 nodata_value   -9999.0
3765   -90.000    -10.000   -30.0 -20.000
3766   -10.000    99.000  -11.000 -30.000
3767   -40.000    66.000   22.000 -50.000
3768   -90.000    -10.000   -30.0 -20.000
3769""")
3770        fid.close()
3771       
3772        vcur_dir =  tempfile.mkdtemp()
3773        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3774        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3775       
3776        fid = open(vcur_dir_filename1, 'w')
3777        fid.write(""" ncols             4
3778 nrows             4
3779 xllcorner    148.00000
3780 yllcorner    -38.00000
3781 cellsize       0.25
3782 nodata_value   -9999.0
3783   -90.000    -10.000   -30.0 -20.000
3784   -10.000    80.000  -20.000 -30.000
3785   -40.000    50.000   10.000 -50.000
3786   -90.000    -10.000   -30.0 -20.000
3787""")
3788        fid.close()
3789        fid = open(vcur_dir_filename2, 'w')
3790        fid.write(""" ncols             4
3791 nrows             4
3792 xllcorner    148.00000
3793 yllcorner    -38.00000
3794 cellsize       0.25
3795 nodata_value   -9999.0
3796   -90.000    -10.000   -30.0 -20.000
3797   -10.000    88.000  -22.000 -30.000
3798   -40.000    55.000   11.000 -50.000
3799   -90.000    -10.000   -30.0 -20.000
3800""")
3801        fid.close()
3802       
3803        sww_file = tempfile.mktemp(".sww")
3804        #sww_file = 'a_test.sww'
3805        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3806                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3807                      mean_stage = 100,
3808                       minlat = -37.6, maxlat = -37.6,
3809                  minlon = 148.3, maxlon = 148.3
3810                      #,verbose = True
3811                      )
3812
3813        # check the sww file
3814       
3815        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3816        x = fid.variables['x'][:]
3817        y = fid.variables['y'][:]
3818        z = fid.variables['z'][:]
3819        stage = fid.variables['stage'][:]
3820        xmomentum = fid.variables['xmomentum'][:]
3821        ymomentum = fid.variables['ymomentum'][:]
3822        geo_ref = Geo_reference(NetCDFObject=fid)
3823        #print "geo_ref",geo_ref
3824        x_ref = geo_ref.get_xllcorner()
3825        y_ref = geo_ref.get_yllcorner()
3826        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3827
3828        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3829        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3830        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
3831
3832        #Easting:  632145.632  Northing: 5820863.269
3833        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3834       
3835        #print "x",x
3836        #print "y",y
3837        self.failUnless(len(x) == 4,'failed') # 2*2
3838        self.failUnless(len(x) == 4,'failed') # 2*2
3839       
3840        #Zone:   55
3841        #Easting:  632145.632  Northing: 5820863.269
3842        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3843        # magic number - y is close enough for me.
3844        assert allclose(x[3], 632145.63 - x_ref)
3845        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
3846
3847        assert allclose(z[0],9000.0 ) #z is elevation info
3848        #print "z",z
3849        # 2 time steps, 4 points
3850        self.failUnless(xmomentum.shape == (2,4), 'failed') 
3851        self.failUnless(ymomentum.shape == (2,4), 'failed') 
3852       
3853        #(100.0 - -1000.000)*10
3854        #assert allclose(xmomentum[0][5], 11000.0 )
3855       
3856        fid.close()
3857       
3858        # is the sww file readable?
3859        #Lets see if we can convert it to a dem!
3860        #print "sww_file",sww_file
3861        #dem_file = tempfile.mktemp(".dem")
3862        domain = sww2domain(sww_file) ###, dem_file)
3863        domain.check_integrity()
3864
3865        #tidy up
3866        os.remove(bath_dir_filename)
3867        os.rmdir(bath_dir)
3868
3869        os.remove(elevation_dir_filename1)
3870        os.remove(elevation_dir_filename2)
3871        os.rmdir(elevation_dir)
3872       
3873        os.remove(ucur_dir_filename1)
3874        os.remove(ucur_dir_filename2)
3875        os.rmdir(ucur_dir)
3876       
3877        os.remove(vcur_dir_filename1)
3878        os.remove(vcur_dir_filename2)
3879        os.rmdir(vcur_dir)
3880
3881
3882     
3883       
3884        # remove sww file
3885        os.remove(sww_file)
3886       
3887        # remove dem file
3888        #os.remove(dem_file)
3889
3890    def test_get_min_max_indexes(self):
3891        latitudes = [3,2,1,0]
3892        longitudes = [0,10,20,30]
3893
3894        # k - lat
3895        # l - lon
3896        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3897            latitudes,longitudes,
3898            -10,4,-10,31)
3899
3900        #print "kmin",kmin;print "kmax",kmax
3901        #print "lmin",lmin;print "lmax",lmax
3902        latitudes_new = latitudes[kmin:kmax]
3903        longitudes_news = longitudes[lmin:lmax]
3904        #print "latitudes_new", latitudes_new
3905        #print "longitudes_news",longitudes_news
3906        self.failUnless(latitudes == latitudes_new and \
3907                        longitudes == longitudes_news,
3908                         'failed')
3909       
3910        ## 2nd test
3911        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3912            latitudes,longitudes,
3913            0.5,2.5,5,25)
3914        #print "kmin",kmin;print "kmax",kmax
3915        #print "lmin",lmin;print "lmax",lmax
3916        latitudes_new = latitudes[kmin:kmax]
3917        longitudes_news = longitudes[lmin:lmax]
3918        #print "latitudes_new", latitudes_new
3919        #print "longitudes_news",longitudes_news
3920       
3921        self.failUnless(latitudes == latitudes_new and \
3922                        longitudes == longitudes_news,
3923                         'failed')
3924       
3925        ## 3rd test
3926        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
3927                                                                    longitudes,
3928                                                      1.1,1.9,12,17)
3929        #print "kmin",kmin;print "kmax",kmax
3930        #print "lmin",lmin;print "lmax",lmax
3931        latitudes_new = latitudes[kmin:kmax]
3932        longitudes_news = longitudes[lmin:lmax]
3933        #print "latitudes_new", latitudes_new
3934        #print "longitudes_news",longitudes_news
3935       
3936        self.failUnless(latitudes_new == [2, 1] and \
3937                        longitudes_news == [10, 20],
3938                         'failed')
3939
3940       
3941        ## 4th test
3942        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3943            latitudes,longitudes,
3944                                                      -0.1,1.9,-2,17)
3945        #print "kmin",kmin;print "kmax",kmax
3946        #print "lmin",lmin;print "lmax",lmax
3947        latitudes_new = latitudes[kmin:kmax]
3948        longitudes_news = longitudes[lmin:lmax]
3949        #print "latitudes_new", latitudes_new
3950        #print "longitudes_news",longitudes_news
3951       
3952        self.failUnless(latitudes_new == [2, 1, 0] and \
3953                        longitudes_news == [0, 10, 20],
3954                         'failed')
3955        ## 5th test
3956        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3957            latitudes,longitudes,
3958            0.1,1.9,2,17)
3959        #print "kmin",kmin;print "kmax",kmax
3960        #print "lmin",lmin;print "lmax",lmax
3961        latitudes_new = latitudes[kmin:kmax]
3962        longitudes_news = longitudes[lmin:lmax]
3963        #print "latitudes_new", latitudes_new
3964        #print "longitudes_news",longitudes_news
3965       
3966        self.failUnless(latitudes_new == [2, 1, 0] and \
3967                        longitudes_news == [0, 10, 20],
3968                         'failed')
3969       
3970        ## 6th test
3971       
3972        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3973            latitudes,longitudes,
3974            1.5,4,18,32)
3975        #print "kmin",kmin;print "kmax",kmax
3976        #print "lmin",lmin;print "lmax",lmax
3977        latitudes_new = latitudes[kmin:kmax]
3978        longitudes_news = longitudes[lmin:lmax]
3979        #print "latitudes_new", latitudes_new
3980        #print "longitudes_news",longitudes_news
3981       
3982        self.failUnless(latitudes_new == [3, 2, 1] and \
3983                        longitudes_news == [10, 20, 30],
3984                         'failed')
3985       
3986       
3987        ## 7th test
3988        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3989        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3990            latitudes,longitudes,
3991            1.5,1.5,15,15)
3992        #print "kmin",kmin;print "kmax",kmax
3993        #print "lmin",lmin;print "lmax",lmax
3994        latitudes_new = latitudes[kmin:kmax]
3995        longitudes_news = longitudes[lmin:lmax]
3996        m2d = m2d[kmin:kmax,lmin:lmax]
3997        #print "m2d", m2d
3998        #print "latitudes_new", latitudes_new
3999        #print "longitudes_news",longitudes_news
4000       
4001        self.failUnless(latitudes_new == [2, 1] and \
4002                        longitudes_news == [10, 20],
4003                         'failed')
4004       
4005        self.failUnless(m2d == [[5,6],[9,10]],
4006                         'failed')
4007       
4008    def test_get_min_max_indexes2(self):
4009        latitudes = [-30,-35,-40,-45]
4010        longitudes = [148,149,150,151]
4011
4012        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4013       
4014        # k - lat
4015        # l - lon
4016        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4017            latitudes,longitudes,
4018            -37,-27,147,149.5)
4019
4020        #print "kmin",kmin;print "kmax",kmax
4021        #print "lmin",lmin;print "lmax",lmax
4022        #print "m2d", m2d
4023        #print "latitudes", latitudes
4024        #print "longitudes",longitudes
4025        #print "latitudes[kmax]", latitudes[kmax]
4026        latitudes_new = latitudes[kmin:kmax]
4027        longitudes_new = longitudes[lmin:lmax]
4028        m2d = m2d[kmin:kmax,lmin:lmax]
4029        #print "m2d", m2d
4030        #print "latitudes_new", latitudes_new
4031        #print "longitudes_new",longitudes_new
4032       
4033        self.failUnless(latitudes_new == [-30, -35, -40] and \
4034                        longitudes_new == [148, 149,150],
4035                         'failed')
4036        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
4037                         'failed')
4038       
4039    def test_get_min_max_indexes3(self):
4040        latitudes = [-30,-35,-40,-45,-50,-55,-60]
4041        longitudes = [148,149,150,151]
4042
4043        # k - lat
4044        # l - lon
4045        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4046            latitudes,longitudes,
4047            -43,-37,148.5,149.5)
4048
4049       
4050        #print "kmin",kmin;print "kmax",kmax
4051        #print "lmin",lmin;print "lmax",lmax
4052        #print "latitudes", latitudes
4053        #print "longitudes",longitudes
4054        latitudes_new = latitudes[kmin:kmax]
4055        longitudes_news = longitudes[lmin:lmax]
4056        #print "latitudes_new", latitudes_new
4057        #print "longitudes_news",longitudes_news
4058       
4059        self.failUnless(latitudes_new == [-35, -40, -45] and \
4060                        longitudes_news == [148, 149,150],
4061                         'failed')
4062       
4063    def test_get_min_max_indexes4(self):
4064        latitudes = [-30,-35,-40,-45,-50,-55,-60]
4065        longitudes = [148,149,150,151]
4066
4067        # k - lat
4068        # l - lon
4069        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4070            latitudes,longitudes)
4071       
4072        #print "kmin",kmin;print "kmax",kmax
4073        #print "lmin",lmin;print "lmax",lmax
4074        #print "latitudes", latitudes
4075        #print "longitudes",longitudes
4076        latitudes_new = latitudes[kmin:kmax]
4077        longitudes_news = longitudes[lmin:lmax]
4078        #print "latitudes_new", latitudes_new
4079        #print "longitudes_news",longitudes_news
4080       
4081        self.failUnless(latitudes_new == latitudes  and \
4082                        longitudes_news == longitudes,
4083                         'failed')
4084       
4085    def test_tsh2sww(self):
4086        import os
4087        import tempfile
4088       
4089        tsh_file = tempfile.mktemp(".tsh")
4090        file = open(tsh_file,"w")
4091        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
40920 0.0 0.0 0.0 0.0 0.01 \n \
40931 1.0 0.0 10.0 10.0 0.02  \n \
40942 0.0 1.0 0.0 10.0 0.03  \n \
40953 0.5 0.25 8.0 12.0 0.04  \n \
4096# Vert att title  \n \
4097elevation  \n \
4098stage  \n \
4099friction  \n \
41002 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
41010 0 3 2 -1  -1  1 dsg\n\
41021 0 1 3 -1  0 -1   ole nielsen\n\
41034 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
41040 1 0 2 \n\
41051 0 2 3 \n\
41062 2 3 \n\
41073 3 1 1 \n\
41083 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
41090 216.0 -86.0 \n \
41101 160.0 -167.0 \n \
41112 114.0 -91.0 \n \
41123 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
41130 0 1 0 \n \
41141 1 2 0 \n \
41152 2 0 0 \n \
41160 # <x> <y> ...Mesh Holes... \n \
41170 # <x> <y> <attribute>...Mesh Regions... \n \
41180 # <x> <y> <attribute>...Mesh Regions, area... \n\
4119#Geo reference \n \
412056 \n \
4121140 \n \
4122120 \n")
4123        file.close()
4124
4125        #sww_file = tempfile.mktemp(".sww")
4126        #print "sww_file",sww_file
4127        #print "sww_file",tsh_file
4128        tsh2sww(tsh_file)
4129
4130        os.remove(tsh_file)
4131        os.remove(tsh_file[:-4] + '.sww')
4132#-------------------------------------------------------------
4133if __name__ == "__main__":
4134    #suite = unittest.makeSuite(Test_Data_Manager,'test_tsh2sww')
4135    suite = unittest.makeSuite(Test_Data_Manager,'test')
4136    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
4137    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')   
4138    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
4139    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
4140    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
4141    runner = unittest.TextTestRunner()
4142    runner.run(suite)
4143 
4144
4145
4146
4147    def test_dem2pts(self):
4148        """Test conversion from dem in ascii format to native NetCDF xya format
4149        """
4150
4151        import time, os
4152        from Numeric import array, zeros, allclose, Float, concatenate
4153        from Scientific.IO.NetCDF import NetCDFFile
4154
4155        #Write test asc file
4156        root = 'demtest'
4157
4158        filename = root+'.asc'
4159        fid = open(filename, 'w')
4160        fid.write("""ncols         5
4161nrows         6
4162xllcorner     2000.5
4163yllcorner     3000.5
4164cellsize      25
4165NODATA_value  -9999
4166""")
4167        #Create linear function
4168
4169        ref_points = []
4170        ref_elevation = []
4171        for i in range(6):
4172            y = (6-i)*25.0
4173            for j in range(5):
4174                x = j*25.0
4175                z = x+2*y
4176
4177                ref_points.append( [x,y] )
4178                ref_elevation.append(z)
4179                fid.write('%f ' %z)
4180            fid.write('\n')
4181
4182        fid.close()
4183
4184        #Write prj file with metadata
4185        metafilename = root+'.prj'
4186        fid = open(metafilename, 'w')
4187
4188
4189        fid.write("""Projection UTM
4190Zone 56
4191Datum WGS84
4192Zunits NO
4193Units METERS
4194Spheroid WGS84
4195Xshift 0.0000000000
4196Yshift 10000000.0000000000
4197Parameters
4198""")
4199        fid.close()
4200
4201        #Convert to NetCDF pts
4202        convert_dem_from_ascii2netcdf(root)
4203        dem2pts(root)
4204
4205        #Check contents
4206        #Get NetCDF
4207        fid = NetCDFFile(root+'.pts', 'r')
4208
4209        # Get the variables
4210        #print fid.variables.keys()
4211        points = fid.variables['points']
4212        elevation = fid.variables['elevation']
4213
4214        #Check values
4215
4216        #print points[:]
4217        #print ref_points
4218        assert allclose(points, ref_points)
4219
4220        #print attributes[:]
4221        #print ref_elevation
4222        assert allclose(elevation, ref_elevation)
4223
4224        #Cleanup
4225        fid.close()
4226
4227
4228        os.remove(root + '.pts')
4229        os.remove(root + '.dem')
4230        os.remove(root + '.asc')
4231        os.remove(root + '.prj')
4232
4233
4234
4235    def test_dem2pts_bounding_box(self):
4236        """Test conversion from dem in ascii format to native NetCDF xya format
4237        """
4238
4239        import time, os
4240        from Numeric import array, zeros, allclose, Float, concatenate
4241        from Scientific.IO.NetCDF import NetCDFFile
4242
4243        #Write test asc file
4244        root = 'demtest'
4245
4246        filename = root+'.asc'
4247        fid = open(filename, 'w')
4248        fid.write("""ncols         5
4249nrows         6
4250xllcorner     2000.5
4251yllcorner     3000.5
4252cellsize      25
4253NODATA_value  -9999
4254""")
4255        #Create linear function
4256
4257        ref_points = []
4258        ref_elevation = []
4259        for i in range(6):
4260            y = (6-i)*25.0
4261            for j in range(5):
4262                x = j*25.0
4263                z = x+2*y
4264
4265                ref_points.append( [x,y] )
4266                ref_elevation.append(z)
4267                fid.write('%f ' %z)
4268            fid.write('\n')
4269
4270        fid.close()
4271
4272        #Write prj file with metadata
4273        metafilename = root+'.prj'
4274        fid = open(metafilename, 'w')
4275
4276
4277        fid.write("""Projection UTM
4278Zone 56
4279Datum WGS84
4280Zunits NO
4281Units METERS
4282Spheroid WGS84
4283Xshift 0.0000000000
4284Yshift 10000000.0000000000
4285Parameters
4286""")
4287        fid.close()
4288
4289        #Convert to NetCDF pts
4290        convert_dem_from_ascii2netcdf(root)
4291        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4292                northing_min=3035.0, northing_max=3125.5)
4293
4294        #Check contents
4295        #Get NetCDF
4296        fid = NetCDFFile(root+'.pts', 'r')
4297
4298        # Get the variables
4299        #print fid.variables.keys()
4300        points = fid.variables['points']
4301        elevation = fid.variables['elevation']
4302
4303        #Check values
4304        assert fid.xllcorner[0] == 2010.0
4305        assert fid.yllcorner[0] == 3035.0
4306
4307        #create new reference points
4308        ref_points = []
4309        ref_elevation = []
4310        for i in range(4):
4311            y = (4-i)*25.0 + 25.0
4312            y_new = y + 3000.5 - 3035.0
4313            for j in range(4):
4314                x = j*25.0 + 25.0
4315                x_new = x + 2000.5 - 2010.0
4316                z = x+2*y
4317
4318                ref_points.append( [x_new,y_new] )
4319                ref_elevation.append(z)
4320
4321        #print points[:]
4322        #print ref_points
4323        assert allclose(points, ref_points)
4324
4325        #print attributes[:]
4326        #print ref_elevation
4327        assert allclose(elevation, ref_elevation)
4328
4329        #Cleanup
4330        fid.close()
4331
4332
4333        os.remove(root + '.pts')
4334        os.remove(root + '.dem')
4335        os.remove(root + '.asc')
4336        os.remove(root + '.prj')
4337
4338
4339
4340    def test_dem2pts_remove_Nullvalues(self):
4341        """Test conversion from dem in ascii format to native NetCDF xya format
4342        """
4343
4344        import time, os
4345        from Numeric import array, zeros, allclose, Float, concatenate
4346        from Scientific.IO.NetCDF import NetCDFFile
4347
4348        #Write test asc file
4349        root = 'demtest'
4350
4351        filename = root+'.asc'
4352        fid = open(filename, 'w')
4353        fid.write("""ncols         5
4354nrows         6
4355xllcorner     2000.5
4356yllcorner     3000.5
4357cellsize      25
4358NODATA_value  -9999
4359""")
4360        #Create linear function
4361        # ref_ will write all the values
4362        # new_ref_ will write the values except for NODATA_values
4363        ref_points = []
4364        ref_elevation = []
4365        new_ref_pts = []
4366        new_ref_elev = []
4367        NODATA_value = -9999
4368        for i in range(6):
4369            y = (6-i)*25.0
4370            for j in range(5):
4371                x = j*25.0
4372                z = x+2*y
4373                if j == 4: z = NODATA_value # column
4374                if i == 2 and j == 2: z = NODATA_value # random
4375                if i == 5 and j == 1: z = NODATA_value
4376                if i == 1: z = NODATA_value # row
4377                if i == 3 and j == 1: z = NODATA_value # two pts/row
4378                if i == 3 and j == 3: z = NODATA_value
4379
4380               
4381                if z <> NODATA_value:
4382                    new_ref_elev.append(z)
4383                    new_ref_pts.append( [x,y] )
4384               
4385                ref_points.append( [x,y] )
4386                ref_elevation.append(z)
4387
4388                fid.write('%f ' %z)
4389            fid.write('\n')
4390
4391        fid.close()
4392
4393       
4394        #Write prj file with metadata
4395        metafilename = root+'.prj'
4396        fid = open(metafilename, 'w')
4397
4398
4399        fid.write("""Projection UTM
4400Zone 56
4401Datum WGS84
4402Zunits NO
4403Units METERS
4404Spheroid WGS84
4405Xshift 0.0000000000
4406Yshift 10000000.0000000000
4407Parameters
4408""")
4409        fid.close()
4410
4411        #Convert to NetCDF pts
4412        convert_dem_from_ascii2netcdf(root)
4413        dem2pts(root)
4414
4415        #Check contents
4416        #Get NetCDF
4417        fid = NetCDFFile(root+'.pts', 'r')
4418
4419        # Get the variables
4420        #print fid.variables.keys()
4421        points = fid.variables['points']
4422        elevation = fid.variables['elevation']
4423
4424        #Check values
4425        #print 'points', points[:]
4426        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
4427        assert allclose(points, new_ref_pts), 'points do not align'
4428
4429        #print 'elevation', elevation[:]
4430        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
4431        assert allclose(elevation, new_ref_elev), 'elevations do not align'
4432
4433        #Cleanup
4434        fid.close()
4435
4436
4437        os.remove(root + '.pts')
4438        os.remove(root + '.dem')
4439        os.remove(root + '.asc')
4440        os.remove(root + '.prj')
4441
4442    def test_dem2pts_bounding_box_Nullvalues(self):
4443        """Test conversion from dem in ascii format to native NetCDF xya format
4444        """
4445
4446        import time, os
4447        from Numeric import array, zeros, allclose, Float, concatenate
4448        from Scientific.IO.NetCDF import NetCDFFile
4449
4450        #Write test asc file
4451        root = 'demtest'
4452
4453        filename = root+'.asc'
4454        fid = open(filename, 'w')
4455        fid.write("""ncols         5
4456nrows         6
4457xllcorner     2000.5
4458yllcorner     3000.5
4459cellsize      25
4460NODATA_value  -9999
4461""")
4462        #Create linear function
4463
4464        ref_points = []
4465        ref_elevation = []
4466        new_ref_pts1 = []
4467        new_ref_elev1 = []
4468        NODATA_value = -9999
4469        for i in range(6):
4470            y = (6-i)*25.0
4471            for j in range(5):
4472                x = j*25.0
4473                z = x+2*y
4474                if j == 4: z = NODATA_value # column
4475                if i == 2 and j == 2: z = NODATA_value # random
4476                if i == 5 and j == 1: z = NODATA_value
4477                if i == 1: z = NODATA_value # row
4478                if i == 3 and j == 1: z = NODATA_value # two pts/row
4479                if i == 3 and j == 3: z = NODATA_value
4480
4481                if z <> NODATA_value:
4482                    new_ref_elev1.append(z)
4483                    new_ref_pts1.append( [x,y] )
4484                   
4485                ref_points.append( [x,y] )
4486                ref_elevation.append(z)
4487                fid.write('%f ' %z)
4488            fid.write('\n')
4489
4490        fid.close()
4491
4492        #Write prj file with metadata
4493        metafilename = root+'.prj'
4494        fid = open(metafilename, 'w')
4495
4496
4497        fid.write("""Projection UTM
4498Zone 56
4499Datum WGS84
4500Zunits NO
4501Units METERS
4502Spheroid WGS84
4503Xshift 0.0000000000
4504Yshift 10000000.0000000000
4505Parameters
4506""")
4507        fid.close()
4508       
4509        #Convert to NetCDF pts
4510        convert_dem_from_ascii2netcdf(root)
4511        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4512                northing_min=3035.0, northing_max=3125.5)
4513
4514        #Check contents
4515        #Get NetCDF
4516        fid = NetCDFFile(root+'.pts', 'r')
4517
4518        # Get the variables
4519        #print fid.variables.keys()
4520        points = fid.variables['points']
4521        elevation = fid.variables['elevation']
4522
4523        #Check values
4524        assert fid.xllcorner[0] == 2010.0
4525        assert fid.yllcorner[0] == 3035.0
4526
4527        #create new reference points
4528        ref_points = []
4529        ref_elevation = []
4530        new_ref_pts2 = []
4531        new_ref_elev2 = []
4532        for i in range(4):
4533            y = (4-i)*25.0 + 25.0
4534            y_new = y + 3000.5 - 3035.0
4535            for j in range(4):
4536                x = j*25.0 + 25.0
4537                x_new = x + 2000.5 - 2010.0
4538                z = x+2*y
4539
4540                if j == 3: z = NODATA_value # column
4541                if i == 1 and j == 1: z = NODATA_value # random
4542                if i == 4 and j == 0: z = NODATA_value
4543                if i == 0: z = NODATA_value # row
4544                if i == 2 and j == 0: z = NODATA_value # two pts/row
4545                if i == 2 and j == 2: z = NODATA_value
4546
4547                if z <> NODATA_value:
4548                    new_ref_elev2.append(z)
4549                    new_ref_pts2.append( [x_new,y_new] )
4550
4551                   
4552                ref_points.append( [x_new,y_new] )
4553                ref_elevation.append(z)
4554
4555        #print points[:]
4556        #print ref_points
4557        #assert allclose(points, ref_points)
4558
4559        #print attributes[:]
4560        #print ref_elevation
4561        #assert allclose(elevation, ref_elevation)
4562
4563       
4564        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
4565        assert allclose(points, new_ref_pts2), 'points do not align'
4566
4567        #print 'elevation', elevation[:]
4568        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
4569        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
4570        #Cleanup
4571        fid.close()
4572
4573
4574        os.remove(root + '.pts')
4575        os.remove(root + '.dem')
4576        os.remove(root + '.asc')
4577        os.remove(root + '.prj')
4578
4579
4580
4581
4582
Note: See TracBrowser for help on using the repository browser.