source: inundation/pyvolution/test_data_manager.py @ 2650

Last change on this file since 2650 was 2650, checked in by steve, 18 years ago
File size: 140.6 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        long_name = 'LON'
2378        lat_name = 'LAT'
2379        time_name = 'TIME'
2380
2381        nx = 3
2382        ny = 2
2383
2384        for fid in [fid1,fid2,fid3]:
2385            fid.createDimension(long_name,nx)
2386            fid.createVariable(long_name,'d',(long_name,))
2387            fid.variables[long_name].point_spacing='uneven'
2388            fid.variables[long_name].units='degrees_east'
2389            fid.variables[long_name].assignValue(h1_list)
2390
2391            fid.createDimension(lat_name,ny)
2392            fid.createVariable(lat_name,'d',(lat_name,))
2393            fid.variables[lat_name].point_spacing='uneven'
2394            fid.variables[lat_name].units='degrees_north'
2395            fid.variables[lat_name].assignValue(h2_list)
2396
2397            fid.createDimension(time_name,2)
2398            fid.createVariable(time_name,'d',(time_name,))
2399            fid.variables[time_name].point_spacing='uneven'
2400            fid.variables[time_name].units='seconds'
2401            fid.variables[time_name].assignValue([0.,1.])
2402            if fid == fid3: break
2403
2404
2405        for fid in [fid4]:
2406            fid.createDimension(long_name,nx)
2407            fid.createVariable(long_name,'d',(long_name,))
2408            fid.variables[long_name].point_spacing='uneven'
2409            fid.variables[long_name].units='degrees_east'
2410            fid.variables[long_name].assignValue(h1_list)
2411
2412            fid.createDimension(lat_name,ny)
2413            fid.createVariable(lat_name,'d',(lat_name,))
2414            fid.variables[lat_name].point_spacing='uneven'
2415            fid.variables[lat_name].units='degrees_north'
2416            fid.variables[lat_name].assignValue(h2_list)
2417
2418        name = {}
2419        name[fid1]='HA'
2420        name[fid2]='UA'
2421        name[fid3]='VA'
2422        name[fid4]='ELEVATION'
2423
2424        units = {}
2425        units[fid1]='cm'
2426        units[fid2]='cm/s'
2427        units[fid3]='cm/s'
2428        units[fid4]='m'
2429
2430        values = {}
2431        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2432        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2433        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2434        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2435
2436        for fid in [fid1,fid2,fid3]:
2437          fid.createVariable(name[fid],'d',(time_name,lat_name,long_name))
2438          fid.variables[name[fid]].point_spacing='uneven'
2439          fid.variables[name[fid]].units=units[fid]
2440          fid.variables[name[fid]].assignValue(values[fid])
2441          fid.variables[name[fid]].missing_value = -99999999.
2442          if fid == fid3: break
2443
2444        for fid in [fid4]:
2445            fid.createVariable(name[fid],'d',(lat_name,long_name))
2446            fid.variables[name[fid]].point_spacing='uneven'
2447            fid.variables[name[fid]].units=units[fid]
2448            fid.variables[name[fid]].assignValue(values[fid])
2449            fid.variables[name[fid]].missing_value = -99999999.
2450
2451
2452        fid1.sync(); fid1.close()
2453        fid2.sync(); fid2.close()
2454        fid3.sync(); fid3.close()
2455        fid4.sync(); fid4.close()
2456
2457        fid1 = NetCDFFile('test_ha.nc','r')
2458        fid2 = NetCDFFile('test_e.nc','r')
2459        fid3 = NetCDFFile('test_va.nc','r')
2460
2461
2462        first_amp = fid1.variables['HA'][:][0,0,0]
2463        third_amp = fid1.variables['HA'][:][0,0,2]
2464        first_elevation = fid2.variables['ELEVATION'][0,0]
2465        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2466        first_speed = fid3.variables['VA'][0,0,0]
2467        third_speed = fid3.variables['VA'][:][0,0,2]
2468
2469        fid1.close()
2470        fid2.close()
2471        fid3.close()
2472
2473        #Call conversion (with zero origin)
2474        ferret2sww('test', verbose=False, origin = (56, 0, 0))
2475
2476        os.remove('test_va.nc')
2477        os.remove('test_ua.nc')
2478        os.remove('test_ha.nc')
2479        os.remove('test_e.nc')
2480
2481        #Read output file 'test.sww'
2482        fid = NetCDFFile('test.sww')
2483
2484
2485        #Check first value
2486        elevation = fid.variables['elevation'][:]
2487        stage = fid.variables['stage'][:]
2488        xmomentum = fid.variables['xmomentum'][:]
2489        ymomentum = fid.variables['ymomentum'][:]
2490
2491        #print ymomentum
2492        first_height = first_amp/100 - first_elevation
2493        third_height = third_amp/100 - third_elevation
2494        first_momentum=first_speed*first_height/100
2495        third_momentum=third_speed*third_height/100
2496
2497        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2498        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2499
2500        fid.close()
2501
2502        #Cleanup
2503        os.remove('test.sww')
2504
2505
2506
2507
2508    def test_ferret2sww_nz_origin(self):
2509        from Scientific.IO.NetCDF import NetCDFFile
2510        from coordinate_transforms.redfearn import redfearn
2511
2512        #Call conversion (with nonzero origin)
2513        ferret2sww(self.test_MOST_file, verbose=False,
2514                   origin = (56, 100000, 200000))
2515
2516
2517        #Work out the UTM coordinates for first point
2518        zone, e, n = redfearn(-34.5, 150.66667)
2519
2520        #Read output file 'small.sww'
2521        #fid = NetCDFFile('small.sww', 'r')
2522        fid = NetCDFFile(self.test_MOST_file + '.sww')
2523
2524        x = fid.variables['x'][:]
2525        y = fid.variables['y'][:]
2526
2527        #Check that first coordinate is correctly represented
2528        assert allclose(x[0], e-100000)
2529        assert allclose(y[0], n-200000)
2530
2531        fid.close()
2532
2533        #Cleanup
2534        os.remove(self.test_MOST_file + '.sww')
2535
2536
2537
2538    def test_sww_extent(self):
2539        """Not a test, rather a look at the sww format
2540        """
2541
2542        import time, os
2543        from Numeric import array, zeros, allclose, Float, concatenate
2544        from Scientific.IO.NetCDF import NetCDFFile
2545
2546        self.domain.filename = 'datatest' + str(id(self))
2547        self.domain.format = 'sww'
2548        self.domain.smooth = True
2549        self.domain.reduction = mean
2550        self.domain.set_datadir('.')
2551
2552
2553        sww = get_dataobject(self.domain)
2554        sww.store_connectivity()
2555        sww.store_timestep('stage')
2556        self.domain.time = 2.
2557
2558        #Modify stage at second timestep
2559        stage = self.domain.quantities['stage'].vertex_values
2560        self.domain.set_quantity('stage', stage/2)
2561
2562        sww.store_timestep('stage')
2563
2564        file_and_extension_name = self.domain.filename + ".sww"
2565        #print "file_and_extension_name",file_and_extension_name
2566        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2567               extent_sww(file_and_extension_name )
2568
2569        assert allclose(xmin, 0.0)
2570        assert allclose(xmax, 1.0)
2571        assert allclose(ymin, 0.0)
2572        assert allclose(ymax, 1.0)
2573        assert allclose(stagemin, -0.85)
2574        assert allclose(stagemax, 0.15)
2575
2576
2577        #Cleanup
2578        os.remove(sww.filename)
2579
2580
2581
2582    def test_sww2domain1(self):
2583        ################################################
2584        #Create a test domain, and evolve and save it.
2585        ################################################
2586        from mesh_factory import rectangular
2587        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2588            Constant_height, Time_boundary, Transmissive_boundary
2589        from Numeric import array
2590
2591        #Create basic mesh
2592
2593        yiel=0.01
2594        points, vertices, boundary = rectangular(10,10)
2595
2596        #Create shallow water domain
2597        domain = Domain(points, vertices, boundary)
2598        domain.geo_reference = Geo_reference(56,11,11)
2599        domain.smooth = False
2600        domain.visualise = False
2601        domain.store = True
2602        domain.filename = 'bedslope'
2603        domain.default_order=2
2604        #Bed-slope and friction
2605        domain.set_quantity('elevation', lambda x,y: -x/3)
2606        domain.set_quantity('friction', 0.1)
2607        # Boundary conditions
2608        from math import sin, pi
2609        Br = Reflective_boundary(domain)
2610        Bt = Transmissive_boundary(domain)
2611        Bd = Dirichlet_boundary([0.2,0.,0.])
2612        Bw = Time_boundary(domain=domain,f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2613
2614        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2615        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2616
2617        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2618        #Initial condition
2619        h = 0.05
2620        elevation = domain.quantities['elevation'].vertex_values
2621        domain.set_quantity('stage', elevation + h)
2622
2623        domain.check_integrity()
2624        #Evolution
2625        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2626            #domain.write_time()
2627            pass
2628
2629
2630        ##########################################
2631        #Import the example's file as a new domain
2632        ##########################################
2633        from data_manager import sww2domain
2634        from Numeric import allclose
2635        import os
2636
2637        filename = domain.datadir+os.sep+domain.filename+'.sww'
2638        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2639        #points, vertices, boundary = rectangular(15,15)
2640        #domain2.boundary = boundary
2641        ###################
2642        ##NOW TEST IT!!!
2643        ###################
2644
2645        #os.remove(domain.filename + '.sww')
2646        os.remove(filename)
2647
2648        bits = ['vertex_coordinates']
2649        for quantity in ['elevation']+domain.quantities_to_be_stored:
2650            bits.append('get_quantity("%s").get_integral()' %quantity)
2651            bits.append('get_quantity("%s").get_values()' %quantity)
2652
2653        for bit in bits:
2654            #print 'testing that domain.'+bit+' has been restored'
2655            #print bit
2656            #print 'done'
2657            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2658
2659        ######################################
2660        #Now evolve them both, just to be sure
2661        ######################################x
2662        visualise = False
2663        #visualise = True
2664        domain.visualise = visualise
2665        domain.time = 0.
2666        from time import sleep
2667
2668        final = .1
2669        domain.set_quantity('friction', 0.1)
2670        domain.store = False
2671        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2672
2673
2674        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2675            if visualise: sleep(1.)
2676            #domain.write_time()
2677            pass
2678
2679        final = final - (domain2.starttime-domain.starttime)
2680        #BUT since domain1 gets time hacked back to 0:
2681        final = final + (domain2.starttime-domain.starttime)
2682
2683        domain2.smooth = False
2684        domain2.visualise = visualise
2685        domain2.store = False
2686        domain2.default_order=2
2687        domain2.set_quantity('friction', 0.1)
2688        #Bed-slope and friction
2689        # Boundary conditions
2690        Bd2=Dirichlet_boundary([0.2,0.,0.])
2691        domain2.boundary = domain.boundary
2692        #print 'domain2.boundary'
2693        #print domain2.boundary
2694        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2695        #domain2.set_boundary({'exterior': Bd})
2696
2697        domain2.check_integrity()
2698
2699        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2700            if visualise: sleep(1.)
2701            #domain2.write_time()
2702            pass
2703
2704        ###################
2705        ##NOW TEST IT!!!
2706        ##################
2707
2708        bits = [ 'vertex_coordinates']
2709
2710        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2711            bits.append('get_quantity("%s").get_integral()' %quantity)
2712            bits.append('get_quantity("%s").get_values()' %quantity)
2713
2714        #print bits
2715        for bit in bits:
2716            #print bit
2717            #print eval('domain.'+bit)
2718            #print eval('domain2.'+bit)
2719            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2720
2721
2722    def test_sww2domain2(self):
2723        ##################################################################
2724        #Same as previous test, but this checks how NaNs are handled.
2725        ##################################################################
2726
2727
2728        from mesh_factory import rectangular
2729        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2730             Constant_height, Time_boundary, Transmissive_boundary
2731        from Numeric import array
2732
2733        #Create basic mesh
2734        points, vertices, boundary = rectangular(2,2)
2735
2736        #Create shallow water domain
2737        domain = Domain(points, vertices, boundary)
2738        domain.smooth = False
2739        domain.visualise = False
2740        domain.store = True
2741        domain.set_name('test_file')
2742        domain.set_datadir('.')
2743        domain.default_order=2
2744        domain.quantities_to_be_stored=['stage']
2745
2746        domain.set_quantity('elevation', lambda x,y: -x/3)
2747        domain.set_quantity('friction', 0.1)
2748
2749        from math import sin, pi
2750        Br = Reflective_boundary(domain)
2751        Bt = Transmissive_boundary(domain)
2752        Bd = Dirichlet_boundary([0.2,0.,0.])
2753        Bw = Time_boundary(domain=domain,
2754                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2755
2756        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2757
2758        h = 0.05
2759        elevation = domain.quantities['elevation'].vertex_values
2760        domain.set_quantity('stage', elevation + h)
2761
2762        domain.check_integrity()
2763
2764        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2765            pass
2766            #domain.write_time()
2767
2768
2769
2770        ##################################
2771        #Import the file as a new domain
2772        ##################################
2773        from data_manager import sww2domain
2774        from Numeric import allclose
2775        import os
2776
2777        filename = domain.datadir+os.sep+domain.filename+'.sww'
2778
2779        #Fail because NaNs are present
2780        try:
2781            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2782        except:
2783            #Now import it, filling NaNs to be 0
2784            filler = 0
2785            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2786
2787        #Clean up
2788        os.remove(filename)
2789
2790
2791        bits = [ 'geo_reference.get_xllcorner()',
2792                'geo_reference.get_yllcorner()',
2793                'vertex_coordinates']
2794
2795        for quantity in ['elevation']+domain.quantities_to_be_stored:
2796            bits.append('get_quantity("%s").get_integral()' %quantity)
2797            bits.append('get_quantity("%s").get_values()' %quantity)
2798
2799        for bit in bits:
2800        #    print 'testing that domain.'+bit+' has been restored'
2801            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2802
2803        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
2804        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
2805        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
2806        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
2807
2808
2809
2810    #def test_weed(self):
2811        from data_manager import weed
2812
2813        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2814        volumes1 = [[0,1,2],[3,4,5]]
2815        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2816        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2817
2818        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2819
2820        assert len(points2)==len(coordinates2)
2821        for i in range(len(coordinates2)):
2822            coordinate = tuple(coordinates2[i])
2823            assert points2.has_key(coordinate)
2824            points2[coordinate]=i
2825
2826        for triangle in volumes1:
2827            for coordinate in triangle:
2828                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2829                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2830
2831
2832    #FIXME This fails - smooth makes the comparism too hard for allclose
2833    def ztest_sww2domain3(self):
2834        ################################################
2835        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2836        ################################################
2837        from mesh_factory import rectangular
2838        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2839             Constant_height, Time_boundary, Transmissive_boundary
2840        from Numeric import array
2841        #Create basic mesh
2842
2843        yiel=0.01
2844        points, vertices, boundary = rectangular(10,10)
2845
2846        #Create shallow water domain
2847        domain = Domain(points, vertices, boundary)
2848        domain.geo_reference = Geo_reference(56,11,11)
2849        domain.smooth = True
2850        domain.visualise = False
2851        domain.store = True
2852        domain.filename = 'bedslope'
2853        domain.default_order=2
2854        #Bed-slope and friction
2855        domain.set_quantity('elevation', lambda x,y: -x/3)
2856        domain.set_quantity('friction', 0.1)
2857        # Boundary conditions
2858        from math import sin, pi
2859        Br = Reflective_boundary(domain)
2860        Bt = Transmissive_boundary(domain)
2861        Bd = Dirichlet_boundary([0.2,0.,0.])
2862        Bw = Time_boundary(domain=domain,
2863                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2864
2865        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2866
2867        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2868        #Initial condition
2869        h = 0.05
2870        elevation = domain.quantities['elevation'].vertex_values
2871        domain.set_quantity('stage', elevation + h)
2872
2873
2874        domain.check_integrity()
2875        #Evolution
2876        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2877        #    domain.write_time()
2878            pass
2879
2880
2881        ##########################################
2882        #Import the example's file as a new domain
2883        ##########################################
2884        from data_manager import sww2domain
2885        from Numeric import allclose
2886        import os
2887
2888        filename = domain.datadir+os.sep+domain.filename+'.sww'
2889        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2890        #points, vertices, boundary = rectangular(15,15)
2891        #domain2.boundary = boundary
2892        ###################
2893        ##NOW TEST IT!!!
2894        ###################
2895
2896        os.remove(domain.filename + '.sww')
2897
2898        #FIXME smooth domain so that they can be compared
2899
2900
2901        bits = []#'vertex_coordinates']
2902        for quantity in ['elevation']+domain.quantities_to_be_stored:
2903            bits.append('quantities["%s"].get_integral()'%quantity)
2904
2905
2906        for bit in bits:
2907            #print 'testing that domain.'+bit+' has been restored'
2908            #print bit
2909            #print 'done'
2910            #print ('domain.'+bit), eval('domain.'+bit)
2911            #print ('domain2.'+bit), eval('domain2.'+bit)
2912            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2913            pass
2914
2915        ######################################
2916        #Now evolve them both, just to be sure
2917        ######################################x
2918        visualise = False
2919        visualise = True
2920        domain.visualise = visualise
2921        domain.time = 0.
2922        from time import sleep
2923
2924        final = .5
2925        domain.set_quantity('friction', 0.1)
2926        domain.store = False
2927        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2928
2929        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2930            if visualise: sleep(.03)
2931            #domain.write_time()
2932            pass
2933
2934        domain2.smooth = True
2935        domain2.visualise = visualise
2936        domain2.store = False
2937        domain2.default_order=2
2938        domain2.set_quantity('friction', 0.1)
2939        #Bed-slope and friction
2940        # Boundary conditions
2941        Bd2=Dirichlet_boundary([0.2,0.,0.])
2942        Br2 = Reflective_boundary(domain2)
2943        domain2.boundary = domain.boundary
2944        #print 'domain2.boundary'
2945        #print domain2.boundary
2946        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2947        #domain2.boundary = domain.boundary
2948        #domain2.set_boundary({'exterior': Bd})
2949
2950        domain2.check_integrity()
2951
2952        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2953            if visualise: sleep(.03)
2954            #domain2.write_time()
2955            pass
2956
2957        ###################
2958        ##NOW TEST IT!!!
2959        ##################
2960
2961        print '><><><><>>'
2962        bits = [ 'vertex_coordinates']
2963
2964        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2965            #bits.append('quantities["%s"].get_integral()'%quantity)
2966            bits.append('get_quantity("%s").get_values()' %quantity)
2967
2968        for bit in bits:
2969            print bit
2970            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2971
2972
2973    def test_decimate_dem(self):
2974        """Test decimation of dem file
2975        """
2976
2977        import os
2978        from Numeric import ones, allclose, Float, arange
2979        from Scientific.IO.NetCDF import NetCDFFile
2980
2981        #Write test dem file
2982        root = 'decdemtest'
2983
2984        filename = root + '.dem'
2985        fid = NetCDFFile(filename, 'w')
2986
2987        fid.institution = 'Geoscience Australia'
2988        fid.description = 'NetCDF DEM format for compact and portable ' +\
2989                          'storage of spatial point data'
2990
2991        nrows = 15
2992        ncols = 18
2993
2994        fid.ncols = ncols
2995        fid.nrows = nrows
2996        fid.xllcorner = 2000.5
2997        fid.yllcorner = 3000.5
2998        fid.cellsize = 25
2999        fid.NODATA_value = -9999
3000
3001        fid.zone = 56
3002        fid.false_easting = 0.0
3003        fid.false_northing = 0.0
3004        fid.projection = 'UTM'
3005        fid.datum = 'WGS84'
3006        fid.units = 'METERS'
3007
3008        fid.createDimension('number_of_points', nrows*ncols)
3009
3010        fid.createVariable('elevation', Float, ('number_of_points',))
3011
3012        elevation = fid.variables['elevation']
3013
3014        elevation[:] = (arange(nrows*ncols))
3015
3016        fid.close()
3017
3018        #generate the elevation values expected in the decimated file
3019        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
3020                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3021                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3022                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3023                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3024                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
3025                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3026                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3027                         (144+145+146+162+163+164+180+181+182) / 9.0,
3028                         (148+149+150+166+167+168+184+185+186) / 9.0,
3029                         (152+153+154+170+171+172+188+189+190) / 9.0,
3030                         (156+157+158+174+175+176+192+193+194) / 9.0,
3031                         (216+217+218+234+235+236+252+253+254) / 9.0,
3032                         (220+221+222+238+239+240+256+257+258) / 9.0,
3033                         (224+225+226+242+243+244+260+261+262) / 9.0,
3034                         (228+229+230+246+247+248+264+265+266) / 9.0]
3035
3036        #generate a stencil for computing the decimated values
3037        stencil = ones((3,3), Float) / 9.0
3038
3039        decimate_dem(root, stencil=stencil, cellsize_new=100)
3040
3041        #Open decimated NetCDF file
3042        fid = NetCDFFile(root + '_100.dem', 'r')
3043
3044        # Get decimated elevation
3045        elevation = fid.variables['elevation']
3046
3047        #Check values
3048        assert allclose(elevation, ref_elevation)
3049
3050        #Cleanup
3051        fid.close()
3052
3053        os.remove(root + '.dem')
3054        os.remove(root + '_100.dem')
3055
3056    def test_decimate_dem_NODATA(self):
3057        """Test decimation of dem file that includes NODATA values
3058        """
3059
3060        import os
3061        from Numeric import ones, allclose, Float, arange, reshape
3062        from Scientific.IO.NetCDF import NetCDFFile
3063
3064        #Write test dem file
3065        root = 'decdemtest'
3066
3067        filename = root + '.dem'
3068        fid = NetCDFFile(filename, 'w')
3069
3070        fid.institution = 'Geoscience Australia'
3071        fid.description = 'NetCDF DEM format for compact and portable ' +\
3072                          'storage of spatial point data'
3073
3074        nrows = 15
3075        ncols = 18
3076        NODATA_value = -9999
3077
3078        fid.ncols = ncols
3079        fid.nrows = nrows
3080        fid.xllcorner = 2000.5
3081        fid.yllcorner = 3000.5
3082        fid.cellsize = 25
3083        fid.NODATA_value = NODATA_value
3084
3085        fid.zone = 56
3086        fid.false_easting = 0.0
3087        fid.false_northing = 0.0
3088        fid.projection = 'UTM'
3089        fid.datum = 'WGS84'
3090        fid.units = 'METERS'
3091
3092        fid.createDimension('number_of_points', nrows*ncols)
3093
3094        fid.createVariable('elevation', Float, ('number_of_points',))
3095
3096        elevation = fid.variables['elevation']
3097
3098        #generate initial elevation values
3099        elevation_tmp = (arange(nrows*ncols))
3100        #add some NODATA values
3101        elevation_tmp[0]   = NODATA_value
3102        elevation_tmp[95]  = NODATA_value
3103        elevation_tmp[188] = NODATA_value
3104        elevation_tmp[189] = NODATA_value
3105        elevation_tmp[190] = NODATA_value
3106        elevation_tmp[209] = NODATA_value
3107        elevation_tmp[252] = NODATA_value
3108
3109        elevation[:] = elevation_tmp
3110
3111        fid.close()
3112
3113        #generate the elevation values expected in the decimated file
3114        ref_elevation = [NODATA_value,
3115                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3116                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3117                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3118                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3119                         NODATA_value,
3120                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3121                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3122                         (144+145+146+162+163+164+180+181+182) / 9.0,
3123                         (148+149+150+166+167+168+184+185+186) / 9.0,
3124                         NODATA_value,
3125                         (156+157+158+174+175+176+192+193+194) / 9.0,
3126                         NODATA_value,
3127                         (220+221+222+238+239+240+256+257+258) / 9.0,
3128                         (224+225+226+242+243+244+260+261+262) / 9.0,
3129                         (228+229+230+246+247+248+264+265+266) / 9.0]
3130
3131        #generate a stencil for computing the decimated values
3132        stencil = ones((3,3), Float) / 9.0
3133
3134        decimate_dem(root, stencil=stencil, cellsize_new=100)
3135
3136        #Open decimated NetCDF file
3137        fid = NetCDFFile(root + '_100.dem', 'r')
3138
3139        # Get decimated elevation
3140        elevation = fid.variables['elevation']
3141
3142        #Check values
3143        assert allclose(elevation, ref_elevation)
3144
3145        #Cleanup
3146        fid.close()
3147
3148        os.remove(root + '.dem')
3149        os.remove(root + '_100.dem')
3150
3151    def xxxtestz_sww2ers_real(self):
3152        """Test that sww information can be converted correctly to asc/prj
3153        format readable by e.g. ArcView
3154        """
3155
3156        import time, os
3157        from Numeric import array, zeros, allclose, Float, concatenate
3158        from Scientific.IO.NetCDF import NetCDFFile
3159
3160        # the memory optimised least squares
3161        #  cellsize = 20,   # this one seems to hang
3162        #  cellsize = 200000, # Ran 1 test in 269.703s
3163                                #Ran 1 test in 267.344s
3164        #  cellsize = 20000,  # Ran 1 test in 460.922s
3165        #  cellsize = 2000   #Ran 1 test in 5340.250s
3166        #  cellsize = 200   #this one seems to hang, building matirx A
3167
3168        # not optimised
3169        # seems to hang
3170        #  cellsize = 2000   # Ran 1 test in 5334.563s
3171        #Export to ascii/prj files
3172        sww2dem('karratha_100m',
3173                quantity = 'depth',
3174                cellsize = 200000,
3175                verbose = True)
3176
3177    def test_read_asc(self):
3178        """Test conversion from dem in ascii format to native NetCDF xya format
3179        """
3180
3181        import time, os
3182        from Numeric import array, zeros, allclose, Float, concatenate
3183        from Scientific.IO.NetCDF import NetCDFFile
3184
3185        import data_manager
3186        #Write test asc file
3187        filename = tempfile.mktemp(".000")
3188        fid = open(filename, 'w')
3189        fid.write("""ncols         7
3190nrows         4
3191xllcorner     2000.5
3192yllcorner     3000.5
3193cellsize      25
3194NODATA_value  -9999
3195    97.921    99.285   125.588   180.830   258.645   342.872   415.836
3196   473.157   514.391   553.893   607.120   678.125   777.283   883.038
3197   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
3198   502.645   516.230   504.739   450.604   388.500   338.097   514.980
3199""")
3200        fid.close()
3201        bath_metadata, grid = \
3202                   data_manager._read_asc(filename, verbose=False)
3203        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
3204        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
3205        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
3206        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
3207        self.failUnless(grid[0][0]  == 97.921,  'Failed')
3208        self.failUnless(grid[3][6]  == 514.980,  'Failed')
3209
3210        os.remove(filename)
3211
3212    def test_asc_csiro2sww(self):
3213        import tempfile
3214
3215        bath_dir = tempfile.mkdtemp()
3216        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3217        #bath_dir = 'bath_data_manager_test'
3218        #print "os.getcwd( )",os.getcwd( )
3219        elevation_dir =  tempfile.mkdtemp()
3220        #elevation_dir = 'elev_expanded'
3221        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3222        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3223
3224        fid = open(bath_dir_filename, 'w')
3225        fid.write(""" ncols             3
3226 nrows             2
3227 xllcorner    148.00000
3228 yllcorner    -38.00000
3229 cellsize       0.25
3230 nodata_value   -9999.0
3231    9000.000    -1000.000    3000.0
3232   -1000.000    9000.000  -1000.000
3233""")
3234        fid.close()
3235
3236        fid = open(elevation_dir_filename1, 'w')
3237        fid.write(""" ncols             3
3238 nrows             2
3239 xllcorner    148.00000
3240 yllcorner    -38.00000
3241 cellsize       0.25
3242 nodata_value   -9999.0
3243    9000.000    0.000    3000.0
3244     0.000     9000.000     0.000
3245""")
3246        fid.close()
3247
3248        fid = open(elevation_dir_filename2, 'w')
3249        fid.write(""" ncols             3
3250 nrows             2
3251 xllcorner    148.00000
3252 yllcorner    -38.00000
3253 cellsize       0.25
3254 nodata_value   -9999.0
3255    9000.000    4000.000    4000.0
3256    4000.000    9000.000    4000.000
3257""")
3258        fid.close()
3259
3260        ucur_dir =  tempfile.mkdtemp()
3261        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3262        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3263
3264        fid = open(ucur_dir_filename1, 'w')
3265        fid.write(""" ncols             3
3266 nrows             2
3267 xllcorner    148.00000
3268 yllcorner    -38.00000
3269 cellsize       0.25
3270 nodata_value   -9999.0
3271    90.000    60.000    30.0
3272    10.000    10.000    10.000
3273""")
3274        fid.close()
3275        fid = open(ucur_dir_filename2, 'w')
3276        fid.write(""" ncols             3
3277 nrows             2
3278 xllcorner    148.00000
3279 yllcorner    -38.00000
3280 cellsize       0.25
3281 nodata_value   -9999.0
3282    90.000    60.000    30.0
3283    10.000    10.000    10.000
3284""")
3285        fid.close()
3286
3287        vcur_dir =  tempfile.mkdtemp()
3288        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3289        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3290
3291        fid = open(vcur_dir_filename1, 'w')
3292        fid.write(""" ncols             3
3293 nrows             2
3294 xllcorner    148.00000
3295 yllcorner    -38.00000
3296 cellsize       0.25
3297 nodata_value   -9999.0
3298    90.000    60.000    30.0
3299    10.000    10.000    10.000
3300""")
3301        fid.close()
3302        fid = open(vcur_dir_filename2, 'w')
3303        fid.write(""" ncols             3
3304 nrows             2
3305 xllcorner    148.00000
3306 yllcorner    -38.00000
3307 cellsize       0.25
3308 nodata_value   -9999.0
3309    90.000    60.000    30.0
3310    10.000    10.000    10.000
3311""")
3312        fid.close()
3313
3314        sww_file = 'a_test.sww'
3315        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
3316
3317        # check the sww file
3318
3319        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3320        x = fid.variables['x'][:]
3321        y = fid.variables['y'][:]
3322        z = fid.variables['z'][:]
3323        stage = fid.variables['stage'][:]
3324        xmomentum = fid.variables['xmomentum'][:]
3325        geo_ref = Geo_reference(NetCDFObject=fid)
3326        #print "geo_ref",geo_ref
3327        x_ref = geo_ref.get_xllcorner()
3328        y_ref = geo_ref.get_yllcorner()
3329        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3330        assert allclose(x_ref, 587798.418) # (-38, 148)
3331        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3332
3333        #Zone:   55
3334        #Easting:  588095.674  Northing: 5821451.722
3335        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3336        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3337
3338        #Zone:   55
3339        #Easting:  632145.632  Northing: 5820863.269
3340        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3341        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3342
3343        #Zone:   55
3344        #Easting:  609748.788  Northing: 5793447.860
3345        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3346        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3347
3348        assert allclose(z[0],9000.0 )
3349        assert allclose(stage[0][1],0.0 )
3350
3351        #(4000+1000)*60
3352        assert allclose(xmomentum[1][1],300000.0 )
3353
3354
3355        fid.close()
3356
3357        #tidy up
3358        os.remove(bath_dir_filename)
3359        os.rmdir(bath_dir)
3360
3361        os.remove(elevation_dir_filename1)
3362        os.remove(elevation_dir_filename2)
3363        os.rmdir(elevation_dir)
3364
3365        os.remove(ucur_dir_filename1)
3366        os.remove(ucur_dir_filename2)
3367        os.rmdir(ucur_dir)
3368
3369        os.remove(vcur_dir_filename1)
3370        os.remove(vcur_dir_filename2)
3371        os.rmdir(vcur_dir)
3372
3373
3374        # remove sww file
3375        os.remove(sww_file)
3376
3377    def test_asc_csiro2sww2(self):
3378        import tempfile
3379
3380        bath_dir = tempfile.mkdtemp()
3381        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3382        #bath_dir = 'bath_data_manager_test'
3383        #print "os.getcwd( )",os.getcwd( )
3384        elevation_dir =  tempfile.mkdtemp()
3385        #elevation_dir = 'elev_expanded'
3386        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3387        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3388
3389        fid = open(bath_dir_filename, 'w')
3390        fid.write(""" ncols             3
3391 nrows             2
3392 xllcorner    148.00000
3393 yllcorner    -38.00000
3394 cellsize       0.25
3395 nodata_value   -9999.0
3396    9000.000    -1000.000    3000.0
3397   -1000.000    9000.000  -1000.000
3398""")
3399        fid.close()
3400
3401        fid = open(elevation_dir_filename1, 'w')
3402        fid.write(""" ncols             3
3403 nrows             2
3404 xllcorner    148.00000
3405 yllcorner    -38.00000
3406 cellsize       0.25
3407 nodata_value   -9999.0
3408    9000.000    0.000    3000.0
3409     0.000     -9999.000     -9999.000
3410""")
3411        fid.close()
3412
3413        fid = open(elevation_dir_filename2, 'w')
3414        fid.write(""" ncols             3
3415 nrows             2
3416 xllcorner    148.00000
3417 yllcorner    -38.00000
3418 cellsize       0.25
3419 nodata_value   -9999.0
3420    9000.000    4000.000    4000.0
3421    4000.000    9000.000    4000.000
3422""")
3423        fid.close()
3424
3425        ucur_dir =  tempfile.mkdtemp()
3426        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3427        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3428
3429        fid = open(ucur_dir_filename1, 'w')
3430        fid.write(""" ncols             3
3431 nrows             2
3432 xllcorner    148.00000
3433 yllcorner    -38.00000
3434 cellsize       0.25
3435 nodata_value   -9999.0
3436    90.000    60.000    30.0
3437    10.000    10.000    10.000
3438""")
3439        fid.close()
3440        fid = open(ucur_dir_filename2, 'w')
3441        fid.write(""" ncols             3
3442 nrows             2
3443 xllcorner    148.00000
3444 yllcorner    -38.00000
3445 cellsize       0.25
3446 nodata_value   -9999.0
3447    90.000    60.000    30.0
3448    10.000    10.000    10.000
3449""")
3450        fid.close()
3451
3452        vcur_dir =  tempfile.mkdtemp()
3453        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3454        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3455
3456        fid = open(vcur_dir_filename1, 'w')
3457        fid.write(""" ncols             3
3458 nrows             2
3459 xllcorner    148.00000
3460 yllcorner    -38.00000
3461 cellsize       0.25
3462 nodata_value   -9999.0
3463    90.000    60.000    30.0
3464    10.000    10.000    10.000
3465""")
3466        fid.close()
3467        fid = open(vcur_dir_filename2, 'w')
3468        fid.write(""" ncols             3
3469 nrows             2
3470 xllcorner    148.00000
3471 yllcorner    -38.00000
3472 cellsize       0.25
3473 nodata_value   -9999.0
3474    90.000    60.000    30.0
3475    10.000    10.000    10.000
3476""")
3477        fid.close()
3478
3479        try:
3480            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3481                          vcur_dir, sww_file)
3482        except:
3483            #tidy up
3484            os.remove(bath_dir_filename)
3485            os.rmdir(bath_dir)
3486
3487            os.remove(elevation_dir_filename1)
3488            os.remove(elevation_dir_filename2)
3489            os.rmdir(elevation_dir)
3490
3491            os.remove(ucur_dir_filename1)
3492            os.remove(ucur_dir_filename2)
3493            os.rmdir(ucur_dir)
3494
3495            os.remove(vcur_dir_filename1)
3496            os.remove(vcur_dir_filename2)
3497            os.rmdir(vcur_dir)
3498        else:
3499            #tidy up
3500            os.remove(bath_dir_filename)
3501            os.rmdir(bath_dir)
3502
3503            os.remove(elevation_dir_filename1)
3504            os.remove(elevation_dir_filename2)
3505            os.rmdir(elevation_dir)
3506            raise 'Should raise exception'
3507
3508            os.remove(ucur_dir_filename1)
3509            os.remove(ucur_dir_filename2)
3510            os.rmdir(ucur_dir)
3511
3512            os.remove(vcur_dir_filename1)
3513            os.remove(vcur_dir_filename2)
3514            os.rmdir(vcur_dir)
3515
3516
3517
3518    def test_asc_csiro2sww3(self):
3519        import tempfile
3520
3521        bath_dir = tempfile.mkdtemp()
3522        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3523        #bath_dir = 'bath_data_manager_test'
3524        #print "os.getcwd( )",os.getcwd( )
3525        elevation_dir =  tempfile.mkdtemp()
3526        #elevation_dir = 'elev_expanded'
3527        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3528        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3529
3530        fid = open(bath_dir_filename, 'w')
3531        fid.write(""" ncols             3
3532 nrows             2
3533 xllcorner    148.00000
3534 yllcorner    -38.00000
3535 cellsize       0.25
3536 nodata_value   -9999.0
3537    9000.000    -1000.000    3000.0
3538   -1000.000    9000.000  -1000.000
3539""")
3540        fid.close()
3541
3542        fid = open(elevation_dir_filename1, 'w')
3543        fid.write(""" ncols             3
3544 nrows             2
3545 xllcorner    148.00000
3546 yllcorner    -38.00000
3547 cellsize       0.25
3548 nodata_value   -9999.0
3549    9000.000    0.000    3000.0
3550     0.000     -9999.000     -9999.000
3551""")
3552        fid.close()
3553
3554        fid = open(elevation_dir_filename2, 'w')
3555        fid.write(""" ncols             3
3556 nrows             2
3557 xllcorner    148.00000
3558 yllcorner    -38.00000
3559 cellsize       0.25
3560 nodata_value   -9999.0
3561    9000.000    4000.000    4000.0
3562    4000.000    9000.000    4000.000
3563""")
3564        fid.close()
3565
3566        ucur_dir =  tempfile.mkdtemp()
3567        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3568        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3569
3570        fid = open(ucur_dir_filename1, 'w')
3571        fid.write(""" ncols             3
3572 nrows             2
3573 xllcorner    148.00000
3574 yllcorner    -38.00000
3575 cellsize       0.25
3576 nodata_value   -9999.0
3577    90.000    60.000    30.0
3578    10.000    10.000    10.000
3579""")
3580        fid.close()
3581        fid = open(ucur_dir_filename2, 'w')
3582        fid.write(""" ncols             3
3583 nrows             2
3584 xllcorner    148.00000
3585 yllcorner    -38.00000
3586 cellsize       0.25
3587 nodata_value   -9999.0
3588    90.000    60.000    30.0
3589    10.000    10.000    10.000
3590""")
3591        fid.close()
3592
3593        vcur_dir =  tempfile.mkdtemp()
3594        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3595        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3596
3597        fid = open(vcur_dir_filename1, 'w')
3598        fid.write(""" ncols             3
3599 nrows             2
3600 xllcorner    148.00000
3601 yllcorner    -38.00000
3602 cellsize       0.25
3603 nodata_value   -9999.0
3604    90.000    60.000    30.0
3605    10.000    10.000    10.000
3606""")
3607        fid.close()
3608        fid = open(vcur_dir_filename2, 'w')
3609        fid.write(""" ncols             3
3610 nrows             2
3611 xllcorner    148.00000
3612 yllcorner    -38.00000
3613 cellsize       0.25
3614 nodata_value   -9999.0
3615    90.000    60.000    30.0
3616    10.000    10.000    10.000
3617""")
3618        fid.close()
3619
3620        sww_file = 'a_test.sww'
3621        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3622                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3623                      mean_stage = 100)
3624
3625        # check the sww file
3626
3627        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3628        x = fid.variables['x'][:]
3629        y = fid.variables['y'][:]
3630        z = fid.variables['z'][:]
3631        stage = fid.variables['stage'][:]
3632        xmomentum = fid.variables['xmomentum'][:]
3633        geo_ref = Geo_reference(NetCDFObject=fid)
3634        #print "geo_ref",geo_ref
3635        x_ref = geo_ref.get_xllcorner()
3636        y_ref = geo_ref.get_yllcorner()
3637        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3638        assert allclose(x_ref, 587798.418) # (-38, 148)
3639        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3640
3641        #Zone:   55
3642        #Easting:  588095.674  Northing: 5821451.722
3643        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3644        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3645
3646        #Zone:   55
3647        #Easting:  632145.632  Northing: 5820863.269
3648        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3649        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3650
3651        #Zone:   55
3652        #Easting:  609748.788  Northing: 5793447.860
3653        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3654        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3655
3656        assert allclose(z[0],9000.0 )
3657        assert allclose(stage[0][4],100.0 )
3658        assert allclose(stage[0][5],100.0 )
3659
3660        #(100.0 - 9000)*10
3661        assert allclose(xmomentum[0][4], -89000.0 )
3662
3663        #(100.0 - -1000.000)*10
3664        assert allclose(xmomentum[0][5], 11000.0 )
3665
3666        fid.close()
3667
3668        #tidy up
3669        os.remove(bath_dir_filename)
3670        os.rmdir(bath_dir)
3671
3672        os.remove(elevation_dir_filename1)
3673        os.remove(elevation_dir_filename2)
3674        os.rmdir(elevation_dir)
3675
3676        os.remove(ucur_dir_filename1)
3677        os.remove(ucur_dir_filename2)
3678        os.rmdir(ucur_dir)
3679
3680        os.remove(vcur_dir_filename1)
3681        os.remove(vcur_dir_filename2)
3682        os.rmdir(vcur_dir)
3683
3684        # remove sww file
3685        os.remove(sww_file)
3686
3687
3688    def test_asc_csiro2sww4(self):
3689        """
3690        Test specifying the extent
3691        """
3692
3693        import tempfile
3694
3695        bath_dir = tempfile.mkdtemp()
3696        bath_dir_filename = bath_dir + os.sep +'ba19940524.000'
3697        #bath_dir = 'bath_data_manager_test'
3698        #print "os.getcwd( )",os.getcwd( )
3699        elevation_dir =  tempfile.mkdtemp()
3700        #elevation_dir = 'elev_expanded'
3701        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000'
3702        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001'
3703
3704        fid = open(bath_dir_filename, 'w')
3705        fid.write(""" ncols             4
3706 nrows             4
3707 xllcorner    148.00000
3708 yllcorner    -38.00000
3709 cellsize       0.25
3710 nodata_value   -9999.0
3711   -9000.000    -1000.000   -3000.0 -2000.000
3712   -1000.000    9000.000  -1000.000 -3000.000
3713   -4000.000    6000.000   2000.000 -5000.000
3714   -9000.000    -1000.000   -3000.0 -2000.000
3715""")
3716        fid.close()
3717
3718        fid = open(elevation_dir_filename1, 'w')
3719        fid.write(""" ncols             4
3720 nrows             4
3721 xllcorner    148.00000
3722 yllcorner    -38.00000
3723 cellsize       0.25
3724 nodata_value   -9999.0
3725   -900.000    -100.000   -300.0 -200.000
3726   -100.000    900.000  -100.000 -300.000
3727   -400.000    600.000   200.000 -500.000
3728   -900.000    -100.000   -300.0 -200.000
3729""")
3730        fid.close()
3731
3732        fid = open(elevation_dir_filename2, 'w')
3733        fid.write(""" ncols             4
3734 nrows             4
3735 xllcorner    148.00000
3736 yllcorner    -38.00000
3737 cellsize       0.25
3738 nodata_value   -9999.0
3739   -990.000    -110.000   -330.0 -220.000
3740   -110.000    990.000  -110.000 -330.000
3741   -440.000    660.000   220.000 -550.000
3742   -990.000    -110.000   -330.0 -220.000
3743""")
3744        fid.close()
3745
3746        ucur_dir =  tempfile.mkdtemp()
3747        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3748        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3749
3750        fid = open(ucur_dir_filename1, 'w')
3751        fid.write(""" ncols             4
3752 nrows             4
3753 xllcorner    148.00000
3754 yllcorner    -38.00000
3755 cellsize       0.25
3756 nodata_value   -9999.0
3757   -90.000    -10.000   -30.0 -20.000
3758   -10.000    90.000  -10.000 -30.000
3759   -40.000    60.000   20.000 -50.000
3760   -90.000    -10.000   -30.0 -20.000
3761""")
3762        fid.close()
3763        fid = open(ucur_dir_filename2, 'w')
3764        fid.write(""" ncols             4
3765 nrows             4
3766 xllcorner    148.00000
3767 yllcorner    -38.00000
3768 cellsize       0.25
3769 nodata_value   -9999.0
3770   -90.000    -10.000   -30.0 -20.000
3771   -10.000    99.000  -11.000 -30.000
3772   -40.000    66.000   22.000 -50.000
3773   -90.000    -10.000   -30.0 -20.000
3774""")
3775        fid.close()
3776
3777        vcur_dir =  tempfile.mkdtemp()
3778        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3779        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3780
3781        fid = open(vcur_dir_filename1, 'w')
3782        fid.write(""" ncols             4
3783 nrows             4
3784 xllcorner    148.00000
3785 yllcorner    -38.00000
3786 cellsize       0.25
3787 nodata_value   -9999.0
3788   -90.000    -10.000   -30.0 -20.000
3789   -10.000    80.000  -20.000 -30.000
3790   -40.000    50.000   10.000 -50.000
3791   -90.000    -10.000   -30.0 -20.000
3792""")
3793        fid.close()
3794        fid = open(vcur_dir_filename2, 'w')
3795        fid.write(""" ncols             4
3796 nrows             4
3797 xllcorner    148.00000
3798 yllcorner    -38.00000
3799 cellsize       0.25
3800 nodata_value   -9999.0
3801   -90.000    -10.000   -30.0 -20.000
3802   -10.000    88.000  -22.000 -30.000
3803   -40.000    55.000   11.000 -50.000
3804   -90.000    -10.000   -30.0 -20.000
3805""")
3806        fid.close()
3807
3808        sww_file = tempfile.mktemp(".sww")
3809        #sww_file = 'a_test.sww'
3810        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3811                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3812                      mean_stage = 100,
3813                       minlat = -37.6, maxlat = -37.6,
3814                  minlon = 148.3, maxlon = 148.3
3815                      #,verbose = True
3816                      )
3817
3818        # check the sww file
3819
3820        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3821        x = fid.variables['x'][:]
3822        y = fid.variables['y'][:]
3823        z = fid.variables['z'][:]
3824        stage = fid.variables['stage'][:]
3825        xmomentum = fid.variables['xmomentum'][:]
3826        ymomentum = fid.variables['ymomentum'][:]
3827        geo_ref = Geo_reference(NetCDFObject=fid)
3828        #print "geo_ref",geo_ref
3829        x_ref = geo_ref.get_xllcorner()
3830        y_ref = geo_ref.get_yllcorner()
3831        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3832
3833        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3834        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3835        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
3836
3837        #Easting:  632145.632  Northing: 5820863.269
3838        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3839
3840        #print "x",x
3841        #print "y",y
3842        self.failUnless(len(x) == 4,'failed') # 2*2
3843        self.failUnless(len(x) == 4,'failed') # 2*2
3844
3845        #Zone:   55
3846        #Easting:  632145.632  Northing: 5820863.269
3847        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3848        # magic number - y is close enough for me.
3849        assert allclose(x[3], 632145.63 - x_ref)
3850        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
3851
3852        assert allclose(z[0],9000.0 ) #z is elevation info
3853        #print "z",z
3854        # 2 time steps, 4 points
3855        self.failUnless(xmomentum.shape == (2,4), 'failed')
3856        self.failUnless(ymomentum.shape == (2,4), 'failed')
3857
3858        #(100.0 - -1000.000)*10
3859        #assert allclose(xmomentum[0][5], 11000.0 )
3860
3861        fid.close()
3862
3863        # is the sww file readable?
3864        #Lets see if we can convert it to a dem!
3865        #print "sww_file",sww_file
3866        #dem_file = tempfile.mktemp(".dem")
3867        domain = sww2domain(sww_file) ###, dem_file)
3868        domain.check_integrity()
3869
3870        #tidy up
3871        os.remove(bath_dir_filename)
3872        os.rmdir(bath_dir)
3873
3874        os.remove(elevation_dir_filename1)
3875        os.remove(elevation_dir_filename2)
3876        os.rmdir(elevation_dir)
3877
3878        os.remove(ucur_dir_filename1)
3879        os.remove(ucur_dir_filename2)
3880        os.rmdir(ucur_dir)
3881
3882        os.remove(vcur_dir_filename1)
3883        os.remove(vcur_dir_filename2)
3884        os.rmdir(vcur_dir)
3885
3886
3887
3888
3889        # remove sww file
3890        os.remove(sww_file)
3891
3892        # remove dem file
3893        #os.remove(dem_file)
3894
3895    def test_get_min_max_indexes(self):
3896        latitudes = [3,2,1,0]
3897        longitudes = [0,10,20,30]
3898
3899        # k - lat
3900        # l - lon
3901        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3902            latitudes,longitudes,
3903            -10,4,-10,31)
3904
3905        #print "kmin",kmin;print "kmax",kmax
3906        #print "lmin",lmin;print "lmax",lmax
3907        latitudes_new = latitudes[kmin:kmax]
3908        longitudes_news = longitudes[lmin:lmax]
3909        #print "latitudes_new", latitudes_new
3910        #print "longitudes_news",longitudes_news
3911        self.failUnless(latitudes == latitudes_new and \
3912                        longitudes == longitudes_news,
3913                         'failed')
3914
3915        ## 2nd test
3916        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3917            latitudes,longitudes,
3918            0.5,2.5,5,25)
3919        #print "kmin",kmin;print "kmax",kmax
3920        #print "lmin",lmin;print "lmax",lmax
3921        latitudes_new = latitudes[kmin:kmax]
3922        longitudes_news = longitudes[lmin:lmax]
3923        #print "latitudes_new", latitudes_new
3924        #print "longitudes_news",longitudes_news
3925
3926        self.failUnless(latitudes == latitudes_new and \
3927                        longitudes == longitudes_news,
3928                         'failed')
3929
3930        ## 3rd test
3931        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
3932                                                                    longitudes,
3933                                                      1.1,1.9,12,17)
3934        #print "kmin",kmin;print "kmax",kmax
3935        #print "lmin",lmin;print "lmax",lmax
3936        latitudes_new = latitudes[kmin:kmax]
3937        longitudes_news = longitudes[lmin:lmax]
3938        #print "latitudes_new", latitudes_new
3939        #print "longitudes_news",longitudes_news
3940
3941        self.failUnless(latitudes_new == [2, 1] and \
3942                        longitudes_news == [10, 20],
3943                         'failed')
3944
3945
3946        ## 4th test
3947        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3948            latitudes,longitudes,
3949                                                      -0.1,1.9,-2,17)
3950        #print "kmin",kmin;print "kmax",kmax
3951        #print "lmin",lmin;print "lmax",lmax
3952        latitudes_new = latitudes[kmin:kmax]
3953        longitudes_news = longitudes[lmin:lmax]
3954        #print "latitudes_new", latitudes_new
3955        #print "longitudes_news",longitudes_news
3956
3957        self.failUnless(latitudes_new == [2, 1, 0] and \
3958                        longitudes_news == [0, 10, 20],
3959                         'failed')
3960        ## 5th test
3961        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3962            latitudes,longitudes,
3963            0.1,1.9,2,17)
3964        #print "kmin",kmin;print "kmax",kmax
3965        #print "lmin",lmin;print "lmax",lmax
3966        latitudes_new = latitudes[kmin:kmax]
3967        longitudes_news = longitudes[lmin:lmax]
3968        #print "latitudes_new", latitudes_new
3969        #print "longitudes_news",longitudes_news
3970
3971        self.failUnless(latitudes_new == [2, 1, 0] and \
3972                        longitudes_news == [0, 10, 20],
3973                         'failed')
3974
3975        ## 6th test
3976
3977        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3978            latitudes,longitudes,
3979            1.5,4,18,32)
3980        #print "kmin",kmin;print "kmax",kmax
3981        #print "lmin",lmin;print "lmax",lmax
3982        latitudes_new = latitudes[kmin:kmax]
3983        longitudes_news = longitudes[lmin:lmax]
3984        #print "latitudes_new", latitudes_new
3985        #print "longitudes_news",longitudes_news
3986
3987        self.failUnless(latitudes_new == [3, 2, 1] and \
3988                        longitudes_news == [10, 20, 30],
3989                         'failed')
3990
3991
3992        ## 7th test
3993        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3994        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3995            latitudes,longitudes,
3996            1.5,1.5,15,15)
3997        #print "kmin",kmin;print "kmax",kmax
3998        #print "lmin",lmin;print "lmax",lmax
3999        latitudes_new = latitudes[kmin:kmax]
4000        longitudes_news = longitudes[lmin:lmax]
4001        m2d = m2d[kmin:kmax,lmin:lmax]
4002        #print "m2d", m2d
4003        #print "latitudes_new", latitudes_new
4004        #print "longitudes_news",longitudes_news
4005
4006        self.failUnless(latitudes_new == [2, 1] and \
4007                        longitudes_news == [10, 20],
4008                         'failed')
4009
4010        self.failUnless(m2d == [[5,6],[9,10]],
4011                         'failed')
4012
4013    def test_get_min_max_indexes2(self):
4014        latitudes = [-30,-35,-40,-45]
4015        longitudes = [148,149,150,151]
4016
4017        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
4018
4019        # k - lat
4020        # l - lon
4021        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4022            latitudes,longitudes,
4023            -37,-27,147,149.5)
4024
4025        #print "kmin",kmin;print "kmax",kmax
4026        #print "lmin",lmin;print "lmax",lmax
4027        #print "m2d", m2d
4028        #print "latitudes", latitudes
4029        #print "longitudes",longitudes
4030        #print "latitudes[kmax]", latitudes[kmax]
4031        latitudes_new = latitudes[kmin:kmax]
4032        longitudes_new = longitudes[lmin:lmax]
4033        m2d = m2d[kmin:kmax,lmin:lmax]
4034        #print "m2d", m2d
4035        #print "latitudes_new", latitudes_new
4036        #print "longitudes_new",longitudes_new
4037
4038        self.failUnless(latitudes_new == [-30, -35, -40] and \
4039                        longitudes_new == [148, 149,150],
4040                         'failed')
4041        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
4042                         'failed')
4043
4044    def test_get_min_max_indexes3(self):
4045        latitudes = [-30,-35,-40,-45,-50,-55,-60]
4046        longitudes = [148,149,150,151]
4047
4048        # k - lat
4049        # l - lon
4050        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4051            latitudes,longitudes,
4052            -43,-37,148.5,149.5)
4053
4054
4055        #print "kmin",kmin;print "kmax",kmax
4056        #print "lmin",lmin;print "lmax",lmax
4057        #print "latitudes", latitudes
4058        #print "longitudes",longitudes
4059        latitudes_new = latitudes[kmin:kmax]
4060        longitudes_news = longitudes[lmin:lmax]
4061        #print "latitudes_new", latitudes_new
4062        #print "longitudes_news",longitudes_news
4063
4064        self.failUnless(latitudes_new == [-35, -40, -45] and \
4065                        longitudes_news == [148, 149,150],
4066                         'failed')
4067
4068    def test_get_min_max_indexes4(self):
4069        latitudes = [-30,-35,-40,-45,-50,-55,-60]
4070        longitudes = [148,149,150,151]
4071
4072        # k - lat
4073        # l - lon
4074        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
4075            latitudes,longitudes)
4076
4077        #print "kmin",kmin;print "kmax",kmax
4078        #print "lmin",lmin;print "lmax",lmax
4079        #print "latitudes", latitudes
4080        #print "longitudes",longitudes
4081        latitudes_new = latitudes[kmin:kmax]
4082        longitudes_news = longitudes[lmin:lmax]
4083        #print "latitudes_new", latitudes_new
4084        #print "longitudes_news",longitudes_news
4085
4086        self.failUnless(latitudes_new == latitudes  and \
4087                        longitudes_news == longitudes,
4088                         'failed')
4089
4090    def test_tsh2sww(self):
4091        import os
4092        import tempfile
4093
4094        tsh_file = tempfile.mktemp(".tsh")
4095        file = open(tsh_file,"w")
4096        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
40970 0.0 0.0 0.0 0.0 0.01 \n \
40981 1.0 0.0 10.0 10.0 0.02  \n \
40992 0.0 1.0 0.0 10.0 0.03  \n \
41003 0.5 0.25 8.0 12.0 0.04  \n \
4101# Vert att title  \n \
4102elevation  \n \
4103stage  \n \
4104friction  \n \
41052 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
41060 0 3 2 -1  -1  1 dsg\n\
41071 0 1 3 -1  0 -1   ole nielsen\n\
41084 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
41090 1 0 2 \n\
41101 0 2 3 \n\
41112 2 3 \n\
41123 3 1 1 \n\
41133 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
41140 216.0 -86.0 \n \
41151 160.0 -167.0 \n \
41162 114.0 -91.0 \n \
41173 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
41180 0 1 0 \n \
41191 1 2 0 \n \
41202 2 0 0 \n \
41210 # <x> <y> ...Mesh Holes... \n \
41220 # <x> <y> <attribute>...Mesh Regions... \n \
41230 # <x> <y> <attribute>...Mesh Regions, area... \n\
4124#Geo reference \n \
412556 \n \
4126140 \n \
4127120 \n")
4128        file.close()
4129
4130        #sww_file = tempfile.mktemp(".sww")
4131        #print "sww_file",sww_file
4132        #print "sww_file",tsh_file
4133        tsh2sww(tsh_file)
4134
4135        os.remove(tsh_file)
4136        os.remove(tsh_file[:-4] + '.sww')
4137#-------------------------------------------------------------
4138if __name__ == "__main__":
4139    #suite = unittest.makeSuite(Test_Data_Manager,'test_tsh2sww')
4140    suite = unittest.makeSuite(Test_Data_Manager,'test')
4141    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
4142    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')
4143    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
4144    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
4145    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
4146    runner = unittest.TextTestRunner()
4147    runner.run(suite)
4148
4149
4150
4151
4152    def test_dem2pts(self):
4153        """Test conversion from dem in ascii format to native NetCDF xya format
4154        """
4155
4156        import time, os
4157        from Numeric import array, zeros, allclose, Float, concatenate
4158        from Scientific.IO.NetCDF import NetCDFFile
4159
4160        #Write test asc file
4161        root = 'demtest'
4162
4163        filename = root+'.asc'
4164        fid = open(filename, 'w')
4165        fid.write("""ncols         5
4166nrows         6
4167xllcorner     2000.5
4168yllcorner     3000.5
4169cellsize      25
4170NODATA_value  -9999
4171""")
4172        #Create linear function
4173
4174        ref_points = []
4175        ref_elevation = []
4176        for i in range(6):
4177            y = (6-i)*25.0
4178            for j in range(5):
4179                x = j*25.0
4180                z = x+2*y
4181
4182                ref_points.append( [x,y] )
4183                ref_elevation.append(z)
4184                fid.write('%f ' %z)
4185            fid.write('\n')
4186
4187        fid.close()
4188
4189        #Write prj file with metadata
4190        metafilename = root+'.prj'
4191        fid = open(metafilename, 'w')
4192
4193
4194        fid.write("""Projection UTM
4195Zone 56
4196Datum WGS84
4197Zunits NO
4198Units METERS
4199Spheroid WGS84
4200Xshift 0.0000000000
4201Yshift 10000000.0000000000
4202Parameters
4203""")
4204        fid.close()
4205
4206        #Convert to NetCDF pts
4207        convert_dem_from_ascii2netcdf(root)
4208        dem2pts(root)
4209
4210        #Check contents
4211        #Get NetCDF
4212        fid = NetCDFFile(root+'.pts', 'r')
4213
4214        # Get the variables
4215        #print fid.variables.keys()
4216        points = fid.variables['points']
4217        elevation = fid.variables['elevation']
4218
4219        #Check values
4220
4221        #print points[:]
4222        #print ref_points
4223        assert allclose(points, ref_points)
4224
4225        #print attributes[:]
4226        #print ref_elevation
4227        assert allclose(elevation, ref_elevation)
4228
4229        #Cleanup
4230        fid.close()
4231
4232
4233        os.remove(root + '.pts')
4234        os.remove(root + '.dem')
4235        os.remove(root + '.asc')
4236        os.remove(root + '.prj')
4237
4238
4239
4240    def test_dem2pts_bounding_box(self):
4241        """Test conversion from dem in ascii format to native NetCDF xya format
4242        """
4243
4244        import time, os
4245        from Numeric import array, zeros, allclose, Float, concatenate
4246        from Scientific.IO.NetCDF import NetCDFFile
4247
4248        #Write test asc file
4249        root = 'demtest'
4250
4251        filename = root+'.asc'
4252        fid = open(filename, 'w')
4253        fid.write("""ncols         5
4254nrows         6
4255xllcorner     2000.5
4256yllcorner     3000.5
4257cellsize      25
4258NODATA_value  -9999
4259""")
4260        #Create linear function
4261
4262        ref_points = []
4263        ref_elevation = []
4264        for i in range(6):
4265            y = (6-i)*25.0
4266            for j in range(5):
4267                x = j*25.0
4268                z = x+2*y
4269
4270                ref_points.append( [x,y] )
4271                ref_elevation.append(z)
4272                fid.write('%f ' %z)
4273            fid.write('\n')
4274
4275        fid.close()
4276
4277        #Write prj file with metadata
4278        metafilename = root+'.prj'
4279        fid = open(metafilename, 'w')
4280
4281
4282        fid.write("""Projection UTM
4283Zone 56
4284Datum WGS84
4285Zunits NO
4286Units METERS
4287Spheroid WGS84
4288Xshift 0.0000000000
4289Yshift 10000000.0000000000
4290Parameters
4291""")
4292        fid.close()
4293
4294        #Convert to NetCDF pts
4295        convert_dem_from_ascii2netcdf(root)
4296        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4297                northing_min=3035.0, northing_max=3125.5)
4298
4299        #Check contents
4300        #Get NetCDF
4301        fid = NetCDFFile(root+'.pts', 'r')
4302
4303        # Get the variables
4304        #print fid.variables.keys()
4305        points = fid.variables['points']
4306        elevation = fid.variables['elevation']
4307
4308        #Check values
4309        assert fid.xllcorner[0] == 2010.0
4310        assert fid.yllcorner[0] == 3035.0
4311
4312        #create new reference points
4313        ref_points = []
4314        ref_elevation = []
4315        for i in range(4):
4316            y = (4-i)*25.0 + 25.0
4317            y_new = y + 3000.5 - 3035.0
4318            for j in range(4):
4319                x = j*25.0 + 25.0
4320                x_new = x + 2000.5 - 2010.0
4321                z = x+2*y
4322
4323                ref_points.append( [x_new,y_new] )
4324                ref_elevation.append(z)
4325
4326        #print points[:]
4327        #print ref_points
4328        assert allclose(points, ref_points)
4329
4330        #print attributes[:]
4331        #print ref_elevation
4332        assert allclose(elevation, ref_elevation)
4333
4334        #Cleanup
4335        fid.close()
4336
4337
4338        os.remove(root + '.pts')
4339        os.remove(root + '.dem')
4340        os.remove(root + '.asc')
4341        os.remove(root + '.prj')
4342
4343
4344
4345    def test_dem2pts_remove_Nullvalues(self):
4346        """Test conversion from dem in ascii format to native NetCDF xya format
4347        """
4348
4349        import time, os
4350        from Numeric import array, zeros, allclose, Float, concatenate
4351        from Scientific.IO.NetCDF import NetCDFFile
4352
4353        #Write test asc file
4354        root = 'demtest'
4355
4356        filename = root+'.asc'
4357        fid = open(filename, 'w')
4358        fid.write("""ncols         5
4359nrows         6
4360xllcorner     2000.5
4361yllcorner     3000.5
4362cellsize      25
4363NODATA_value  -9999
4364""")
4365        #Create linear function
4366        # ref_ will write all the values
4367        # new_ref_ will write the values except for NODATA_values
4368        ref_points = []
4369        ref_elevation = []
4370        new_ref_pts = []
4371        new_ref_elev = []
4372        NODATA_value = -9999
4373        for i in range(6):
4374            y = (6-i)*25.0
4375            for j in range(5):
4376                x = j*25.0
4377                z = x+2*y
4378                if j == 4: z = NODATA_value # column
4379                if i == 2 and j == 2: z = NODATA_value # random
4380                if i == 5 and j == 1: z = NODATA_value
4381                if i == 1: z = NODATA_value # row
4382                if i == 3 and j == 1: z = NODATA_value # two pts/row
4383                if i == 3 and j == 3: z = NODATA_value
4384
4385
4386                if z <> NODATA_value:
4387                    new_ref_elev.append(z)
4388                    new_ref_pts.append( [x,y] )
4389
4390                ref_points.append( [x,y] )
4391                ref_elevation.append(z)
4392
4393                fid.write('%f ' %z)
4394            fid.write('\n')
4395
4396        fid.close()
4397
4398
4399        #Write prj file with metadata
4400        metafilename = root+'.prj'
4401        fid = open(metafilename, 'w')
4402
4403
4404        fid.write("""Projection UTM
4405Zone 56
4406Datum WGS84
4407Zunits NO
4408Units METERS
4409Spheroid WGS84
4410Xshift 0.0000000000
4411Yshift 10000000.0000000000
4412Parameters
4413""")
4414        fid.close()
4415
4416        #Convert to NetCDF pts
4417        convert_dem_from_ascii2netcdf(root)
4418        dem2pts(root)
4419
4420        #Check contents
4421        #Get NetCDF
4422        fid = NetCDFFile(root+'.pts', 'r')
4423
4424        # Get the variables
4425        #print fid.variables.keys()
4426        points = fid.variables['points']
4427        elevation = fid.variables['elevation']
4428
4429        #Check values
4430        #print 'points', points[:]
4431        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
4432        assert allclose(points, new_ref_pts), 'points do not align'
4433
4434        #print 'elevation', elevation[:]
4435        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
4436        assert allclose(elevation, new_ref_elev), 'elevations do not align'
4437
4438        #Cleanup
4439        fid.close()
4440
4441
4442        os.remove(root + '.pts')
4443        os.remove(root + '.dem')
4444        os.remove(root + '.asc')
4445        os.remove(root + '.prj')
4446
4447    def test_dem2pts_bounding_box_Nullvalues(self):
4448        """Test conversion from dem in ascii format to native NetCDF xya format
4449        """
4450
4451        import time, os
4452        from Numeric import array, zeros, allclose, Float, concatenate
4453        from Scientific.IO.NetCDF import NetCDFFile
4454
4455        #Write test asc file
4456        root = 'demtest'
4457
4458        filename = root+'.asc'
4459        fid = open(filename, 'w')
4460        fid.write("""ncols         5
4461nrows         6
4462xllcorner     2000.5
4463yllcorner     3000.5
4464cellsize      25
4465NODATA_value  -9999
4466""")
4467        #Create linear function
4468
4469        ref_points = []
4470        ref_elevation = []
4471        new_ref_pts1 = []
4472        new_ref_elev1 = []
4473        NODATA_value = -9999
4474        for i in range(6):
4475            y = (6-i)*25.0
4476            for j in range(5):
4477                x = j*25.0
4478                z = x+2*y
4479                if j == 4: z = NODATA_value # column
4480                if i == 2 and j == 2: z = NODATA_value # random
4481                if i == 5 and j == 1: z = NODATA_value
4482                if i == 1: z = NODATA_value # row
4483                if i == 3 and j == 1: z = NODATA_value # two pts/row
4484                if i == 3 and j == 3: z = NODATA_value
4485
4486                if z <> NODATA_value:
4487                    new_ref_elev1.append(z)
4488                    new_ref_pts1.append( [x,y] )
4489
4490                ref_points.append( [x,y] )
4491                ref_elevation.append(z)
4492                fid.write('%f ' %z)
4493            fid.write('\n')
4494
4495        fid.close()
4496
4497        #Write prj file with metadata
4498        metafilename = root+'.prj'
4499        fid = open(metafilename, 'w')
4500
4501
4502        fid.write("""Projection UTM
4503Zone 56
4504Datum WGS84
4505Zunits NO
4506Units METERS
4507Spheroid WGS84
4508Xshift 0.0000000000
4509Yshift 10000000.0000000000
4510Parameters
4511""")
4512        fid.close()
4513
4514        #Convert to NetCDF pts
4515        convert_dem_from_ascii2netcdf(root)
4516        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
4517                northing_min=3035.0, northing_max=3125.5)
4518
4519        #Check contents
4520        #Get NetCDF
4521        fid = NetCDFFile(root+'.pts', 'r')
4522
4523        # Get the variables
4524        #print fid.variables.keys()
4525        points = fid.variables['points']
4526        elevation = fid.variables['elevation']
4527
4528        #Check values
4529        assert fid.xllcorner[0] == 2010.0
4530        assert fid.yllcorner[0] == 3035.0
4531
4532        #create new reference points
4533        ref_points = []
4534        ref_elevation = []
4535        new_ref_pts2 = []
4536        new_ref_elev2 = []
4537        for i in range(4):
4538            y = (4-i)*25.0 + 25.0
4539            y_new = y + 3000.5 - 3035.0
4540            for j in range(4):
4541                x = j*25.0 + 25.0
4542                x_new = x + 2000.5 - 2010.0
4543                z = x+2*y
4544
4545                if j == 3: z = NODATA_value # column
4546                if i == 1 and j == 1: z = NODATA_value # random
4547                if i == 4 and j == 0: z = NODATA_value
4548                if i == 0: z = NODATA_value # row
4549                if i == 2 and j == 0: z = NODATA_value # two pts/row
4550                if i == 2 and j == 2: z = NODATA_value
4551
4552                if z <> NODATA_value:
4553                    new_ref_elev2.append(z)
4554                    new_ref_pts2.append( [x_new,y_new] )
4555
4556
4557                ref_points.append( [x_new,y_new] )
4558                ref_elevation.append(z)
4559
4560        #print points[:]
4561        #print ref_points
4562        #assert allclose(points, ref_points)
4563
4564        #print attributes[:]
4565        #print ref_elevation
4566        #assert allclose(elevation, ref_elevation)
4567
4568
4569        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
4570        assert allclose(points, new_ref_pts2), 'points do not align'
4571
4572        #print 'elevation', elevation[:]
4573        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
4574        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
4575        #Cleanup
4576        fid.close()
4577
4578
4579        os.remove(root + '.pts')
4580        os.remove(root + '.dem')
4581        os.remove(root + '.asc')
4582        os.remove(root + '.prj')
4583
4584
4585
Note: See TracBrowser for help on using the repository browser.