source: inundation-numpy-branch/pyvolution/test_data_manager.py @ 3514

Last change on this file since 3514 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 121.5 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import copy
6from Numeric import zeros, array, allclose, Float
7from anuga.utilities.numerical_tools import mean
8import tempfile
9import os
10from Scientific.IO.NetCDF import NetCDFFile
11
12from anuga.pyvolution.data_manager import *
13from shallow_water import *
14from config import epsilon
15import anuga.pyvolution.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
574
575    def test_dem2pts(self):
576        """Test conversion from dem in ascii format to native NetCDF xya format
577        """
578
579        import time, os
580        from Numeric import array, zeros, allclose, Float, concatenate
581        from Scientific.IO.NetCDF import NetCDFFile
582
583        #Write test asc file
584        root = 'demtest'
585
586        filename = root+'.asc'
587        fid = open(filename, 'w')
588        fid.write("""ncols         5
589nrows         6
590xllcorner     2000.5
591yllcorner     3000.5
592cellsize      25
593NODATA_value  -9999
594""")
595        #Create linear function
596
597        ref_points = []
598        ref_elevation = []
599        for i in range(6):
600            y = (6-i)*25.0
601            for j in range(5):
602                x = j*25.0
603                z = x+2*y
604
605                ref_points.append( [x,y] )
606                ref_elevation.append(z)
607                fid.write('%f ' %z)
608            fid.write('\n')
609
610        fid.close()
611
612        #Write prj file with metadata
613        metafilename = root+'.prj'
614        fid = open(metafilename, 'w')
615
616
617        fid.write("""Projection UTM
618Zone 56
619Datum WGS84
620Zunits NO
621Units METERS
622Spheroid WGS84
623Xshift 0.0000000000
624Yshift 10000000.0000000000
625Parameters
626""")
627        fid.close()
628
629        #Convert to NetCDF pts
630        convert_dem_from_ascii2netcdf(root)
631        dem2pts(root)
632
633        #Check contents
634        #Get NetCDF
635        fid = NetCDFFile(root+'.pts', 'r')
636
637        # Get the variables
638        #print fid.variables.keys()
639        points = fid.variables['points']
640        elevation = fid.variables['elevation']
641
642        #Check values
643
644        #print points[:]
645        #print ref_points
646        assert allclose(points, ref_points)
647
648        #print attributes[:]
649        #print ref_elevation
650        assert allclose(elevation, ref_elevation)
651
652        #Cleanup
653        fid.close()
654
655
656        os.remove(root + '.pts')
657        os.remove(root + '.dem')
658        os.remove(root + '.asc')
659        os.remove(root + '.prj')
660
661
662
663    def test_dem2pts_bounding_box(self):
664        """Test conversion from dem in ascii format to native NetCDF xya format
665        """
666
667        import time, os
668        from Numeric import array, zeros, allclose, Float, concatenate
669        from Scientific.IO.NetCDF import NetCDFFile
670
671        #Write test asc file
672        root = 'demtest'
673
674        filename = root+'.asc'
675        fid = open(filename, 'w')
676        fid.write("""ncols         5
677nrows         6
678xllcorner     2000.5
679yllcorner     3000.5
680cellsize      25
681NODATA_value  -9999
682""")
683        #Create linear function
684
685        ref_points = []
686        ref_elevation = []
687        for i in range(6):
688            y = (6-i)*25.0
689            for j in range(5):
690                x = j*25.0
691                z = x+2*y
692
693                ref_points.append( [x,y] )
694                ref_elevation.append(z)
695                fid.write('%f ' %z)
696            fid.write('\n')
697
698        fid.close()
699
700        #Write prj file with metadata
701        metafilename = root+'.prj'
702        fid = open(metafilename, 'w')
703
704
705        fid.write("""Projection UTM
706Zone 56
707Datum WGS84
708Zunits NO
709Units METERS
710Spheroid WGS84
711Xshift 0.0000000000
712Yshift 10000000.0000000000
713Parameters
714""")
715        fid.close()
716
717        #Convert to NetCDF pts
718        convert_dem_from_ascii2netcdf(root)
719        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
720                northing_min=3035.0, northing_max=3125.5)
721
722        #Check contents
723        #Get NetCDF
724        fid = NetCDFFile(root+'.pts', 'r')
725
726        # Get the variables
727        #print fid.variables.keys()
728        points = fid.variables['points']
729        elevation = fid.variables['elevation']
730
731        #Check values
732        assert fid.xllcorner[0] == 2010.0
733        assert fid.yllcorner[0] == 3035.0
734
735        #create new reference points
736        ref_points = []
737        ref_elevation = []
738        for i in range(4):
739            y = (4-i)*25.0 + 25.0
740            y_new = y + 3000.5 - 3035.0
741            for j in range(4):
742                x = j*25.0 + 25.0
743                x_new = x + 2000.5 - 2010.0
744                z = x+2*y
745
746                ref_points.append( [x_new,y_new] )
747                ref_elevation.append(z)
748
749        #print points[:]
750        #print ref_points
751        assert allclose(points, ref_points)
752
753        #print attributes[:]
754        #print ref_elevation
755        assert allclose(elevation, ref_elevation)
756
757        #Cleanup
758        fid.close()
759
760
761        os.remove(root + '.pts')
762        os.remove(root + '.dem')
763        os.remove(root + '.asc')
764        os.remove(root + '.prj')
765
766    def test_dem2pts_remove_Nullvalues(self):
767        """Test conversion from dem in ascii format to native NetCDF xya format
768        """
769
770        import time, os
771        from Numeric import array, zeros, allclose, Float, concatenate
772        from Scientific.IO.NetCDF import NetCDFFile
773
774        #Write test asc file
775        root = 'demtest'
776
777        filename = root+'.asc'
778        fid = open(filename, 'w')
779        fid.write("""ncols         5
780nrows         6
781xllcorner     2000.5
782yllcorner     3000.5
783cellsize      25
784NODATA_value  -9999
785""")
786        #Create linear function
787        # ref_ will write all the values
788        # new_ref_ will write the values except for NODATA_values
789        ref_points = []
790        ref_elevation = []
791        new_ref_pts = []
792        new_ref_elev = []
793        NODATA_value = -9999
794        for i in range(6):
795            y = (6-i)*25.0
796            for j in range(5):
797                x = j*25.0
798                z = x+2*y
799                if j == 4: z = NODATA_value # column
800                if i == 2 and j == 2: z = NODATA_value # random
801                if i == 5 and j == 1: z = NODATA_value
802                if i == 1: z = NODATA_value # row
803                if i == 3 and j == 1: z = NODATA_value # two pts/row
804                if i == 3 and j == 3: z = NODATA_value
805
806               
807                if z <> NODATA_value:
808                    new_ref_elev.append(z)
809                    new_ref_pts.append( [x,y] )
810               
811                ref_points.append( [x,y] )
812                ref_elevation.append(z)
813
814                fid.write('%f ' %z)
815            fid.write('\n')
816
817        fid.close()
818
819       
820        #Write prj file with metadata
821        metafilename = root+'.prj'
822        fid = open(metafilename, 'w')
823
824
825        fid.write("""Projection UTM
826Zone 56
827Datum WGS84
828Zunits NO
829Units METERS
830Spheroid WGS84
831Xshift 0.0000000000
832Yshift 10000000.0000000000
833Parameters
834""")
835        fid.close()
836
837        #Convert to NetCDF pts
838        convert_dem_from_ascii2netcdf(root)
839        dem2pts(root)
840
841        #Check contents
842        #Get NetCDF
843        fid = NetCDFFile(root+'.pts', 'r')
844
845        # Get the variables
846        #print fid.variables.keys()
847        points = fid.variables['points']
848        elevation = fid.variables['elevation']
849
850        #Check values
851        #print 'points', points[:]
852        assert len(points) == len(new_ref_pts), 'length of returned points not correct'
853        assert allclose(points, new_ref_pts), 'points do not align'
854
855        #print 'elevation', elevation[:]
856        assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct'
857        assert allclose(elevation, new_ref_elev), 'elevations do not align'
858
859        #Cleanup
860        fid.close()
861
862
863        os.remove(root + '.pts')
864        os.remove(root + '.dem')
865        os.remove(root + '.asc')
866        os.remove(root + '.prj')
867
868    def test_dem2pts_bounding_box_Nullvalues(self):
869        """Test conversion from dem in ascii format to native NetCDF xya format
870        """
871
872        import time, os
873        from Numeric import array, zeros, allclose, Float, concatenate
874        from Scientific.IO.NetCDF import NetCDFFile
875
876        #Write test asc file
877        root = 'demtest'
878
879        filename = root+'.asc'
880        fid = open(filename, 'w')
881        fid.write("""ncols         5
882nrows         6
883xllcorner     2000.5
884yllcorner     3000.5
885cellsize      25
886NODATA_value  -9999
887""")
888        #Create linear function
889
890        ref_points = []
891        ref_elevation = []
892        new_ref_pts1 = []
893        new_ref_elev1 = []
894        NODATA_value = -9999
895        for i in range(6):
896            y = (6-i)*25.0
897            for j in range(5):
898                x = j*25.0
899                z = x+2*y
900                if j == 4: z = NODATA_value # column
901                if i == 2 and j == 2: z = NODATA_value # random
902                if i == 5 and j == 1: z = NODATA_value
903                if i == 1: z = NODATA_value # row
904                if i == 3 and j == 1: z = NODATA_value # two pts/row
905                if i == 3 and j == 3: z = NODATA_value
906
907                if z <> NODATA_value:
908                    new_ref_elev1.append(z)
909                    new_ref_pts1.append( [x,y] )
910                   
911                ref_points.append( [x,y] )
912                ref_elevation.append(z)
913                fid.write('%f ' %z)
914            fid.write('\n')
915
916        fid.close()
917
918        #Write prj file with metadata
919        metafilename = root+'.prj'
920        fid = open(metafilename, 'w')
921
922
923        fid.write("""Projection UTM
924Zone 56
925Datum WGS84
926Zunits NO
927Units METERS
928Spheroid WGS84
929Xshift 0.0000000000
930Yshift 10000000.0000000000
931Parameters
932""")
933        fid.close()
934       
935        #Convert to NetCDF pts
936        convert_dem_from_ascii2netcdf(root)
937        dem2pts(root, easting_min=2010.0, easting_max=2110.0,
938                northing_min=3035.0, northing_max=3125.5)
939
940        #Check contents
941        #Get NetCDF
942        fid = NetCDFFile(root+'.pts', 'r')
943
944        # Get the variables
945        #print fid.variables.keys()
946        points = fid.variables['points']
947        elevation = fid.variables['elevation']
948
949        #Check values
950        assert fid.xllcorner[0] == 2010.0
951        assert fid.yllcorner[0] == 3035.0
952
953        #create new reference points
954        ref_points = []
955        ref_elevation = []
956        new_ref_pts2 = []
957        new_ref_elev2 = []
958        for i in range(4):
959            y = (4-i)*25.0 + 25.0
960            y_new = y + 3000.5 - 3035.0
961            for j in range(4):
962                x = j*25.0 + 25.0
963                x_new = x + 2000.5 - 2010.0
964                z = x+2*y
965
966                if j == 3: z = NODATA_value # column
967                if i == 1 and j == 1: z = NODATA_value # random
968                if i == 4 and j == 0: z = NODATA_value
969                if i == 0: z = NODATA_value # row
970                if i == 2 and j == 0: z = NODATA_value # two pts/row
971                if i == 2 and j == 2: z = NODATA_value
972
973                if z <> NODATA_value:
974                    new_ref_elev2.append(z)
975                    new_ref_pts2.append( [x_new,y_new] )
976
977                   
978                ref_points.append( [x_new,y_new] )
979                ref_elevation.append(z)
980
981        #print points[:]
982        #print ref_points
983        #assert allclose(points, ref_points)
984
985        #print attributes[:]
986        #print ref_elevation
987        #assert allclose(elevation, ref_elevation)
988
989       
990        assert len(points) == len(new_ref_pts2), 'length of returned points not correct'
991        assert allclose(points, new_ref_pts2), 'points do not align'
992
993        #print 'elevation', elevation[:]
994        assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct'
995        assert allclose(elevation, new_ref_elev2), 'elevations do not align'
996        #Cleanup
997        fid.close()
998
999
1000        os.remove(root + '.pts')
1001        os.remove(root + '.dem')
1002        os.remove(root + '.asc')
1003        os.remove(root + '.prj')
1004
1005
1006
1007
1008    def test_hecras_cross_sections2pts(self):
1009        """Test conversion from HECRAS cross sections in ascii format
1010        to native NetCDF pts format
1011        """
1012
1013        import time, os
1014        from Numeric import array, zeros, allclose, Float, concatenate
1015        from Scientific.IO.NetCDF import NetCDFFile
1016
1017        #Write test asc file
1018        root = 'hecrastest'
1019
1020        filename = root+'.sdf'
1021        fid = open(filename, 'w')
1022        fid.write("""
1023# RAS export file created on Mon 15Aug2005 11:42
1024# by HEC-RAS Version 3.1.1
1025
1026BEGIN HEADER:
1027  UNITS: METRIC
1028  DTM TYPE: TIN
1029  DTM: v:\1\cit\perth_topo\river_tin
1030  STREAM LAYER: c:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
1031  CROSS-SECTION LAYER: c:\\x_local\hecras\21_02_03\up_can_xs3d.shp
1032  MAP PROJECTION: UTM
1033  PROJECTION ZONE: 50
1034  DATUM: AGD66
1035  VERTICAL DATUM:
1036  NUMBER OF REACHES:  19
1037  NUMBER OF CROSS-SECTIONS:  2
1038END HEADER:
1039
1040
1041BEGIN CROSS-SECTIONS:
1042
1043  CROSS-SECTION:
1044    STREAM ID:Southern-Wungong
1045    REACH ID:Southern-Wungong
1046    STATION:21410   
1047    CUT LINE:
1048      407546.08 , 6437277.542
1049      407329.32 , 6437489.482
1050      407283.11 , 6437541.232
1051    SURFACE LINE:
1052     407546.08,   6437277.54,   52.14
1053     407538.88,   6437284.58,   51.07
1054     407531.68,   6437291.62,   50.56
1055     407524.48,   6437298.66,   49.58
1056     407517.28,   6437305.70,   49.09
1057     407510.08,   6437312.74,   48.76
1058  END:
1059
1060  CROSS-SECTION:
1061    STREAM ID:Swan River     
1062    REACH ID:Swan Mouth     
1063    STATION:840.*   
1064    CUT LINE:
1065      381178.0855 , 6452559.0685
1066      380485.4755 , 6453169.272
1067    SURFACE LINE:
1068     381178.09,   6452559.07,   4.17
1069     381169.49,   6452566.64,   4.26
1070     381157.78,   6452576.96,   4.34
1071     381155.97,   6452578.56,   4.35
1072     381143.72,   6452589.35,   4.43
1073     381136.69,   6452595.54,   4.58
1074     381114.74,   6452614.88,   4.41
1075     381075.53,   6452649.43,   4.17
1076     381071.47,   6452653.00,   3.99
1077     381063.46,   6452660.06,   3.67
1078     381054.41,   6452668.03,   3.67
1079  END: 
1080END CROSS-SECTIONS:
1081""")
1082
1083        fid.close()
1084
1085
1086        #Convert to NetCDF pts
1087        hecras_cross_sections2pts(root)
1088
1089        #Check contents
1090        #Get NetCDF
1091        fid = NetCDFFile(root+'.pts', 'r')
1092
1093        # Get the variables
1094        #print fid.variables.keys()
1095        points = fid.variables['points']
1096        elevation = fid.variables['elevation']
1097
1098        #Check values
1099        ref_points = [[407546.08, 6437277.54],
1100                      [407538.88, 6437284.58],
1101                      [407531.68, 6437291.62],
1102                      [407524.48, 6437298.66],
1103                      [407517.28, 6437305.70],
1104                      [407510.08, 6437312.74]]
1105
1106        ref_points += [[381178.09, 6452559.07], 
1107                       [381169.49, 6452566.64], 
1108                       [381157.78, 6452576.96], 
1109                       [381155.97, 6452578.56], 
1110                       [381143.72, 6452589.35], 
1111                       [381136.69, 6452595.54], 
1112                       [381114.74, 6452614.88], 
1113                       [381075.53, 6452649.43], 
1114                       [381071.47, 6452653.00], 
1115                       [381063.46, 6452660.06], 
1116                       [381054.41, 6452668.03]] 
1117       
1118                     
1119        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
1120        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
1121
1122        #print points[:]
1123        #print ref_points
1124        assert allclose(points, ref_points)
1125
1126        #print attributes[:]
1127        #print ref_elevation
1128        assert allclose(elevation, ref_elevation)
1129
1130        #Cleanup
1131        fid.close()
1132
1133
1134        os.remove(root + '.sdf')
1135        os.remove(root + '.pts')
1136
1137
1138
1139    def test_sww2dem_asc_elevation(self):
1140        """Test that sww information can be converted correctly to asc/prj
1141        format readable by e.g. ArcView
1142        """
1143
1144        import time, os
1145        from Numeric import array, zeros, allclose, Float, concatenate
1146        from Scientific.IO.NetCDF import NetCDFFile
1147
1148        #Setup
1149        self.domain.filename = 'datatest'
1150
1151        prjfile = self.domain.filename + '_elevation.prj'
1152        ascfile = self.domain.filename + '_elevation.asc'
1153        swwfile = self.domain.filename + '.sww'
1154
1155        self.domain.set_datadir('.')
1156        self.domain.format = 'sww'
1157        self.domain.smooth = True
1158        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1159
1160        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1161
1162        sww = get_dataobject(self.domain)
1163        sww.store_connectivity()
1164        sww.store_timestep('stage')
1165
1166        self.domain.evolve_to_end(finaltime = 0.01)
1167        sww.store_timestep('stage')
1168
1169        cellsize = 0.25
1170        #Check contents
1171        #Get NetCDF
1172
1173        fid = NetCDFFile(sww.filename, 'r')
1174
1175        # Get the variables
1176        x = fid.variables['x'][:]
1177        y = fid.variables['y'][:]
1178        z = fid.variables['elevation'][:]
1179        time = fid.variables['time'][:]
1180        stage = fid.variables['stage'][:]
1181
1182
1183        #Export to ascii/prj files
1184        sww2dem(self.domain.filename,
1185                quantity = 'elevation',
1186                cellsize = cellsize,
1187                verbose = False,
1188                format = 'asc')
1189               
1190        #Check prj (meta data)
1191        prjid = open(prjfile)
1192        lines = prjid.readlines()
1193        prjid.close()
1194
1195        L = lines[0].strip().split()
1196        assert L[0].strip().lower() == 'projection'
1197        assert L[1].strip().lower() == 'utm'
1198
1199        L = lines[1].strip().split()
1200        assert L[0].strip().lower() == 'zone'
1201        assert L[1].strip().lower() == '56'
1202
1203        L = lines[2].strip().split()
1204        assert L[0].strip().lower() == 'datum'
1205        assert L[1].strip().lower() == 'wgs84'
1206
1207        L = lines[3].strip().split()
1208        assert L[0].strip().lower() == 'zunits'
1209        assert L[1].strip().lower() == 'no'
1210
1211        L = lines[4].strip().split()
1212        assert L[0].strip().lower() == 'units'
1213        assert L[1].strip().lower() == 'meters'
1214
1215        L = lines[5].strip().split()
1216        assert L[0].strip().lower() == 'spheroid'
1217        assert L[1].strip().lower() == 'wgs84'
1218
1219        L = lines[6].strip().split()
1220        assert L[0].strip().lower() == 'xshift'
1221        assert L[1].strip().lower() == '500000'
1222
1223        L = lines[7].strip().split()
1224        assert L[0].strip().lower() == 'yshift'
1225        assert L[1].strip().lower() == '10000000'
1226
1227        L = lines[8].strip().split()
1228        assert L[0].strip().lower() == 'parameters'
1229
1230
1231        #Check asc file
1232        ascid = open(ascfile)
1233        lines = ascid.readlines()
1234        ascid.close()
1235
1236        L = lines[0].strip().split()
1237        assert L[0].strip().lower() == 'ncols'
1238        assert L[1].strip().lower() == '5'
1239
1240        L = lines[1].strip().split()
1241        assert L[0].strip().lower() == 'nrows'
1242        assert L[1].strip().lower() == '5'
1243
1244        L = lines[2].strip().split()
1245        assert L[0].strip().lower() == 'xllcorner'
1246        assert allclose(float(L[1].strip().lower()), 308500)
1247
1248        L = lines[3].strip().split()
1249        assert L[0].strip().lower() == 'yllcorner'
1250        assert allclose(float(L[1].strip().lower()), 6189000)
1251
1252        L = lines[4].strip().split()
1253        assert L[0].strip().lower() == 'cellsize'
1254        assert allclose(float(L[1].strip().lower()), cellsize)
1255
1256        L = lines[5].strip().split()
1257        assert L[0].strip() == 'NODATA_value'
1258        assert L[1].strip().lower() == '-9999'
1259
1260        #Check grid values
1261        for j in range(5):
1262            L = lines[6+j].strip().split()
1263            y = (4-j) * cellsize
1264            for i in range(5):
1265                assert allclose(float(L[i]), -i*cellsize - y)
1266
1267
1268        fid.close()
1269
1270        #Cleanup
1271        os.remove(prjfile)
1272        os.remove(ascfile)
1273        os.remove(swwfile)
1274
1275
1276
1277    def test_sww2dem_larger(self):
1278        """Test that sww information can be converted correctly to asc/prj
1279        format readable by e.g. ArcView. Here:
1280
1281        ncols         11
1282        nrows         11
1283        xllcorner     308500
1284        yllcorner     6189000
1285        cellsize      10.000000
1286        NODATA_value  -9999
1287        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
1288         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
1289         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
1290         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
1291         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
1292         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
1293         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
1294         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
1295         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
1296         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
1297           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
1298       
1299        """
1300
1301        import time, os
1302        from Numeric import array, zeros, allclose, Float, concatenate
1303        from Scientific.IO.NetCDF import NetCDFFile
1304
1305        #Setup
1306
1307        from mesh_factory import rectangular
1308
1309        #Create basic mesh (100m x 100m)
1310        points, vertices, boundary = rectangular(2, 2, 100, 100)
1311
1312        #Create shallow water domain
1313        domain = Domain(points, vertices, boundary)
1314        domain.default_order = 2
1315
1316        domain.filename = 'datatest'
1317
1318        prjfile = domain.filename + '_elevation.prj'
1319        ascfile = domain.filename + '_elevation.asc'
1320        swwfile = domain.filename + '.sww'
1321
1322        domain.set_datadir('.')
1323        domain.format = 'sww'
1324        domain.smooth = True
1325        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1326
1327        #
1328        domain.set_quantity('elevation', lambda x,y: -x-y)
1329        domain.set_quantity('stage', 0)       
1330
1331        B = Transmissive_boundary(domain)
1332        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1333
1334
1335        #
1336        sww = get_dataobject(domain)
1337        sww.store_connectivity()
1338        sww.store_timestep('stage')
1339
1340        domain.evolve_to_end(finaltime = 0.01)
1341        sww.store_timestep('stage')
1342
1343        cellsize = 10  #10m grid
1344
1345       
1346        #Check contents
1347        #Get NetCDF
1348
1349        fid = NetCDFFile(sww.filename, 'r')
1350
1351        # Get the variables
1352        x = fid.variables['x'][:]
1353        y = fid.variables['y'][:]
1354        z = fid.variables['elevation'][:]
1355        time = fid.variables['time'][:]
1356        stage = fid.variables['stage'][:]
1357
1358
1359        #Export to ascii/prj files
1360        sww2dem(domain.filename,
1361                quantity = 'elevation',
1362                cellsize = cellsize,
1363                verbose = False,
1364                format = 'asc')
1365       
1366               
1367        #Check prj (meta data)
1368        prjid = open(prjfile)
1369        lines = prjid.readlines()
1370        prjid.close()
1371
1372        L = lines[0].strip().split()
1373        assert L[0].strip().lower() == 'projection'
1374        assert L[1].strip().lower() == 'utm'
1375
1376        L = lines[1].strip().split()
1377        assert L[0].strip().lower() == 'zone'
1378        assert L[1].strip().lower() == '56'
1379
1380        L = lines[2].strip().split()
1381        assert L[0].strip().lower() == 'datum'
1382        assert L[1].strip().lower() == 'wgs84'
1383
1384        L = lines[3].strip().split()
1385        assert L[0].strip().lower() == 'zunits'
1386        assert L[1].strip().lower() == 'no'
1387
1388        L = lines[4].strip().split()
1389        assert L[0].strip().lower() == 'units'
1390        assert L[1].strip().lower() == 'meters'
1391
1392        L = lines[5].strip().split()
1393        assert L[0].strip().lower() == 'spheroid'
1394        assert L[1].strip().lower() == 'wgs84'
1395
1396        L = lines[6].strip().split()
1397        assert L[0].strip().lower() == 'xshift'
1398        assert L[1].strip().lower() == '500000'
1399
1400        L = lines[7].strip().split()
1401        assert L[0].strip().lower() == 'yshift'
1402        assert L[1].strip().lower() == '10000000'
1403
1404        L = lines[8].strip().split()
1405        assert L[0].strip().lower() == 'parameters'
1406
1407
1408        #Check asc file
1409        ascid = open(ascfile)
1410        lines = ascid.readlines()
1411        ascid.close()
1412
1413        L = lines[0].strip().split()
1414        assert L[0].strip().lower() == 'ncols'
1415        assert L[1].strip().lower() == '11'
1416
1417        L = lines[1].strip().split()
1418        assert L[0].strip().lower() == 'nrows'
1419        assert L[1].strip().lower() == '11'
1420
1421        L = lines[2].strip().split()
1422        assert L[0].strip().lower() == 'xllcorner'
1423        assert allclose(float(L[1].strip().lower()), 308500)
1424
1425        L = lines[3].strip().split()
1426        assert L[0].strip().lower() == 'yllcorner'
1427        assert allclose(float(L[1].strip().lower()), 6189000)
1428
1429        L = lines[4].strip().split()
1430        assert L[0].strip().lower() == 'cellsize'
1431        assert allclose(float(L[1].strip().lower()), cellsize)
1432
1433        L = lines[5].strip().split()
1434        assert L[0].strip() == 'NODATA_value'
1435        assert L[1].strip().lower() == '-9999'
1436
1437        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
1438        for i, line in enumerate(lines[6:]):
1439            for j, value in enumerate( line.split() ):
1440                #assert allclose(float(value), -(10-i+j)*cellsize)
1441                assert float(value) == -(10-i+j)*cellsize
1442       
1443
1444        fid.close()
1445
1446        #Cleanup
1447        os.remove(prjfile)
1448        os.remove(ascfile)
1449        os.remove(swwfile)
1450
1451
1452
1453    def test_sww2dem_boundingbox(self):
1454        """Test that sww information can be converted correctly to asc/prj
1455        format readable by e.g. ArcView.
1456        This will test that mesh can be restricted by bounding box
1457
1458        Original extent is 100m x 100m:
1459
1460        Eastings:   308500 -  308600
1461        Northings: 6189000 - 6189100
1462
1463        Bounding box changes this to the 50m x 50m square defined by
1464
1465        Eastings:   308530 -  308570
1466        Northings: 6189050 - 6189100
1467
1468        The cropped values should be
1469
1470         -130 -140 -150 -160 -170 
1471         -120 -130 -140 -150 -160 
1472         -110 -120 -130 -140 -150 
1473         -100 -110 -120 -130 -140 
1474          -90 -100 -110 -120 -130 
1475          -80  -90 -100 -110 -120
1476
1477        and the new lower reference point should be   
1478        Eastings:   308530
1479        Northings: 6189050
1480       
1481        Original dataset is the same as in test_sww2dem_larger()
1482       
1483        """
1484
1485        import time, os
1486        from Numeric import array, zeros, allclose, Float, concatenate
1487        from Scientific.IO.NetCDF import NetCDFFile
1488
1489        #Setup
1490
1491        from mesh_factory import rectangular
1492
1493        #Create basic mesh (100m x 100m)
1494        points, vertices, boundary = rectangular(2, 2, 100, 100)
1495
1496        #Create shallow water domain
1497        domain = Domain(points, vertices, boundary)
1498        domain.default_order = 2
1499
1500        domain.filename = 'datatest'
1501
1502        prjfile = domain.filename + '_elevation.prj'
1503        ascfile = domain.filename + '_elevation.asc'
1504        swwfile = domain.filename + '.sww'
1505
1506        domain.set_datadir('.')
1507        domain.format = 'sww'
1508        domain.smooth = True
1509        domain.geo_reference = Geo_reference(56, 308500, 6189000)
1510
1511        #
1512        domain.set_quantity('elevation', lambda x,y: -x-y)
1513        domain.set_quantity('stage', 0)       
1514
1515        B = Transmissive_boundary(domain)
1516        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
1517
1518
1519        #
1520        sww = get_dataobject(domain)
1521        sww.store_connectivity()
1522        sww.store_timestep('stage')
1523
1524        domain.evolve_to_end(finaltime = 0.01)
1525        sww.store_timestep('stage')
1526
1527        cellsize = 10  #10m grid
1528
1529       
1530        #Check contents
1531        #Get NetCDF
1532
1533        fid = NetCDFFile(sww.filename, 'r')
1534
1535        # Get the variables
1536        x = fid.variables['x'][:]
1537        y = fid.variables['y'][:]
1538        z = fid.variables['elevation'][:]
1539        time = fid.variables['time'][:]
1540        stage = fid.variables['stage'][:]
1541
1542
1543        #Export to ascii/prj files
1544        sww2dem(domain.filename,
1545                quantity = 'elevation',
1546                cellsize = cellsize,
1547                easting_min = 308530,
1548                easting_max = 308570,
1549                northing_min = 6189050, 
1550                northing_max = 6189100,                               
1551                verbose = False,
1552                format = 'asc')
1553       
1554        fid.close()
1555
1556       
1557        #Check prj (meta data)
1558        prjid = open(prjfile)
1559        lines = prjid.readlines()
1560        prjid.close()
1561
1562        L = lines[0].strip().split()
1563        assert L[0].strip().lower() == 'projection'
1564        assert L[1].strip().lower() == 'utm'
1565
1566        L = lines[1].strip().split()
1567        assert L[0].strip().lower() == 'zone'
1568        assert L[1].strip().lower() == '56'
1569
1570        L = lines[2].strip().split()
1571        assert L[0].strip().lower() == 'datum'
1572        assert L[1].strip().lower() == 'wgs84'
1573
1574        L = lines[3].strip().split()
1575        assert L[0].strip().lower() == 'zunits'
1576        assert L[1].strip().lower() == 'no'
1577
1578        L = lines[4].strip().split()
1579        assert L[0].strip().lower() == 'units'
1580        assert L[1].strip().lower() == 'meters'
1581
1582        L = lines[5].strip().split()
1583        assert L[0].strip().lower() == 'spheroid'
1584        assert L[1].strip().lower() == 'wgs84'
1585
1586        L = lines[6].strip().split()
1587        assert L[0].strip().lower() == 'xshift'
1588        assert L[1].strip().lower() == '500000'
1589
1590        L = lines[7].strip().split()
1591        assert L[0].strip().lower() == 'yshift'
1592        assert L[1].strip().lower() == '10000000'
1593
1594        L = lines[8].strip().split()
1595        assert L[0].strip().lower() == 'parameters'
1596
1597
1598        #Check asc file
1599        ascid = open(ascfile)
1600        lines = ascid.readlines()
1601        ascid.close()
1602
1603        L = lines[0].strip().split()
1604        assert L[0].strip().lower() == 'ncols'
1605        assert L[1].strip().lower() == '5'
1606
1607        L = lines[1].strip().split()
1608        assert L[0].strip().lower() == 'nrows'
1609        assert L[1].strip().lower() == '6'
1610
1611        L = lines[2].strip().split()
1612        assert L[0].strip().lower() == 'xllcorner'
1613        assert allclose(float(L[1].strip().lower()), 308530)
1614
1615        L = lines[3].strip().split()
1616        assert L[0].strip().lower() == 'yllcorner'
1617        assert allclose(float(L[1].strip().lower()), 6189050)
1618
1619        L = lines[4].strip().split()
1620        assert L[0].strip().lower() == 'cellsize'
1621        assert allclose(float(L[1].strip().lower()), cellsize)
1622
1623        L = lines[5].strip().split()
1624        assert L[0].strip() == 'NODATA_value'
1625        assert L[1].strip().lower() == '-9999'
1626
1627        #Check grid values
1628        for i, line in enumerate(lines[6:]):
1629            for j, value in enumerate( line.split() ):
1630                #assert float(value) == -(10-i+j)*cellsize               
1631                assert float(value) == -(10-i+j+3)*cellsize
1632       
1633
1634
1635        #Cleanup
1636        os.remove(prjfile)
1637        os.remove(ascfile)
1638        os.remove(swwfile)
1639
1640
1641
1642    def test_sww2dem_asc_stage_reduction(self):
1643        """Test that sww information can be converted correctly to asc/prj
1644        format readable by e.g. ArcView
1645
1646        This tests the reduction of quantity stage using min
1647        """
1648
1649        import time, os
1650        from Numeric import array, zeros, allclose, Float, concatenate
1651        from Scientific.IO.NetCDF import NetCDFFile
1652
1653        #Setup
1654        self.domain.filename = 'datatest'
1655
1656        prjfile = self.domain.filename + '_stage.prj'
1657        ascfile = self.domain.filename + '_stage.asc'
1658        swwfile = self.domain.filename + '.sww'
1659
1660        self.domain.set_datadir('.')
1661        self.domain.format = 'sww'
1662        self.domain.smooth = True
1663        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1664
1665        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1666
1667
1668        sww = get_dataobject(self.domain)
1669        sww.store_connectivity()
1670        sww.store_timestep('stage')
1671
1672        self.domain.evolve_to_end(finaltime = 0.01)
1673        sww.store_timestep('stage')
1674
1675        cellsize = 0.25
1676        #Check contents
1677        #Get NetCDF
1678
1679        fid = NetCDFFile(sww.filename, 'r')
1680
1681        # Get the variables
1682        x = fid.variables['x'][:]
1683        y = fid.variables['y'][:]
1684        z = fid.variables['elevation'][:]
1685        time = fid.variables['time'][:]
1686        stage = fid.variables['stage'][:]
1687
1688
1689        #Export to ascii/prj files
1690        sww2dem(self.domain.filename,
1691                quantity = 'stage',
1692                cellsize = cellsize,
1693                reduction = min,
1694                format = 'asc')
1695
1696
1697        #Check asc file
1698        ascid = open(ascfile)
1699        lines = ascid.readlines()
1700        ascid.close()
1701
1702        L = lines[0].strip().split()
1703        assert L[0].strip().lower() == 'ncols'
1704        assert L[1].strip().lower() == '5'
1705
1706        L = lines[1].strip().split()
1707        assert L[0].strip().lower() == 'nrows'
1708        assert L[1].strip().lower() == '5'
1709
1710        L = lines[2].strip().split()
1711        assert L[0].strip().lower() == 'xllcorner'
1712        assert allclose(float(L[1].strip().lower()), 308500)
1713
1714        L = lines[3].strip().split()
1715        assert L[0].strip().lower() == 'yllcorner'
1716        assert allclose(float(L[1].strip().lower()), 6189000)
1717
1718        L = lines[4].strip().split()
1719        assert L[0].strip().lower() == 'cellsize'
1720        assert allclose(float(L[1].strip().lower()), cellsize)
1721
1722        L = lines[5].strip().split()
1723        assert L[0].strip() == 'NODATA_value'
1724        assert L[1].strip().lower() == '-9999'
1725
1726
1727        #Check grid values (where applicable)
1728        for j in range(5):
1729            if j%2 == 0:
1730                L = lines[6+j].strip().split()
1731                jj = 4-j
1732                for i in range(5):
1733                    if i%2 == 0:
1734                        index = jj/2 + i/2*3
1735                        val0 = stage[0,index]
1736                        val1 = stage[1,index]
1737
1738                        #print i, j, index, ':', L[i], val0, val1
1739                        assert allclose(float(L[i]), min(val0, val1))
1740
1741
1742        fid.close()
1743
1744        #Cleanup
1745        os.remove(prjfile)
1746        os.remove(ascfile)
1747        #os.remove(swwfile)
1748
1749
1750
1751    def test_sww2dem_asc_derived_quantity(self):
1752        """Test that sww information can be converted correctly to asc/prj
1753        format readable by e.g. ArcView
1754
1755        This tests the use of derived quantities
1756        """
1757
1758        import time, os
1759        from Numeric import array, zeros, allclose, Float, concatenate
1760        from Scientific.IO.NetCDF import NetCDFFile
1761
1762        #Setup
1763        self.domain.filename = 'datatest'
1764
1765        prjfile = self.domain.filename + '_depth.prj'
1766        ascfile = self.domain.filename + '_depth.asc'
1767        swwfile = self.domain.filename + '.sww'
1768
1769        self.domain.set_datadir('.')
1770        self.domain.format = 'sww'
1771        self.domain.smooth = True
1772        self.domain.set_quantity('elevation', lambda x,y: -x-y)
1773        self.domain.set_quantity('stage', 0.0)
1774
1775        self.domain.geo_reference = Geo_reference(56,308500,6189000)
1776
1777
1778        sww = get_dataobject(self.domain)
1779        sww.store_connectivity()
1780        sww.store_timestep('stage')
1781
1782        self.domain.evolve_to_end(finaltime = 0.01)
1783        sww.store_timestep('stage')
1784
1785        cellsize = 0.25
1786        #Check contents
1787        #Get NetCDF
1788
1789        fid = NetCDFFile(sww.filename, 'r')
1790
1791        # Get the variables
1792        x = fid.variables['x'][:]
1793        y = fid.variables['y'][:]
1794        z = fid.variables['elevation'][:]
1795        time = fid.variables['time'][:]
1796        stage = fid.variables['stage'][:]
1797
1798
1799        #Export to ascii/prj files
1800        sww2dem(self.domain.filename,
1801                basename_out = 'datatest_depth',
1802                quantity = 'stage - elevation',
1803                cellsize = cellsize,
1804                reduction = min,
1805                format = 'asc',
1806                verbose = False)
1807
1808
1809        #Check asc file
1810        ascid = open(ascfile)
1811        lines = ascid.readlines()
1812        ascid.close()
1813
1814        L = lines[0].strip().split()
1815        assert L[0].strip().lower() == 'ncols'
1816        assert L[1].strip().lower() == '5'
1817
1818        L = lines[1].strip().split()
1819        assert L[0].strip().lower() == 'nrows'
1820        assert L[1].strip().lower() == '5'
1821
1822        L = lines[2].strip().split()
1823        assert L[0].strip().lower() == 'xllcorner'
1824        assert allclose(float(L[1].strip().lower()), 308500)
1825
1826        L = lines[3].strip().split()
1827        assert L[0].strip().lower() == 'yllcorner'
1828        assert allclose(float(L[1].strip().lower()), 6189000)
1829
1830        L = lines[4].strip().split()
1831        assert L[0].strip().lower() == 'cellsize'
1832        assert allclose(float(L[1].strip().lower()), cellsize)
1833
1834        L = lines[5].strip().split()
1835        assert L[0].strip() == 'NODATA_value'
1836        assert L[1].strip().lower() == '-9999'
1837
1838
1839        #Check grid values (where applicable)
1840        for j in range(5):
1841            if j%2 == 0:
1842                L = lines[6+j].strip().split()
1843                jj = 4-j
1844                for i in range(5):
1845                    if i%2 == 0:
1846                        index = jj/2 + i/2*3
1847                        val0 = stage[0,index] - z[index]
1848                        val1 = stage[1,index] - z[index]
1849
1850                        #print i, j, index, ':', L[i], val0, val1
1851                        assert allclose(float(L[i]), min(val0, val1))
1852
1853
1854        fid.close()
1855
1856        #Cleanup
1857        os.remove(prjfile)
1858        os.remove(ascfile)
1859        #os.remove(swwfile)
1860
1861
1862
1863
1864
1865    def test_sww2dem_asc_missing_points(self):
1866        """Test that sww information can be converted correctly to asc/prj
1867        format readable by e.g. ArcView
1868
1869        This test includes the writing of missing values
1870        """
1871
1872        import time, os
1873        from Numeric import array, zeros, allclose, Float, concatenate
1874        from Scientific.IO.NetCDF import NetCDFFile
1875
1876        #Setup mesh not coinciding with rectangle.
1877        #This will cause missing values to occur in gridded data
1878
1879
1880        points = [                        [1.0, 1.0],
1881                              [0.5, 0.5], [1.0, 0.5],
1882                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
1883
1884        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
1885
1886        #Create shallow water domain
1887        domain = Domain(points, vertices)
1888        domain.default_order=2
1889
1890
1891        #Set some field values
1892        domain.set_quantity('elevation', lambda x,y: -x-y)
1893        domain.set_quantity('friction', 0.03)
1894
1895
1896        ######################
1897        # Boundary conditions
1898        B = Transmissive_boundary(domain)
1899        domain.set_boundary( {'exterior': B} )
1900
1901
1902        ######################
1903        #Initial condition - with jumps
1904
1905        bed = domain.quantities['elevation'].vertex_values
1906        stage = zeros(bed.shape, Float)
1907
1908        h = 0.3
1909        for i in range(stage.shape[0]):
1910            if i % 2 == 0:
1911                stage[i,:] = bed[i,:] + h
1912            else:
1913                stage[i,:] = bed[i,:]
1914
1915        domain.set_quantity('stage', stage)
1916        domain.distribute_to_vertices_and_edges()
1917
1918        domain.filename = 'datatest'
1919
1920        prjfile = domain.filename + '_elevation.prj'
1921        ascfile = domain.filename + '_elevation.asc'
1922        swwfile = domain.filename + '.sww'
1923
1924        domain.set_datadir('.')
1925        domain.format = 'sww'
1926        domain.smooth = True
1927
1928        domain.geo_reference = Geo_reference(56,308500,6189000)
1929
1930        sww = get_dataobject(domain)
1931        sww.store_connectivity()
1932        sww.store_timestep('stage')
1933
1934        cellsize = 0.25
1935        #Check contents
1936        #Get NetCDF
1937
1938        fid = NetCDFFile(swwfile, 'r')
1939
1940        # Get the variables
1941        x = fid.variables['x'][:]
1942        y = fid.variables['y'][:]
1943        z = fid.variables['elevation'][:]
1944        time = fid.variables['time'][:]
1945
1946        try:
1947            geo_reference = Geo_reference(NetCDFObject=fid)
1948        except AttributeError, e:
1949            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
1950
1951        #Export to ascii/prj files
1952        sww2dem(domain.filename,
1953                quantity = 'elevation',
1954                cellsize = cellsize,
1955                verbose = False,
1956                format = 'asc')
1957
1958
1959        #Check asc file
1960        ascid = open(ascfile)
1961        lines = ascid.readlines()
1962        ascid.close()
1963
1964        L = lines[0].strip().split()
1965        assert L[0].strip().lower() == 'ncols'
1966        assert L[1].strip().lower() == '5'
1967
1968        L = lines[1].strip().split()
1969        assert L[0].strip().lower() == 'nrows'
1970        assert L[1].strip().lower() == '5'
1971
1972        L = lines[2].strip().split()
1973        assert L[0].strip().lower() == 'xllcorner'
1974        assert allclose(float(L[1].strip().lower()), 308500)
1975
1976        L = lines[3].strip().split()
1977        assert L[0].strip().lower() == 'yllcorner'
1978        assert allclose(float(L[1].strip().lower()), 6189000)
1979
1980        L = lines[4].strip().split()
1981        assert L[0].strip().lower() == 'cellsize'
1982        assert allclose(float(L[1].strip().lower()), cellsize)
1983
1984        L = lines[5].strip().split()
1985        assert L[0].strip() == 'NODATA_value'
1986        assert L[1].strip().lower() == '-9999'
1987
1988        #Check grid values
1989        for j in range(5):
1990            L = lines[6+j].strip().split()
1991            assert len(L) == 5
1992            y = (4-j) * cellsize
1993
1994            for i in range(5):
1995                #print i
1996                if i+j >= 4:
1997                    assert allclose(float(L[i]), -i*cellsize - y)
1998                else:
1999                    #Missing values
2000                    assert allclose(float(L[i]), -9999)
2001
2002
2003
2004        fid.close()
2005
2006        #Cleanup
2007        os.remove(prjfile)
2008        os.remove(ascfile)
2009        os.remove(swwfile)
2010
2011    def test_sww2ers_simple(self):
2012        """Test that sww information can be converted correctly to asc/prj
2013        format readable by e.g. ArcView
2014        """
2015
2016        import time, os
2017        from Numeric import array, zeros, allclose, Float, concatenate
2018        from Scientific.IO.NetCDF import NetCDFFile
2019
2020
2021        NODATA_value = 1758323
2022
2023        #Setup
2024        self.domain.filename = 'datatest'
2025
2026        headerfile = self.domain.filename + '.ers'
2027        swwfile = self.domain.filename + '.sww'
2028
2029        self.domain.set_datadir('.')
2030        self.domain.format = 'sww'
2031        self.domain.smooth = True
2032        self.domain.set_quantity('elevation', lambda x,y: -x-y)
2033
2034        self.domain.geo_reference = Geo_reference(56,308500,6189000)
2035
2036        sww = get_dataobject(self.domain)
2037        sww.store_connectivity()
2038        sww.store_timestep('stage')
2039
2040        self.domain.evolve_to_end(finaltime = 0.01)
2041        sww.store_timestep('stage')
2042
2043        cellsize = 0.25
2044        #Check contents
2045        #Get NetCDF
2046
2047        fid = NetCDFFile(sww.filename, 'r')
2048
2049        # Get the variables
2050        x = fid.variables['x'][:]
2051        y = fid.variables['y'][:]
2052        z = fid.variables['elevation'][:]
2053        time = fid.variables['time'][:]
2054        stage = fid.variables['stage'][:]
2055
2056
2057        #Export to ers files
2058        #sww2ers(self.domain.filename,
2059        #        quantity = 'elevation',
2060        #        cellsize = cellsize,
2061        #        verbose = False)
2062               
2063        sww2dem(self.domain.filename,
2064                quantity = 'elevation',
2065                cellsize = cellsize,
2066                NODATA_value = NODATA_value,
2067                verbose = False,
2068                format = 'ers')
2069
2070        #Check header data
2071        from ermapper_grids import read_ermapper_header, read_ermapper_data
2072       
2073        header = read_ermapper_header(self.domain.filename + '_elevation.ers')
2074        #print header
2075        assert header['projection'].lower() == '"utm-56"'
2076        assert header['datum'].lower() == '"wgs84"'
2077        assert header['units'].lower() == '"meters"'   
2078        assert header['value'].lower() == '"elevation"'         
2079        assert header['xdimension'] == '0.25'
2080        assert header['ydimension'] == '0.25'   
2081        assert float(header['eastings']) == 308500.0   #xllcorner
2082        assert float(header['northings']) == 6189000.0 #yllcorner       
2083        assert int(header['nroflines']) == 5
2084        assert int(header['nrofcellsperline']) == 5     
2085        assert int(header['nullcellvalue']) == NODATA_value
2086        #FIXME - there is more in the header                   
2087
2088               
2089        #Check grid data               
2090        grid = read_ermapper_data(self.domain.filename + '_elevation') 
2091       
2092        #FIXME (Ole): Why is this the desired reference grid for -x-y?
2093        ref_grid = [NODATA_value, NODATA_value, NODATA_value, NODATA_value, NODATA_value,
2094                    -1,    -1.25, -1.5,  -1.75, -2.0,
2095                    -0.75, -1.0,  -1.25, -1.5,  -1.75,             
2096                    -0.5,  -0.75, -1.0,  -1.25, -1.5,
2097                    -0.25, -0.5,  -0.75, -1.0,  -1.25]             
2098                                         
2099
2100        #print grid
2101        assert allclose(grid, ref_grid)
2102
2103        fid.close()
2104       
2105        #Cleanup
2106        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
2107        #Done (Ole) - it was because sww2ers didn't close it's sww file
2108        os.remove(sww.filename)
2109        os.remove(self.domain.filename + '_elevation')
2110        os.remove(self.domain.filename + '_elevation.ers')                 
2111
2112
2113
2114    def test_ferret2sww1(self):
2115        """Test that georeferencing etc works when converting from
2116        ferret format (lat/lon) to sww format (UTM)
2117        """
2118        from Scientific.IO.NetCDF import NetCDFFile
2119        import os, sys
2120
2121        #The test file has
2122        # LON = 150.66667, 150.83334, 151, 151.16667
2123        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2124        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2125        #
2126        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2127        # Fourth value (index==3) is -6.50198 cm
2128
2129       
2130
2131        #Read
2132        from coordinate_transforms.redfearn import redfearn
2133        #fid = NetCDFFile(self.test_MOST_file)
2134        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')       
2135        first_value = fid.variables['HA'][:][0,0,0]
2136        fourth_value = fid.variables['HA'][:][0,0,3]
2137        fid.close()
2138
2139
2140        #Call conversion (with zero origin)
2141        #ferret2sww('small', verbose=False,
2142        #           origin = (56, 0, 0))
2143        ferret2sww(self.test_MOST_file, verbose=False,
2144                   origin = (56, 0, 0))
2145
2146        #Work out the UTM coordinates for first point
2147        zone, e, n = redfearn(-34.5, 150.66667)
2148        #print zone, e, n
2149
2150        #Read output file 'small.sww'
2151        #fid = NetCDFFile('small.sww')
2152        fid = NetCDFFile(self.test_MOST_file + '.sww')       
2153
2154        x = fid.variables['x'][:]
2155        y = fid.variables['y'][:]
2156
2157        #Check that first coordinate is correctly represented
2158        assert allclose(x[0], e)
2159        assert allclose(y[0], n)
2160
2161        #Check first value
2162        stage = fid.variables['stage'][:]
2163        xmomentum = fid.variables['xmomentum'][:]
2164        ymomentum = fid.variables['ymomentum'][:]
2165
2166        #print ymomentum
2167
2168        assert allclose(stage[0,0], first_value/100)  #Meters
2169
2170        #Check fourth value
2171        assert allclose(stage[0,3], fourth_value/100)  #Meters
2172
2173        fid.close()
2174
2175        #Cleanup
2176        import os
2177        os.remove(self.test_MOST_file + '.sww')
2178
2179
2180    def test_ferret2sww_2(self):
2181        """Test that georeferencing etc works when converting from
2182        ferret format (lat/lon) to sww format (UTM)
2183        """
2184        from Scientific.IO.NetCDF import NetCDFFile
2185
2186        #The test file has
2187        # LON = 150.66667, 150.83334, 151, 151.16667
2188        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2189        # TIME = 0, 0.1, 0.6, 1.1, 1.6, 2.1 ;
2190        #
2191        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2192        # Fourth value (index==3) is -6.50198 cm
2193
2194
2195        from coordinate_transforms.redfearn import redfearn
2196
2197        #fid = NetCDFFile('small_ha.nc')
2198        fid = NetCDFFile(self.test_MOST_file + '_ha.nc')
2199
2200        #Pick a coordinate and a value
2201
2202        time_index = 1
2203        lat_index = 0
2204        lon_index = 2
2205
2206        test_value = fid.variables['HA'][:][time_index, lat_index, lon_index]
2207        test_time = fid.variables['TIME'][:][time_index]
2208        test_lat = fid.variables['LAT'][:][lat_index]
2209        test_lon = fid.variables['LON'][:][lon_index]
2210
2211        linear_point_index = lat_index*4 + lon_index
2212        fid.close()
2213
2214        #Call conversion (with zero origin)
2215        ferret2sww(self.test_MOST_file, verbose=False,
2216                   origin = (56, 0, 0))
2217
2218
2219        #Work out the UTM coordinates for test point
2220        zone, e, n = redfearn(test_lat, test_lon)
2221
2222        #Read output file 'small.sww'
2223        fid = NetCDFFile(self.test_MOST_file + '.sww')
2224
2225        x = fid.variables['x'][:]
2226        y = fid.variables['y'][:]
2227
2228        #Check that test coordinate is correctly represented
2229        assert allclose(x[linear_point_index], e)
2230        assert allclose(y[linear_point_index], n)
2231
2232        #Check test value
2233        stage = fid.variables['stage'][:]
2234
2235        assert allclose(stage[time_index, linear_point_index], test_value/100)
2236
2237        fid.close()
2238
2239        #Cleanup
2240        import os
2241        os.remove(self.test_MOST_file + '.sww')
2242
2243
2244
2245    def test_ferret2sww3(self):
2246        """Elevation included
2247        """
2248        from Scientific.IO.NetCDF import NetCDFFile
2249
2250        #The test file has
2251        # LON = 150.66667, 150.83334, 151, 151.16667
2252        # LAT = -34.5, -34.33333, -34.16667, -34 ;
2253        # ELEVATION = [-1 -2 -3 -4
2254        #             -5 -6 -7 -8
2255        #              ...
2256        #              ...      -16]
2257        # where the top left corner is -1m,
2258        # and the ll corner is -13.0m
2259        #
2260        # First value (index=0) in small_ha.nc is 0.3400644 cm,
2261        # Fourth value (index==3) is -6.50198 cm
2262
2263        from coordinate_transforms.redfearn import redfearn
2264        import os
2265        fid1 = NetCDFFile('test_ha.nc','w')
2266        fid2 = NetCDFFile('test_ua.nc','w')
2267        fid3 = NetCDFFile('test_va.nc','w')
2268        fid4 = NetCDFFile('test_e.nc','w')
2269
2270        h1_list = [150.66667,150.83334,151.]
2271        h2_list = [-34.5,-34.33333]
2272
2273        long_name = 'LON'
2274        lat_name = 'LAT'
2275
2276        nx = 3
2277        ny = 2
2278
2279        for fid in [fid1,fid2,fid3]:
2280            fid.createDimension(long_name,nx)
2281            fid.createVariable(long_name,'d',(long_name,))
2282            fid.variables[long_name].point_spacing='uneven'
2283            fid.variables[long_name].units='degrees_east'
2284            fid.variables[long_name].assignValue(h1_list)
2285
2286            fid.createDimension(lat_name,ny)
2287            fid.createVariable(lat_name,'d',(lat_name,))
2288            fid.variables[lat_name].point_spacing='uneven'
2289            fid.variables[lat_name].units='degrees_north'
2290            fid.variables[lat_name].assignValue(h2_list)
2291
2292            fid.createDimension('TIME',2)
2293            fid.createVariable('TIME','d',('TIME',))
2294            fid.variables['TIME'].point_spacing='uneven'
2295            fid.variables['TIME'].units='seconds'
2296            fid.variables['TIME'].assignValue([0.,1.])
2297            if fid == fid3: break
2298
2299
2300        for fid in [fid4]:
2301            fid.createDimension(long_name,nx)
2302            fid.createVariable(long_name,'d',(long_name,))
2303            fid.variables[long_name].point_spacing='uneven'
2304            fid.variables[long_name].units='degrees_east'
2305            fid.variables[long_name].assignValue(h1_list)
2306
2307            fid.createDimension(lat_name,ny)
2308            fid.createVariable(lat_name,'d',(lat_name,))
2309            fid.variables[lat_name].point_spacing='uneven'
2310            fid.variables[lat_name].units='degrees_north'
2311            fid.variables[lat_name].assignValue(h2_list)
2312
2313        name = {}
2314        name[fid1]='HA'
2315        name[fid2]='UA'
2316        name[fid3]='VA'
2317        name[fid4]='ELEVATION'
2318
2319        units = {}
2320        units[fid1]='cm'
2321        units[fid2]='cm/s'
2322        units[fid3]='cm/s'
2323        units[fid4]='m'
2324
2325        values = {}
2326        values[fid1]=[[[5., 10.,15.], [13.,18.,23.]],[[50.,100.,150.],[130.,180.,230.]]]
2327        values[fid2]=[[[1., 2.,3.], [4.,5.,6.]],[[7.,8.,9.],[10.,11.,12.]]]
2328        values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
2329        values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
2330
2331        for fid in [fid1,fid2,fid3]:
2332          fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
2333          fid.variables[name[fid]].point_spacing='uneven'
2334          fid.variables[name[fid]].units=units[fid]
2335          fid.variables[name[fid]].assignValue(values[fid])
2336          fid.variables[name[fid]].missing_value = -99999999.
2337          if fid == fid3: break
2338
2339        for fid in [fid4]:
2340            fid.createVariable(name[fid],'d',(lat_name,long_name))
2341            fid.variables[name[fid]].point_spacing='uneven'
2342            fid.variables[name[fid]].units=units[fid]
2343            fid.variables[name[fid]].assignValue(values[fid])
2344            fid.variables[name[fid]].missing_value = -99999999.
2345
2346
2347        fid1.sync(); fid1.close()
2348        fid2.sync(); fid2.close()
2349        fid3.sync(); fid3.close()
2350        fid4.sync(); fid4.close()
2351
2352        fid1 = NetCDFFile('test_ha.nc','r')
2353        fid2 = NetCDFFile('test_e.nc','r')
2354        fid3 = NetCDFFile('test_va.nc','r')
2355
2356
2357        first_amp = fid1.variables['HA'][:][0,0,0]
2358        third_amp = fid1.variables['HA'][:][0,0,2]
2359        first_elevation = fid2.variables['ELEVATION'][0,0]
2360        third_elevation= fid2.variables['ELEVATION'][:][0,2]
2361        first_speed = fid3.variables['VA'][0,0,0]
2362        third_speed = fid3.variables['VA'][:][0,0,2]
2363
2364        fid1.close()
2365        fid2.close()
2366        fid3.close()
2367
2368        #Call conversion (with zero origin)
2369        ferret2sww('test', verbose=False,
2370                   origin = (56, 0, 0))
2371
2372        os.remove('test_va.nc')
2373        os.remove('test_ua.nc')
2374        os.remove('test_ha.nc')
2375        os.remove('test_e.nc')
2376
2377        #Read output file 'test.sww'
2378        fid = NetCDFFile('test.sww')
2379
2380
2381        #Check first value
2382        elevation = fid.variables['elevation'][:]
2383        stage = fid.variables['stage'][:]
2384        xmomentum = fid.variables['xmomentum'][:]
2385        ymomentum = fid.variables['ymomentum'][:]
2386
2387        #print ymomentum
2388        first_height = first_amp/100 - first_elevation
2389        third_height = third_amp/100 - third_elevation
2390        first_momentum=first_speed*first_height/100
2391        third_momentum=third_speed*third_height/100
2392
2393        assert allclose(ymomentum[0][0],first_momentum)  #Meters
2394        assert allclose(ymomentum[0][2],third_momentum)  #Meters
2395
2396        fid.close()
2397
2398        #Cleanup
2399        os.remove('test.sww')
2400
2401
2402
2403
2404    def test_ferret2sww_nz_origin(self):
2405        from Scientific.IO.NetCDF import NetCDFFile
2406        from coordinate_transforms.redfearn import redfearn
2407
2408        #Call conversion (with nonzero origin)
2409        ferret2sww(self.test_MOST_file, verbose=False,
2410                   origin = (56, 100000, 200000))
2411
2412
2413        #Work out the UTM coordinates for first point
2414        zone, e, n = redfearn(-34.5, 150.66667)
2415
2416        #Read output file 'small.sww'
2417        #fid = NetCDFFile('small.sww', 'r')
2418        fid = NetCDFFile(self.test_MOST_file + '.sww')
2419       
2420        x = fid.variables['x'][:]
2421        y = fid.variables['y'][:]
2422
2423        #Check that first coordinate is correctly represented
2424        assert allclose(x[0], e-100000)
2425        assert allclose(y[0], n-200000)
2426
2427        fid.close()
2428
2429        #Cleanup
2430        os.remove(self.test_MOST_file + '.sww')
2431
2432
2433
2434    def test_sww_extent(self):
2435        """Not a test, rather a look at the sww format
2436        """
2437
2438        import time, os
2439        from Numeric import array, zeros, allclose, Float, concatenate
2440        from Scientific.IO.NetCDF import NetCDFFile
2441
2442        self.domain.filename = 'datatest' + str(id(self))
2443        self.domain.format = 'sww'
2444        self.domain.smooth = True
2445        self.domain.reduction = mean
2446        self.domain.set_datadir('.')
2447
2448
2449        sww = get_dataobject(self.domain)
2450        sww.store_connectivity()
2451        sww.store_timestep('stage')
2452        self.domain.time = 2.
2453
2454        #Modify stage at second timestep
2455        stage = self.domain.quantities['stage'].vertex_values
2456        self.domain.set_quantity('stage', stage/2)
2457
2458        sww.store_timestep('stage')
2459
2460        file_and_extension_name = self.domain.filename + ".sww"
2461        #print "file_and_extension_name",file_and_extension_name
2462        [xmin, xmax, ymin, ymax, stagemin, stagemax] = \
2463               extent_sww(file_and_extension_name )
2464
2465        assert allclose(xmin, 0.0)
2466        assert allclose(xmax, 1.0)
2467        assert allclose(ymin, 0.0)
2468        assert allclose(ymax, 1.0)
2469        assert allclose(stagemin, -0.85)
2470        assert allclose(stagemax, 0.15)
2471
2472
2473        #Cleanup
2474        os.remove(sww.filename)
2475
2476
2477
2478    def test_sww2domain1(self):
2479        ################################################
2480        #Create a test domain, and evolve and save it.
2481        ################################################
2482        from mesh_factory import rectangular
2483        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2484             Constant_height, Time_boundary, Transmissive_boundary
2485        from Numeric import array
2486
2487        #Create basic mesh
2488
2489        yiel=0.01
2490        points, vertices, boundary = rectangular(10,10)
2491
2492        #Create shallow water domain
2493        domain = Domain(points, vertices, boundary)
2494        domain.geo_reference = Geo_reference(56,11,11)
2495        domain.smooth = False
2496        domain.visualise = False
2497        domain.store = True
2498        domain.filename = 'bedslope'
2499        domain.default_order=2
2500        #Bed-slope and friction
2501        domain.set_quantity('elevation', lambda x,y: -x/3)
2502        domain.set_quantity('friction', 0.1)
2503        # Boundary conditions
2504        from math import sin, pi
2505        Br = Reflective_boundary(domain)
2506        Bt = Transmissive_boundary(domain)
2507        Bd = Dirichlet_boundary([0.2,0.,0.])
2508        Bw = Time_boundary(domain=domain,
2509                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2510
2511        #domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2512        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2513
2514        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2515        #Initial condition
2516        h = 0.05
2517        elevation = domain.quantities['elevation'].vertex_values
2518        domain.set_quantity('stage', elevation + h)
2519
2520        domain.check_integrity()
2521        #Evolution
2522        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2523        #    domain.write_time()
2524            pass
2525
2526
2527        ##########################################
2528        #Import the example's file as a new domain
2529        ##########################################
2530        from anuga.pyvolution.data_manager import sww2domain
2531        from Numeric import allclose
2532        import os
2533
2534        filename = domain.datadir+os.sep+domain.filename+'.sww'
2535        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2536        #points, vertices, boundary = rectangular(15,15)
2537        #domain2.boundary = boundary
2538        ###################
2539        ##NOW TEST IT!!!
2540        ###################
2541
2542        #os.remove(domain.filename + '.sww')
2543        os.remove(filename)       
2544
2545        bits = ['vertex_coordinates']
2546        for quantity in ['elevation']+domain.quantities_to_be_stored:
2547            bits.append('get_quantity("%s").get_integral()' %quantity)
2548            bits.append('get_quantity("%s").get_values()' %quantity)
2549
2550        for bit in bits:
2551            #print 'testing that domain.'+bit+' has been restored'
2552            #print bit
2553        #print 'done'
2554            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2555
2556        ######################################
2557        #Now evolve them both, just to be sure
2558        ######################################x
2559        visualise = False
2560        #visualise = True
2561        domain.visualise = visualise
2562        domain.time = 0.
2563        from time import sleep
2564
2565        final = .1
2566        domain.set_quantity('friction', 0.1)
2567        domain.store = False
2568        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2569
2570
2571        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2572            if visualise: sleep(1.)
2573            #domain.write_time()
2574            pass
2575
2576        final = final - (domain2.starttime-domain.starttime)
2577        #BUT since domain1 gets time hacked back to 0:
2578        final = final + (domain2.starttime-domain.starttime)
2579
2580        domain2.smooth = False
2581        domain2.visualise = visualise
2582        domain2.store = False
2583        domain2.default_order=2
2584        domain2.set_quantity('friction', 0.1)
2585        #Bed-slope and friction
2586        # Boundary conditions
2587        Bd2=Dirichlet_boundary([0.2,0.,0.])
2588        domain2.boundary = domain.boundary
2589        #print 'domain2.boundary'
2590        #print domain2.boundary
2591        domain2.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2592        #domain2.set_boundary({'exterior': Bd})
2593
2594        domain2.check_integrity()
2595
2596        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2597            if visualise: sleep(1.)
2598            #domain2.write_time()
2599            pass
2600
2601        ###################
2602        ##NOW TEST IT!!!
2603        ##################
2604
2605        bits = [ 'vertex_coordinates']
2606
2607        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2608            bits.append('get_quantity("%s").get_integral()' %quantity)
2609            bits.append('get_quantity("%s").get_values()' %quantity)
2610
2611        for bit in bits:
2612            #print bit
2613            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2614
2615
2616    def test_sww2domain2(self):
2617        ##################################################################
2618        #Same as previous test, but this checks how NaNs are handled.
2619        ##################################################################
2620
2621
2622        from mesh_factory import rectangular
2623        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2624             Constant_height, Time_boundary, Transmissive_boundary
2625        from Numeric import array
2626
2627        #Create basic mesh
2628        points, vertices, boundary = rectangular(2,2)
2629
2630        #Create shallow water domain
2631        domain = Domain(points, vertices, boundary)
2632        domain.smooth = False
2633        domain.visualise = False
2634        domain.store = True
2635        domain.set_name('test_file')
2636        domain.set_datadir('.')
2637        domain.default_order=2
2638        domain.quantities_to_be_stored=['stage']
2639
2640        domain.set_quantity('elevation', lambda x,y: -x/3)
2641        domain.set_quantity('friction', 0.1)
2642
2643        from math import sin, pi
2644        Br = Reflective_boundary(domain)
2645        Bt = Transmissive_boundary(domain)
2646        Bd = Dirichlet_boundary([0.2,0.,0.])
2647        Bw = Time_boundary(domain=domain,
2648                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2649
2650        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
2651
2652        h = 0.05
2653        elevation = domain.quantities['elevation'].vertex_values
2654        domain.set_quantity('stage', elevation + h)
2655
2656        domain.check_integrity()
2657
2658        for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
2659            pass
2660            #domain.write_time()
2661
2662
2663
2664        ##################################
2665        #Import the file as a new domain
2666        ##################################
2667        from anuga.pyvolution.data_manager import sww2domain
2668        from Numeric import allclose
2669        import os
2670
2671        filename = domain.datadir+os.sep+domain.filename+'.sww'
2672
2673        #Fail because NaNs are present
2674        try:
2675            domain2 = sww2domain(filename,boundary,fail_if_NaN=True,verbose=False)
2676        except:
2677            #Now import it, filling NaNs to be 0
2678            filler = 0
2679            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
2680
2681        #Clean up
2682        os.remove(filename)
2683           
2684
2685        bits = [ 'geo_reference.get_xllcorner()',
2686                'geo_reference.get_yllcorner()',
2687                'vertex_coordinates']
2688
2689        for quantity in ['elevation']+domain.quantities_to_be_stored:
2690            bits.append('get_quantity("%s").get_integral()' %quantity)
2691            bits.append('get_quantity("%s").get_values()' %quantity)
2692
2693        for bit in bits:
2694        #    print 'testing that domain.'+bit+' has been restored'
2695            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2696
2697        assert max(max(domain2.get_quantity('xmomentum').get_values()))==filler
2698        assert min(min(domain2.get_quantity('xmomentum').get_values()))==filler
2699        assert max(max(domain2.get_quantity('ymomentum').get_values()))==filler
2700        assert min(min(domain2.get_quantity('ymomentum').get_values()))==filler
2701
2702
2703
2704    #def test_weed(self):
2705        from anuga.pyvolution.data_manager import weed
2706
2707        coordinates1 = [[0.,0.],[1.,0.],[1.,1.],[1.,0.],[2.,0.],[1.,1.]]
2708        volumes1 = [[0,1,2],[3,4,5]]
2709        boundary1= {(0,1): 'external',(1,2): 'not external',(2,0): 'external',(3,4): 'external',(4,5): 'external',(5,3): 'not external'}
2710        coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
2711
2712        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
2713
2714        assert len(points2)==len(coordinates2)
2715        for i in range(len(coordinates2)):
2716            coordinate = tuple(coordinates2[i])
2717            assert points2.has_key(coordinate)
2718            points2[coordinate]=i
2719
2720        for triangle in volumes1:
2721            for coordinate in triangle:
2722                assert coordinates2[points2[tuple(coordinates1[coordinate])]][0]==coordinates1[coordinate][0]
2723                assert coordinates2[points2[tuple(coordinates1[coordinate])]][1]==coordinates1[coordinate][1]
2724
2725
2726    #FIXME This fails - smooth makes the comparism too hard for allclose
2727    def ztest_sww2domain3(self):
2728        ################################################
2729        #DOMAIN.SMOOTH = TRUE !!!!!!!!!!!!!!!!!!!!!!!!!!
2730        ################################################
2731        from mesh_factory import rectangular
2732        from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
2733             Constant_height, Time_boundary, Transmissive_boundary
2734        from Numeric import array
2735        #Create basic mesh
2736
2737        yiel=0.01
2738        points, vertices, boundary = rectangular(10,10)
2739
2740        #Create shallow water domain
2741        domain = Domain(points, vertices, boundary)
2742        domain.geo_reference = Geo_reference(56,11,11)
2743        domain.smooth = True
2744        domain.visualise = False
2745        domain.store = True
2746        domain.filename = 'bedslope'
2747        domain.default_order=2
2748        #Bed-slope and friction
2749        domain.set_quantity('elevation', lambda x,y: -x/3)
2750        domain.set_quantity('friction', 0.1)
2751        # Boundary conditions
2752        from math import sin, pi
2753        Br = Reflective_boundary(domain)
2754        Bt = Transmissive_boundary(domain)
2755        Bd = Dirichlet_boundary([0.2,0.,0.])
2756        Bw = Time_boundary(domain=domain,
2757                           f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
2758
2759        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
2760
2761        domain.quantities_to_be_stored.extend(['xmomentum','ymomentum'])
2762        #Initial condition
2763        h = 0.05
2764        elevation = domain.quantities['elevation'].vertex_values
2765        domain.set_quantity('stage', elevation + h)
2766
2767
2768        domain.check_integrity()
2769        #Evolution
2770        for t in domain.evolve(yieldstep = yiel, finaltime = 0.05):
2771        #    domain.write_time()
2772            pass
2773
2774
2775        ##########################################
2776        #Import the example's file as a new domain
2777        ##########################################
2778        from anuga.pyvolution.data_manager import sww2domain
2779        from Numeric import allclose
2780        import os
2781
2782        filename = domain.datadir+os.sep+domain.filename+'.sww'
2783        domain2 = sww2domain(filename,None,fail_if_NaN=False,verbose = False)
2784        #points, vertices, boundary = rectangular(15,15)
2785        #domain2.boundary = boundary
2786        ###################
2787        ##NOW TEST IT!!!
2788        ###################
2789
2790        os.remove(domain.filename + '.sww')
2791
2792        #FIXME smooth domain so that they can be compared
2793
2794
2795        bits = []#'vertex_coordinates']
2796        for quantity in ['elevation']+domain.quantities_to_be_stored:
2797            bits.append('quantities["%s"].get_integral()'%quantity)
2798
2799
2800        for bit in bits:
2801            #print 'testing that domain.'+bit+' has been restored'
2802            #print bit
2803            #print 'done'
2804            #print ('domain.'+bit), eval('domain.'+bit)
2805            #print ('domain2.'+bit), eval('domain2.'+bit)
2806            assert allclose(eval('domain.'+bit),eval('domain2.'+bit),rtol=1.0e-1,atol=1.e-3)
2807            pass
2808
2809        ######################################
2810        #Now evolve them both, just to be sure
2811        ######################################x
2812        visualise = False
2813        visualise = True
2814        domain.visualise = visualise
2815        domain.time = 0.
2816        from time import sleep
2817
2818        final = .5
2819        domain.set_quantity('friction', 0.1)
2820        domain.store = False
2821        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Br})
2822
2823        for t in domain.evolve(yieldstep = yiel, finaltime = final):
2824            if visualise: sleep(.03)
2825            #domain.write_time()
2826            pass
2827
2828        domain2.smooth = True
2829        domain2.visualise = visualise
2830        domain2.store = False
2831        domain2.default_order=2
2832        domain2.set_quantity('friction', 0.1)
2833        #Bed-slope and friction
2834        # Boundary conditions
2835        Bd2=Dirichlet_boundary([0.2,0.,0.])
2836        Br2 = Reflective_boundary(domain2)
2837        domain2.boundary = domain.boundary
2838        #print 'domain2.boundary'
2839        #print domain2.boundary
2840        domain2.set_boundary({'left': Bd2, 'right': Bd2, 'top': Bd2, 'bottom': Br2})
2841        #domain2.boundary = domain.boundary
2842        #domain2.set_boundary({'exterior': Bd})
2843
2844        domain2.check_integrity()
2845
2846        for t in domain2.evolve(yieldstep = yiel, finaltime = final):
2847            if visualise: sleep(.03)
2848            #domain2.write_time()
2849            pass
2850
2851        ###################
2852        ##NOW TEST IT!!!
2853        ##################
2854
2855        print '><><><><>>'
2856        bits = [ 'vertex_coordinates']
2857
2858        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
2859            #bits.append('quantities["%s"].get_integral()'%quantity)
2860            bits.append('get_quantity("%s").get_values()' %quantity)
2861
2862        for bit in bits:
2863            print bit
2864            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
2865
2866
2867    def test_decimate_dem(self):
2868        """Test decimation of dem file
2869        """
2870
2871        import os
2872        from Numeric import ones, allclose, Float, arange
2873        from Scientific.IO.NetCDF import NetCDFFile
2874
2875        #Write test dem file
2876        root = 'decdemtest'
2877
2878        filename = root + '.dem'
2879        fid = NetCDFFile(filename, 'w')
2880
2881        fid.institution = 'Geoscience Australia'
2882        fid.description = 'NetCDF DEM format for compact and portable ' +\
2883                          'storage of spatial point data'
2884
2885        nrows = 15
2886        ncols = 18
2887
2888        fid.ncols = ncols
2889        fid.nrows = nrows
2890        fid.xllcorner = 2000.5
2891        fid.yllcorner = 3000.5
2892        fid.cellsize = 25
2893        fid.NODATA_value = -9999
2894
2895        fid.zone = 56
2896        fid.false_easting = 0.0
2897        fid.false_northing = 0.0
2898        fid.projection = 'UTM'
2899        fid.datum = 'WGS84'
2900        fid.units = 'METERS'
2901
2902        fid.createDimension('number_of_points', nrows*ncols)
2903
2904        fid.createVariable('elevation', Float, ('number_of_points',))
2905
2906        elevation = fid.variables['elevation']
2907
2908        elevation[:] = (arange(nrows*ncols))
2909
2910        fid.close()
2911
2912        #generate the elevation values expected in the decimated file
2913        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
2914                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
2915                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
2916                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
2917                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
2918                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
2919                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
2920                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
2921                         (144+145+146+162+163+164+180+181+182) / 9.0,
2922                         (148+149+150+166+167+168+184+185+186) / 9.0,
2923                         (152+153+154+170+171+172+188+189+190) / 9.0,
2924                         (156+157+158+174+175+176+192+193+194) / 9.0,
2925                         (216+217+218+234+235+236+252+253+254) / 9.0,
2926                         (220+221+222+238+239+240+256+257+258) / 9.0,
2927                         (224+225+226+242+243+244+260+261+262) / 9.0,
2928                         (228+229+230+246+247+248+264+265+266) / 9.0]
2929
2930        #generate a stencil for computing the decimated values
2931        stencil = ones((3,3), Float) / 9.0
2932
2933        decimate_dem(root, stencil=stencil, cellsize_new=100)
2934
2935        #Open decimated NetCDF file
2936        fid = NetCDFFile(root + '_100.dem', 'r')
2937
2938        # Get decimated elevation
2939        elevation = fid.variables['elevation']
2940
2941        #Check values
2942        assert allclose(elevation, ref_elevation)
2943
2944        #Cleanup
2945        fid.close()
2946
2947        os.remove(root + '.dem')
2948        os.remove(root + '_100.dem')
2949
2950    def test_decimate_dem_NODATA(self):
2951        """Test decimation of dem file that includes NODATA values
2952        """
2953
2954        import os
2955        from Numeric import ones, allclose, Float, arange, reshape
2956        from Scientific.IO.NetCDF import NetCDFFile
2957
2958        #Write test dem file
2959        root = 'decdemtest'
2960
2961        filename = root + '.dem'
2962        fid = NetCDFFile(filename, 'w')
2963
2964        fid.institution = 'Geoscience Australia'
2965        fid.description = 'NetCDF DEM format for compact and portable ' +\
2966                          'storage of spatial point data'
2967
2968        nrows = 15
2969        ncols = 18
2970        NODATA_value = -9999
2971
2972        fid.ncols = ncols
2973        fid.nrows = nrows
2974        fid.xllcorner = 2000.5
2975        fid.yllcorner = 3000.5
2976        fid.cellsize = 25
2977        fid.NODATA_value = NODATA_value
2978
2979        fid.zone = 56
2980        fid.false_easting = 0.0
2981        fid.false_northing = 0.0
2982        fid.projection = 'UTM'
2983        fid.datum = 'WGS84'
2984        fid.units = 'METERS'
2985
2986        fid.createDimension('number_of_points', nrows*ncols)
2987
2988        fid.createVariable('elevation', Float, ('number_of_points',))
2989
2990        elevation = fid.variables['elevation']
2991
2992        #generate initial elevation values
2993        elevation_tmp = (arange(nrows*ncols))
2994        #add some NODATA values
2995        elevation_tmp[0]   = NODATA_value
2996        elevation_tmp[95]  = NODATA_value
2997        elevation_tmp[188] = NODATA_value
2998        elevation_tmp[189] = NODATA_value
2999        elevation_tmp[190] = NODATA_value
3000        elevation_tmp[209] = NODATA_value
3001        elevation_tmp[252] = NODATA_value
3002
3003        elevation[:] = elevation_tmp
3004
3005        fid.close()
3006
3007        #generate the elevation values expected in the decimated file
3008        ref_elevation = [NODATA_value,
3009                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
3010                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
3011                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
3012                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
3013                         NODATA_value,
3014                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
3015                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
3016                         (144+145+146+162+163+164+180+181+182) / 9.0,
3017                         (148+149+150+166+167+168+184+185+186) / 9.0,
3018                         NODATA_value,
3019                         (156+157+158+174+175+176+192+193+194) / 9.0,
3020                         NODATA_value,
3021                         (220+221+222+238+239+240+256+257+258) / 9.0,
3022                         (224+225+226+242+243+244+260+261+262) / 9.0,
3023                         (228+229+230+246+247+248+264+265+266) / 9.0]
3024
3025        #generate a stencil for computing the decimated values
3026        stencil = ones((3,3), Float) / 9.0
3027
3028        decimate_dem(root, stencil=stencil, cellsize_new=100)
3029
3030        #Open decimated NetCDF file
3031        fid = NetCDFFile(root + '_100.dem', 'r')
3032
3033        # Get decimated elevation
3034        elevation = fid.variables['elevation']
3035
3036        #Check values
3037        assert allclose(elevation, ref_elevation)
3038
3039        #Cleanup
3040        fid.close()
3041
3042        os.remove(root + '.dem')
3043        os.remove(root + '_100.dem')
3044
3045    def xxxtestz_sww2ers_real(self):
3046        """Test that sww information can be converted correctly to asc/prj
3047        format readable by e.g. ArcView
3048        """
3049
3050        import time, os
3051        from Numeric import array, zeros, allclose, Float, concatenate
3052        from Scientific.IO.NetCDF import NetCDFFile
3053
3054        # the memory optimised least squares
3055        #  cellsize = 20,   # this one seems to hang
3056        #  cellsize = 200000, # Ran 1 test in 269.703s
3057                                #Ran 1 test in 267.344s
3058        #  cellsize = 20000,  # Ran 1 test in 460.922s
3059        #  cellsize = 2000   #Ran 1 test in 5340.250s
3060        #  cellsize = 200   #this one seems to hang, building matirx A
3061
3062        # not optimised
3063        # seems to hang
3064        #  cellsize = 2000   # Ran 1 test in 5334.563s
3065        #Export to ascii/prj files
3066        sww2dem('karratha_100m',
3067                quantity = 'depth',
3068                cellsize = 200000,
3069                verbose = True)
3070
3071    def test_read_asc(self):
3072        """Test conversion from dem in ascii format to native NetCDF xya format
3073        """
3074
3075        import time, os
3076        from Numeric import array, zeros, allclose, Float, concatenate
3077        from Scientific.IO.NetCDF import NetCDFFile
3078
3079        import anuga.pyvolution.data_manager
3080        #Write test asc file
3081        filename = tempfile.mktemp(".000")
3082        fid = open(filename, 'w')
3083        fid.write("""ncols         7
3084nrows         4
3085xllcorner     2000.5
3086yllcorner     3000.5
3087cellsize      25
3088NODATA_value  -9999
3089    97.921    99.285   125.588   180.830   258.645   342.872   415.836
3090   473.157   514.391   553.893   607.120   678.125   777.283   883.038
3091   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
3092   502.645   516.230   504.739   450.604   388.500   338.097   514.980
3093""")
3094        fid.close()
3095        bath_metadata, grid = \
3096                   data_manager._read_asc(filename, verbose=False)
3097        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
3098        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
3099        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
3100        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
3101        self.failUnless(grid[0][0]  == 97.921,  'Failed')
3102        self.failUnless(grid[3][6]  == 514.980,  'Failed')
3103
3104        os.remove(filename)
3105       
3106    def test_asc_csiro2sww(self):
3107        import tempfile
3108
3109        bath_dir = tempfile.mkdtemp()
3110        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3111        #bath_dir = 'bath_data_manager_test'
3112        #print "os.getcwd( )",os.getcwd( )
3113        elevation_dir =  tempfile.mkdtemp()
3114        #elevation_dir = 'elev_expanded'
3115        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3116        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3117       
3118        fid = open(bath_dir_filename, 'w')
3119        fid.write(""" ncols             3
3120 nrows             2
3121 xllcorner    148.00000
3122 yllcorner    -38.00000
3123 cellsize       0.25
3124 nodata_value   -9999.0
3125    9000.000    -1000.000    3000.0
3126   -1000.000    9000.000  -1000.000
3127""")
3128        fid.close()
3129       
3130        fid = open(elevation_dir_filename1, 'w')
3131        fid.write(""" ncols             3
3132 nrows             2
3133 xllcorner    148.00000
3134 yllcorner    -38.00000
3135 cellsize       0.25
3136 nodata_value   -9999.0
3137    9000.000    0.000    3000.0
3138     0.000     9000.000     0.000
3139""")
3140        fid.close()
3141
3142        fid = open(elevation_dir_filename2, 'w')
3143        fid.write(""" ncols             3
3144 nrows             2
3145 xllcorner    148.00000
3146 yllcorner    -38.00000
3147 cellsize       0.25
3148 nodata_value   -9999.0
3149    9000.000    4000.000    4000.0
3150    4000.000    9000.000    4000.000
3151""")
3152        fid.close()
3153         
3154        ucur_dir =  tempfile.mkdtemp()
3155        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3156        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3157       
3158        fid = open(ucur_dir_filename1, 'w')
3159        fid.write(""" ncols             3
3160 nrows             2
3161 xllcorner    148.00000
3162 yllcorner    -38.00000
3163 cellsize       0.25
3164 nodata_value   -9999.0
3165    90.000    60.000    30.0
3166    10.000    10.000    10.000
3167""")
3168        fid.close()
3169        fid = open(ucur_dir_filename2, 'w')
3170        fid.write(""" ncols             3
3171 nrows             2
3172 xllcorner    148.00000
3173 yllcorner    -38.00000
3174 cellsize       0.25
3175 nodata_value   -9999.0
3176    90.000    60.000    30.0
3177    10.000    10.000    10.000
3178""")
3179        fid.close()
3180       
3181        vcur_dir =  tempfile.mkdtemp()
3182        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3183        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3184       
3185        fid = open(vcur_dir_filename1, 'w')
3186        fid.write(""" ncols             3
3187 nrows             2
3188 xllcorner    148.00000
3189 yllcorner    -38.00000
3190 cellsize       0.25
3191 nodata_value   -9999.0
3192    90.000    60.000    30.0
3193    10.000    10.000    10.000
3194""")
3195        fid.close()
3196        fid = open(vcur_dir_filename2, 'w')
3197        fid.write(""" ncols             3
3198 nrows             2
3199 xllcorner    148.00000
3200 yllcorner    -38.00000
3201 cellsize       0.25
3202 nodata_value   -9999.0
3203    90.000    60.000    30.0
3204    10.000    10.000    10.000
3205""")
3206        fid.close()
3207       
3208        sww_file = 'a_test.sww'
3209        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
3210
3211        # check the sww file
3212       
3213        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3214        x = fid.variables['x'][:]
3215        y = fid.variables['y'][:]
3216        z = fid.variables['z'][:]
3217        stage = fid.variables['stage'][:]
3218        xmomentum = fid.variables['xmomentum'][:]
3219        geo_ref = Geo_reference(NetCDFObject=fid)
3220        #print "geo_ref",geo_ref
3221        x_ref = geo_ref.get_xllcorner()
3222        y_ref = geo_ref.get_yllcorner()
3223        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3224        assert allclose(x_ref, 587798.418) # (-38, 148)
3225        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3226       
3227        #Zone:   55   
3228        #Easting:  588095.674  Northing: 5821451.722
3229        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3230        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3231
3232        #Zone:   55   
3233        #Easting:  632145.632  Northing: 5820863.269
3234        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3235        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3236
3237        #Zone:   55   
3238        #Easting:  609748.788  Northing: 5793447.860
3239        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3240        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3241
3242        assert allclose(z[0],9000.0 )
3243        assert allclose(stage[0][1],0.0 )
3244
3245        #(4000+1000)*60
3246        assert allclose(xmomentum[1][1],300000.0 )
3247       
3248       
3249        fid.close()
3250   
3251        #tidy up
3252        os.remove(bath_dir_filename)
3253        os.rmdir(bath_dir)
3254
3255        os.remove(elevation_dir_filename1)
3256        os.remove(elevation_dir_filename2)
3257        os.rmdir(elevation_dir)
3258       
3259        os.remove(ucur_dir_filename1)
3260        os.remove(ucur_dir_filename2)
3261        os.rmdir(ucur_dir)
3262       
3263        os.remove(vcur_dir_filename1)
3264        os.remove(vcur_dir_filename2)
3265        os.rmdir(vcur_dir)
3266
3267
3268        # remove sww file
3269        os.remove(sww_file)
3270       
3271    def test_asc_csiro2sww2(self):
3272        import tempfile
3273
3274        bath_dir = tempfile.mkdtemp()
3275        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3276        #bath_dir = 'bath_data_manager_test'
3277        #print "os.getcwd( )",os.getcwd( )
3278        elevation_dir =  tempfile.mkdtemp()
3279        #elevation_dir = 'elev_expanded'
3280        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3281        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3282       
3283        fid = open(bath_dir_filename, 'w')
3284        fid.write(""" ncols             3
3285 nrows             2
3286 xllcorner    148.00000
3287 yllcorner    -38.00000
3288 cellsize       0.25
3289 nodata_value   -9999.0
3290    9000.000    -1000.000    3000.0
3291   -1000.000    9000.000  -1000.000
3292""")
3293        fid.close()
3294       
3295        fid = open(elevation_dir_filename1, 'w')
3296        fid.write(""" ncols             3
3297 nrows             2
3298 xllcorner    148.00000
3299 yllcorner    -38.00000
3300 cellsize       0.25
3301 nodata_value   -9999.0
3302    9000.000    0.000    3000.0
3303     0.000     -9999.000     -9999.000
3304""")
3305        fid.close()
3306
3307        fid = open(elevation_dir_filename2, 'w')
3308        fid.write(""" ncols             3
3309 nrows             2
3310 xllcorner    148.00000
3311 yllcorner    -38.00000
3312 cellsize       0.25
3313 nodata_value   -9999.0
3314    9000.000    4000.000    4000.0
3315    4000.000    9000.000    4000.000
3316""")
3317        fid.close()
3318         
3319        ucur_dir =  tempfile.mkdtemp()
3320        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3321        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3322       
3323        fid = open(ucur_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    90.000    60.000    30.0
3331    10.000    10.000    10.000
3332""")
3333        fid.close()
3334        fid = open(ucur_dir_filename2, 'w')
3335        fid.write(""" ncols             3
3336 nrows             2
3337 xllcorner    148.00000
3338 yllcorner    -38.00000
3339 cellsize       0.25
3340 nodata_value   -9999.0
3341    90.000    60.000    30.0
3342    10.000    10.000    10.000
3343""")
3344        fid.close()
3345       
3346        vcur_dir =  tempfile.mkdtemp()
3347        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3348        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3349       
3350        fid = open(vcur_dir_filename1, 'w')
3351        fid.write(""" ncols             3
3352 nrows             2
3353 xllcorner    148.00000
3354 yllcorner    -38.00000
3355 cellsize       0.25
3356 nodata_value   -9999.0
3357    90.000    60.000    30.0
3358    10.000    10.000    10.000
3359""")
3360        fid.close()
3361        fid = open(vcur_dir_filename2, 'w')
3362        fid.write(""" ncols             3
3363 nrows             2
3364 xllcorner    148.00000
3365 yllcorner    -38.00000
3366 cellsize       0.25
3367 nodata_value   -9999.0
3368    90.000    60.000    30.0
3369    10.000    10.000    10.000
3370""")
3371        fid.close()
3372       
3373        try:
3374            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
3375                          vcur_dir, sww_file)
3376        except:
3377            #tidy up
3378            os.remove(bath_dir_filename)
3379            os.rmdir(bath_dir)
3380
3381            os.remove(elevation_dir_filename1)
3382            os.remove(elevation_dir_filename2)
3383            os.rmdir(elevation_dir)
3384       
3385            os.remove(ucur_dir_filename1)
3386            os.remove(ucur_dir_filename2)
3387            os.rmdir(ucur_dir)
3388       
3389            os.remove(vcur_dir_filename1)
3390            os.remove(vcur_dir_filename2)
3391            os.rmdir(vcur_dir)
3392        else:
3393            #tidy up
3394            os.remove(bath_dir_filename)
3395            os.rmdir(bath_dir)
3396
3397            os.remove(elevation_dir_filename1)
3398            os.remove(elevation_dir_filename2)
3399            os.rmdir(elevation_dir)
3400            raise 'Should raise exception'
3401       
3402            os.remove(ucur_dir_filename1)
3403            os.remove(ucur_dir_filename2)
3404            os.rmdir(ucur_dir)
3405       
3406            os.remove(vcur_dir_filename1)
3407            os.remove(vcur_dir_filename2)
3408            os.rmdir(vcur_dir)
3409
3410       
3411     
3412    def test_asc_csiro2sww3(self):
3413        import tempfile
3414
3415        bath_dir = tempfile.mkdtemp()
3416        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3417        #bath_dir = 'bath_data_manager_test'
3418        #print "os.getcwd( )",os.getcwd( )
3419        elevation_dir =  tempfile.mkdtemp()
3420        #elevation_dir = 'elev_expanded'
3421        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3422        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3423       
3424        fid = open(bath_dir_filename, 'w')
3425        fid.write(""" ncols             3
3426 nrows             2
3427 xllcorner    148.00000
3428 yllcorner    -38.00000
3429 cellsize       0.25
3430 nodata_value   -9999.0
3431    9000.000    -1000.000    3000.0
3432   -1000.000    9000.000  -1000.000
3433""")
3434        fid.close()
3435       
3436        fid = open(elevation_dir_filename1, 'w')
3437        fid.write(""" ncols             3
3438 nrows             2
3439 xllcorner    148.00000
3440 yllcorner    -38.00000
3441 cellsize       0.25
3442 nodata_value   -9999.0
3443    9000.000    0.000    3000.0
3444     0.000     -9999.000     -9999.000
3445""")
3446        fid.close()
3447
3448        fid = open(elevation_dir_filename2, 'w')
3449        fid.write(""" ncols             3
3450 nrows             2
3451 xllcorner    148.00000
3452 yllcorner    -38.00000
3453 cellsize       0.25
3454 nodata_value   -9999.0
3455    9000.000    4000.000    4000.0
3456    4000.000    9000.000    4000.000
3457""")
3458        fid.close()
3459         
3460        ucur_dir =  tempfile.mkdtemp()
3461        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3462        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3463       
3464        fid = open(ucur_dir_filename1, 'w')
3465        fid.write(""" ncols             3
3466 nrows             2
3467 xllcorner    148.00000
3468 yllcorner    -38.00000
3469 cellsize       0.25
3470 nodata_value   -9999.0
3471    90.000    60.000    30.0
3472    10.000    10.000    10.000
3473""")
3474        fid.close()
3475        fid = open(ucur_dir_filename2, 'w')
3476        fid.write(""" ncols             3
3477 nrows             2
3478 xllcorner    148.00000
3479 yllcorner    -38.00000
3480 cellsize       0.25
3481 nodata_value   -9999.0
3482    90.000    60.000    30.0
3483    10.000    10.000    10.000
3484""")
3485        fid.close()
3486       
3487        vcur_dir =  tempfile.mkdtemp()
3488        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3489        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3490       
3491        fid = open(vcur_dir_filename1, 'w')
3492        fid.write(""" ncols             3
3493 nrows             2
3494 xllcorner    148.00000
3495 yllcorner    -38.00000
3496 cellsize       0.25
3497 nodata_value   -9999.0
3498    90.000    60.000    30.0
3499    10.000    10.000    10.000
3500""")
3501        fid.close()
3502        fid = open(vcur_dir_filename2, 'w')
3503        fid.write(""" ncols             3
3504 nrows             2
3505 xllcorner    148.00000
3506 yllcorner    -38.00000
3507 cellsize       0.25
3508 nodata_value   -9999.0
3509    90.000    60.000    30.0
3510    10.000    10.000    10.000
3511""")
3512        fid.close()
3513       
3514        sww_file = 'a_test.sww'
3515        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3516                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3517                      mean_stage = 100)
3518
3519        # check the sww file
3520       
3521        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
3522        x = fid.variables['x'][:]
3523        y = fid.variables['y'][:]
3524        z = fid.variables['z'][:]
3525        stage = fid.variables['stage'][:]
3526        xmomentum = fid.variables['xmomentum'][:]
3527        geo_ref = Geo_reference(NetCDFObject=fid)
3528        #print "geo_ref",geo_ref
3529        x_ref = geo_ref.get_xllcorner()
3530        y_ref = geo_ref.get_yllcorner()
3531        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3532        assert allclose(x_ref, 587798.418) # (-38, 148)
3533        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
3534       
3535        #Zone:   55   
3536        #Easting:  588095.674  Northing: 5821451.722
3537        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
3538        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
3539
3540        #Zone:   55   
3541        #Easting:  632145.632  Northing: 5820863.269
3542        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3543        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
3544
3545        #Zone:   55   
3546        #Easting:  609748.788  Northing: 5793447.860
3547        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 ''
3548        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
3549
3550        assert allclose(z[0],9000.0 )
3551        assert allclose(stage[0][4],100.0 )
3552        assert allclose(stage[0][5],100.0 )
3553
3554        #(100.0 - 9000)*10
3555        assert allclose(xmomentum[0][4], -89000.0 )
3556       
3557        #(100.0 - -1000.000)*10
3558        assert allclose(xmomentum[0][5], 11000.0 )
3559       
3560        fid.close()
3561   
3562        #tidy up
3563        os.remove(bath_dir_filename)
3564        os.rmdir(bath_dir)
3565
3566        os.remove(elevation_dir_filename1)
3567        os.remove(elevation_dir_filename2)
3568        os.rmdir(elevation_dir)
3569       
3570        os.remove(ucur_dir_filename1)
3571        os.remove(ucur_dir_filename2)
3572        os.rmdir(ucur_dir)
3573       
3574        os.remove(vcur_dir_filename1)
3575        os.remove(vcur_dir_filename2)
3576        os.rmdir(vcur_dir)
3577
3578        # remove sww file
3579        os.remove(sww_file)
3580   
3581     
3582    def test_asc_csiro2sww4(self):
3583        """
3584        Test specifying the extent
3585        """
3586       
3587        import tempfile
3588
3589        bath_dir = tempfile.mkdtemp()
3590        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
3591        #bath_dir = 'bath_data_manager_test'
3592        #print "os.getcwd( )",os.getcwd( )
3593        elevation_dir =  tempfile.mkdtemp()
3594        #elevation_dir = 'elev_expanded'
3595        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
3596        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
3597       
3598        fid = open(bath_dir_filename, 'w')
3599        fid.write(""" ncols             4
3600 nrows             4
3601 xllcorner    148.00000
3602 yllcorner    -38.00000
3603 cellsize       0.25
3604 nodata_value   -9999.0
3605   -9000.000    -1000.000   -3000.0 -2000.000
3606   -1000.000    9000.000  -1000.000 -3000.000
3607   -4000.000    6000.000   2000.000 -5000.000
3608   -9000.000    -1000.000   -3000.0 -2000.000
3609""")
3610        fid.close()
3611       
3612        fid = open(elevation_dir_filename1, 'w')
3613        fid.write(""" ncols             4
3614 nrows             4
3615 xllcorner    148.00000
3616 yllcorner    -38.00000
3617 cellsize       0.25
3618 nodata_value   -9999.0
3619   -900.000    -100.000   -300.0 -200.000
3620   -100.000    900.000  -100.000 -300.000
3621   -400.000    600.000   200.000 -500.000
3622   -900.000    -100.000   -300.0 -200.000
3623""")
3624        fid.close()
3625
3626        fid = open(elevation_dir_filename2, 'w')
3627        fid.write(""" ncols             4
3628 nrows             4
3629 xllcorner    148.00000
3630 yllcorner    -38.00000
3631 cellsize       0.25
3632 nodata_value   -9999.0
3633   -990.000    -110.000   -330.0 -220.000
3634   -110.000    990.000  -110.000 -330.000
3635   -440.000    660.000   220.000 -550.000
3636   -990.000    -110.000   -330.0 -220.000
3637""")
3638        fid.close()
3639         
3640        ucur_dir =  tempfile.mkdtemp()
3641        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
3642        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
3643       
3644        fid = open(ucur_dir_filename1, 'w')
3645        fid.write(""" ncols             4
3646 nrows             4
3647 xllcorner    148.00000
3648 yllcorner    -38.00000
3649 cellsize       0.25
3650 nodata_value   -9999.0
3651   -90.000    -10.000   -30.0 -20.000
3652   -10.000    90.000  -10.000 -30.000
3653   -40.000    60.000   20.000 -50.000
3654   -90.000    -10.000   -30.0 -20.000
3655""")
3656        fid.close()
3657        fid = open(ucur_dir_filename2, 'w')
3658        fid.write(""" ncols             4
3659 nrows             4
3660 xllcorner    148.00000
3661 yllcorner    -38.00000
3662 cellsize       0.25
3663 nodata_value   -9999.0
3664   -90.000    -10.000   -30.0 -20.000
3665   -10.000    99.000  -11.000 -30.000
3666   -40.000    66.000   22.000 -50.000
3667   -90.000    -10.000   -30.0 -20.000
3668""")
3669        fid.close()
3670       
3671        vcur_dir =  tempfile.mkdtemp()
3672        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
3673        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
3674       
3675        fid = open(vcur_dir_filename1, 'w')
3676        fid.write(""" ncols             4
3677 nrows             4
3678 xllcorner    148.00000
3679 yllcorner    -38.00000
3680 cellsize       0.25
3681 nodata_value   -9999.0
3682   -90.000    -10.000   -30.0 -20.000
3683   -10.000    80.000  -20.000 -30.000
3684   -40.000    50.000   10.000 -50.000
3685   -90.000    -10.000   -30.0 -20.000
3686""")
3687        fid.close()
3688        fid = open(vcur_dir_filename2, 'w')
3689        fid.write(""" ncols             4
3690 nrows             4
3691 xllcorner    148.00000
3692 yllcorner    -38.00000
3693 cellsize       0.25
3694 nodata_value   -9999.0
3695   -90.000    -10.000   -30.0 -20.000
3696   -10.000    88.000  -22.000 -30.000
3697   -40.000    55.000   11.000 -50.000
3698   -90.000    -10.000   -30.0 -20.000
3699""")
3700        fid.close()
3701       
3702        sww_file = tempfile.mktemp(".sww")
3703        #sww_file = 'a_test.sww'
3704        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
3705                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
3706                      mean_stage = 100,
3707                       minlat = -37.6, maxlat = -37.6,
3708                  minlon = 148.3, maxlon = 148.3
3709                      #,verbose = True
3710                      )
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        ymomentum = fid.variables['ymomentum'][:]
3721        geo_ref = Geo_reference(NetCDFObject=fid)
3722        #print "geo_ref",geo_ref
3723        x_ref = geo_ref.get_xllcorner()
3724        y_ref = geo_ref.get_yllcorner()
3725        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
3726
3727        assert allclose(fid.starttime, 0.0) # (-37.45, 148.25)
3728        assert allclose(x_ref, 610120.388) # (-37.45, 148.25)
3729        assert allclose(y_ref,  5820863.269 )# (-37.45, 148.5)
3730
3731        #Easting:  632145.632  Northing: 5820863.269
3732        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3733       
3734        #print "x",x
3735        #print "y",y
3736        self.failUnless(len(x) == 4,'failed') # 2*2
3737        self.failUnless(len(x) == 4,'failed') # 2*2
3738       
3739        #Zone:   55
3740        #Easting:  632145.632  Northing: 5820863.269
3741        #Latitude:   -37 45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 ''
3742        # magic number - y is close enough for me.
3743        assert allclose(x[3], 632145.63 - x_ref)
3744        assert allclose(y[3], 5820863.269  - y_ref + 5.22155314684e-005)
3745
3746        assert allclose(z[0],9000.0 ) #z is elevation info
3747        #print "z",z
3748        # 2 time steps, 4 points
3749        self.failUnless(xmomentum.shape == (2,4), 'failed') 
3750        self.failUnless(ymomentum.shape == (2,4), 'failed') 
3751       
3752        #(100.0 - -1000.000)*10
3753        #assert allclose(xmomentum[0][5], 11000.0 )
3754       
3755        fid.close()
3756       
3757        # is the sww file readable?
3758        #Lets see if we can convert it to a dem!
3759        #print "sww_file",sww_file
3760        #dem_file = tempfile.mktemp(".dem")
3761        domain = sww2domain(sww_file) ###, dem_file)
3762        domain.check_integrity()
3763
3764        #tidy up
3765        os.remove(bath_dir_filename)
3766        os.rmdir(bath_dir)
3767
3768        os.remove(elevation_dir_filename1)
3769        os.remove(elevation_dir_filename2)
3770        os.rmdir(elevation_dir)
3771       
3772        os.remove(ucur_dir_filename1)
3773        os.remove(ucur_dir_filename2)
3774        os.rmdir(ucur_dir)
3775       
3776        os.remove(vcur_dir_filename1)
3777        os.remove(vcur_dir_filename2)
3778        os.rmdir(vcur_dir)
3779
3780
3781     
3782       
3783        # remove sww file
3784        os.remove(sww_file)
3785       
3786        # remove dem file
3787        #os.remove(dem_file)
3788
3789    def test_get_min_max_indexes(self):
3790        latitudes = [3,2,1,0]
3791        longitudes = [0,10,20,30]
3792
3793        # k - lat
3794        # l - lon
3795        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3796            latitudes,longitudes,
3797            -10,4,-10,31)
3798
3799        #print "kmin",kmin;print "kmax",kmax
3800        #print "lmin",lmin;print "lmax",lmax
3801        latitudes_new = latitudes[kmin:kmax]
3802        longitudes_news = longitudes[lmin:lmax]
3803        #print "latitudes_new", latitudes_new
3804        #print "longitudes_news",longitudes_news
3805        self.failUnless(latitudes == latitudes_new and \
3806                        longitudes == longitudes_news,
3807                         'failed')
3808       
3809        ## 2nd test
3810        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3811            latitudes,longitudes,
3812            0.5,2.5,5,25)
3813        #print "kmin",kmin;print "kmax",kmax
3814        #print "lmin",lmin;print "lmax",lmax
3815        latitudes_new = latitudes[kmin:kmax]
3816        longitudes_news = longitudes[lmin:lmax]
3817        #print "latitudes_new", latitudes_new
3818        #print "longitudes_news",longitudes_news
3819       
3820        self.failUnless(latitudes == latitudes_new and \
3821                        longitudes == longitudes_news,
3822                         'failed')
3823       
3824        ## 3rd test
3825        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(latitudes,
3826                                                                    longitudes,
3827                                                      1.1,1.9,12,17)
3828        #print "kmin",kmin;print "kmax",kmax
3829        #print "lmin",lmin;print "lmax",lmax
3830        latitudes_new = latitudes[kmin:kmax]
3831        longitudes_news = longitudes[lmin:lmax]
3832        #print "latitudes_new", latitudes_new
3833        #print "longitudes_news",longitudes_news
3834       
3835        self.failUnless(latitudes_new == [2, 1] and \
3836                        longitudes_news == [10, 20],
3837                         'failed')
3838
3839       
3840        ## 4th test
3841        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3842            latitudes,longitudes,
3843                                                      -0.1,1.9,-2,17)
3844        #print "kmin",kmin;print "kmax",kmax
3845        #print "lmin",lmin;print "lmax",lmax
3846        latitudes_new = latitudes[kmin:kmax]
3847        longitudes_news = longitudes[lmin:lmax]
3848        #print "latitudes_new", latitudes_new
3849        #print "longitudes_news",longitudes_news
3850       
3851        self.failUnless(latitudes_new == [2, 1, 0] and \
3852                        longitudes_news == [0, 10, 20],
3853                         'failed')
3854        ## 5th test
3855        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3856            latitudes,longitudes,
3857            0.1,1.9,2,17)
3858        #print "kmin",kmin;print "kmax",kmax
3859        #print "lmin",lmin;print "lmax",lmax
3860        latitudes_new = latitudes[kmin:kmax]
3861        longitudes_news = longitudes[lmin:lmax]
3862        #print "latitudes_new", latitudes_new
3863        #print "longitudes_news",longitudes_news
3864       
3865        self.failUnless(latitudes_new == [2, 1, 0] and \
3866                        longitudes_news == [0, 10, 20],
3867                         'failed')
3868       
3869        ## 6th test
3870       
3871        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3872            latitudes,longitudes,
3873            1.5,4,18,32)
3874        #print "kmin",kmin;print "kmax",kmax
3875        #print "lmin",lmin;print "lmax",lmax
3876        latitudes_new = latitudes[kmin:kmax]
3877        longitudes_news = longitudes[lmin:lmax]
3878        #print "latitudes_new", latitudes_new
3879        #print "longitudes_news",longitudes_news
3880       
3881        self.failUnless(latitudes_new == [3, 2, 1] and \
3882                        longitudes_news == [10, 20, 30],
3883                         'failed')
3884       
3885       
3886        ## 7th test
3887        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3888        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3889            latitudes,longitudes,
3890            1.5,1.5,15,15)
3891        #print "kmin",kmin;print "kmax",kmax
3892        #print "lmin",lmin;print "lmax",lmax
3893        latitudes_new = latitudes[kmin:kmax]
3894        longitudes_news = longitudes[lmin:lmax]
3895        m2d = m2d[kmin:kmax,lmin:lmax]
3896        #print "m2d", m2d
3897        #print "latitudes_new", latitudes_new
3898        #print "longitudes_news",longitudes_news
3899       
3900        self.failUnless(latitudes_new == [2, 1] and \
3901                        longitudes_news == [10, 20],
3902                         'failed')
3903       
3904        self.failUnless(m2d == [[5,6],[9,10]],
3905                         'failed')
3906       
3907    def test_get_min_max_indexes2(self):
3908        latitudes = [-30,-35,-40,-45]
3909        longitudes = [148,149,150,151]
3910
3911        m2d = array([[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15]])
3912       
3913        # k - lat
3914        # l - lon
3915        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3916            latitudes,longitudes,
3917            -37,-27,147,149.5)
3918
3919        #print "kmin",kmin;print "kmax",kmax
3920        #print "lmin",lmin;print "lmax",lmax
3921        #print "m2d", m2d
3922        #print "latitudes", latitudes
3923        #print "longitudes",longitudes
3924        #print "latitudes[kmax]", latitudes[kmax]
3925        latitudes_new = latitudes[kmin:kmax]
3926        longitudes_new = longitudes[lmin:lmax]
3927        m2d = m2d[kmin:kmax,lmin:lmax]
3928        #print "m2d", m2d
3929        #print "latitudes_new", latitudes_new
3930        #print "longitudes_new",longitudes_new
3931       
3932        self.failUnless(latitudes_new == [-30, -35, -40] and \
3933                        longitudes_new == [148, 149,150],
3934                         'failed')
3935        self.failUnless(m2d == [[0,1,2],[4,5,6],[8,9,10]],
3936                         'failed')
3937       
3938    def test_get_min_max_indexes3(self):
3939        latitudes = [-30,-35,-40,-45,-50,-55,-60]
3940        longitudes = [148,149,150,151]
3941
3942        # k - lat
3943        # l - lon
3944        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3945            latitudes,longitudes,
3946            -43,-37,148.5,149.5)
3947
3948       
3949        #print "kmin",kmin;print "kmax",kmax
3950        #print "lmin",lmin;print "lmax",lmax
3951        #print "latitudes", latitudes
3952        #print "longitudes",longitudes
3953        latitudes_new = latitudes[kmin:kmax]
3954        longitudes_news = longitudes[lmin:lmax]
3955        #print "latitudes_new", latitudes_new
3956        #print "longitudes_news",longitudes_news
3957       
3958        self.failUnless(latitudes_new == [-35, -40, -45] and \
3959                        longitudes_news == [148, 149,150],
3960                         'failed')
3961       
3962    def test_get_min_max_indexes4(self):
3963        latitudes = [-30,-35,-40,-45,-50,-55,-60]
3964        longitudes = [148,149,150,151]
3965
3966        # k - lat
3967        # l - lon
3968        kmin, kmax, lmin, lmax = data_manager._get_min_max_indexes(
3969            latitudes,longitudes)
3970       
3971        #print "kmin",kmin;print "kmax",kmax
3972        #print "lmin",lmin;print "lmax",lmax
3973        #print "latitudes", latitudes
3974        #print "longitudes",longitudes
3975        latitudes_new = latitudes[kmin:kmax]
3976        longitudes_news = longitudes[lmin:lmax]
3977        #print "latitudes_new", latitudes_new
3978        #print "longitudes_news",longitudes_news
3979       
3980        self.failUnless(latitudes_new == latitudes  and \
3981                        longitudes_news == longitudes,
3982                         'failed')
3983       
3984    def test_tsh2sww(self):
3985        import os
3986        import tempfile
3987       
3988        tsh_file = tempfile.mktemp(".tsh")
3989        file = open(tsh_file,"w")
3990        file.write("4 3 # <vertex #> <x> <y> [attributes]\n \
39910 0.0 0.0 0.0 0.0 0.01 \n \
39921 1.0 0.0 10.0 10.0 0.02  \n \
39932 0.0 1.0 0.0 10.0 0.03  \n \
39943 0.5 0.25 8.0 12.0 0.04  \n \
3995# Vert att title  \n \
3996elevation  \n \
3997stage  \n \
3998friction  \n \
39992 # <triangle #> [<vertex #>] [<neigbouring triangle #>]  \n\
40000 0 3 2 -1  -1  1 dsg\n\
40011 0 1 3 -1  0 -1   ole nielsen\n\
40024 # <segment #> <vertex #>  <vertex #> [boundary tag] \n\
40030 1 0 2 \n\
40041 0 2 3 \n\
40052 2 3 \n\
40063 3 1 1 \n\
40073 0 # <x> <y> [attributes] ...Mesh Vertices... \n \
40080 216.0 -86.0 \n \
40091 160.0 -167.0 \n \
40102 114.0 -91.0 \n \
40113 # <vertex #>  <vertex #> [boundary tag] ...Mesh Segments... \n \
40120 0 1 0 \n \
40131 1 2 0 \n \
40142 2 0 0 \n \
40150 # <x> <y> ...Mesh Holes... \n \
40160 # <x> <y> <attribute>...Mesh Regions... \n \
40170 # <x> <y> <attribute>...Mesh Regions, area... \n\
4018#Geo reference \n \
401956 \n \
4020140 \n \
4021120 \n")
4022        file.close()
4023
4024        #sww_file = tempfile.mktemp(".sww")
4025        #print "sww_file",sww_file
4026        #print "sww_file",tsh_file
4027        tsh2sww(tsh_file)
4028
4029        os.remove(tsh_file)
4030        os.remove(tsh_file[:-4] + '.sww')
4031#-------------------------------------------------------------
4032if __name__ == "__main__":
4033    #suite = unittest.makeSuite(Test_Data_Manager,'test_tsh2sww')
4034    suite = unittest.makeSuite(Test_Data_Manager,'test')
4035    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
4036    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')   
4037    #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
4038    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')
4039    #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem_NODATA')
4040    runner = unittest.TextTestRunner()
4041    runner.run(suite)
4042 
Note: See TracBrowser for help on using the repository browser.