source: inundation/pyvolution/test_data_manager.py @ 3086

Last change on this file since 3086 was 2891, checked in by ole, 19 years ago

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