source: inundation/pyvolution/test_data_manager.py @ 2555

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

More testing of dem2pts (It still works :-)

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