source: trunk/anuga_core/source/anuga/file_conversion/test_dem2pts.py @ 8780

Last change on this file since 8780 was 8780, checked in by steve, 11 years ago

Some changes to allow netcdf4 use
File size: 10.4 KB
Line 
1import sys
2import unittest
3import numpy as num
4import copy
5import os
6
7# ANUGA modules
8from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
9                            netcdf_float
10
11from anuga.file.netcdf import NetCDFFile
12
13from dem2pts import dem2pts
14from asc2dem import asc2dem
15
16class Test_Dem2Pts(unittest.TestCase):
17    """ A suite of tests to test file conversion functions.
18        These tests are quite coarse-grained: converting a file
19        and checking that its headers and some of its contents
20        are correct.
21    """ 
22 
23    def test_dem2pts_bounding_box_v2(self):
24        """Test conversion from dem in ascii format to native NetCDF format
25        """
26
27        import time, os
28
29
30        #Write test asc file
31        root = 'demtest'
32
33        filename = root+'.asc'
34        fid = open(filename, 'w')
35        fid.write("""ncols         10
36nrows         10
37xllcorner     2000
38yllcorner     3000
39cellsize      1
40NODATA_value  -9999
41""")
42        #Create linear function
43        ref_points = []
44        ref_elevation = []
45        x0 = 2000
46        y = 3010
47        yvec = range(10)
48        xvec = range(10)
49        z = -1
50        for i in range(10):
51            y = y - 1
52            for j in range(10):
53                x = x0 + xvec[j]
54                z += 1
55                ref_points.append ([x,y])
56                ref_elevation.append(z)
57                fid.write('%f ' %z)
58            fid.write('\n')
59
60        fid.close()
61
62        #print 'sending pts', ref_points
63        #print 'sending elev', ref_elevation
64
65        #Write prj file with metadata
66        metafilename = root+'.prj'
67        fid = open(metafilename, 'w')
68
69
70        fid.write("""Projection UTM
71Zone 56
72Datum WGS84
73Zunits NO
74Units METERS
75Spheroid WGS84
76Xshift 0.0000000000
77Yshift 10000000.0000000000
78Parameters
79""")
80        fid.close()
81
82        #Convert to NetCDF pts
83        dem2pts(filename, easting_min=2002.0, easting_max=2007.0,
84                northing_min=3003.0, northing_max=3006.0,
85                verbose=False)
86
87        #Check contents
88        #Get NetCDF
89        fid = NetCDFFile(root+'.pts', netcdf_mode_r)
90
91        # Get the variables
92        #print fid.variables.keys()
93        points = fid.variables['points']
94        elevation = fid.variables['elevation']
95
96        #Check values
97        assert fid.xllcorner == 2002.0
98        assert fid.yllcorner == 3003.0
99
100        #create new reference points
101        newz = []
102        newz[0:5] = ref_elevation[32:38]
103        newz[6:11] = ref_elevation[42:48]
104        newz[12:17] = ref_elevation[52:58]
105        newz[18:23] = ref_elevation[62:68]
106        ref_elevation = []
107        ref_elevation = newz
108        ref_points = []
109        x0 = 2002
110        y = 3007
111        yvec = range(4)
112        xvec = range(6)
113        for i in range(4):
114            y = y - 1
115            ynew = y - 3003.0
116            for j in range(6):
117                x = x0 + xvec[j]
118                xnew = x - 2002.0
119                ref_points.append ([xnew,ynew]) #Relative point values
120
121        assert num.allclose(points, ref_points)
122
123        assert num.allclose(elevation, ref_elevation)
124
125        #Cleanup
126        fid.close()
127
128
129        os.remove(root + '.pts')
130        os.remove(root + '.dem')
131        os.remove(root + '.asc')
132        os.remove(root + '.prj')
133
134
135    def test_dem2pts_bounding_box_removeNullvalues_v2(self):
136        """Test conversion from dem in ascii format to native NetCDF format
137        """
138
139        import time, os
140
141        #Write test asc file
142        root = 'demtest2'
143
144        filename = root+'.asc'
145        fid = open(filename, 'w')
146        fid.write("""ncols         10
147nrows         10
148xllcorner     2000
149yllcorner     3000
150cellsize      1
151NODATA_value  -9999
152""")
153        #Create linear function
154        ref_points = []
155        ref_elevation = []
156        x0 = 2000
157        y = 3010
158        yvec = range(10)
159        xvec = range(10)
160        #z = range(100)
161        z = num.zeros(100, num.int)     #array default#
162        NODATA_value = -9999
163        count = -1
164        for i in range(10):
165            y = y - 1
166            for j in range(10):
167                x = x0 + xvec[j]
168                ref_points.append ([x,y])
169                count += 1
170                z[count] = (4*i - 3*j)%13
171                if j == 4: z[count] = NODATA_value #column inside clipping region
172                if j == 8: z[count] = NODATA_value #column outside clipping region
173                if i == 9: z[count] = NODATA_value #row outside clipping region
174                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
175                ref_elevation.append( z[count] )
176                fid.write('%f ' %z[count])
177            fid.write('\n')
178
179        fid.close()
180
181        #print 'sending elev', ref_elevation
182
183        #Write prj file with metadata
184        metafilename = root+'.prj'
185        fid = open(metafilename, 'w')
186
187
188        fid.write("""Projection UTM
189Zone 56
190Datum WGS84
191Zunits NO
192Units METERS
193Spheroid WGS84
194Xshift 0.0000000000
195Yshift 10000000.0000000000
196Parameters
197""")
198        fid.close()
199
200        #Convert to NetCDF pts
201        dem2pts(filename, easting_min=2002.0, easting_max=2007.0,
202                northing_min=3003.0, northing_max=3006.0)
203
204        #Check contents
205        #Get NetCDF
206        fid = NetCDFFile(root+'.pts', netcdf_mode_r)
207
208        # Get the variables
209        #print fid.variables.keys()
210        points = fid.variables['points']
211        elevation = fid.variables['elevation']
212
213        #Check values
214        assert fid.xllcorner == 2002.0
215        assert fid.yllcorner == 3003.0
216
217        #create new reference points
218        newz = num.zeros(19, num.int)       #array default#
219        newz[0:2] = ref_elevation[32:34]
220        newz[2:5] = ref_elevation[35:38]
221        newz[5:7] = ref_elevation[42:44]
222        newz[7] = ref_elevation[45]
223        newz[8] = ref_elevation[47]
224        newz[9:11] = ref_elevation[52:54]
225        newz[11:14] = ref_elevation[55:58]
226        newz[14:16] = ref_elevation[62:64]
227        newz[16:19] = ref_elevation[65:68]
228
229
230        ref_elevation = newz
231        ref_points = []
232        new_ref_points = []
233        x0 = 2002
234        y = 3007
235        yvec = range(4)
236        xvec = range(6)
237        for i in range(4):
238            y = y - 1
239            ynew = y - 3003.0
240            for j in range(6):
241                x = x0 + xvec[j]
242                xnew = x - 2002.0
243                if j <> 2 and (i<>1 or j<>4):
244                    ref_points.append([x,y])
245                    new_ref_points.append ([xnew,ynew])
246
247
248        assert num.allclose(points, new_ref_points)
249        assert num.allclose(elevation, ref_elevation)
250
251        #Cleanup
252        fid.close()
253
254
255        os.remove(root + '.pts')
256        os.remove(root + '.dem')
257        os.remove(root + '.asc')
258        os.remove(root + '.prj')
259
260
261    def test_dem2pts_bounding_box_removeNullvalues_v3(self):
262        """Test conversion from dem in ascii format to native NetCDF format
263        Check missing values on clipping boundary
264        """
265
266        import time, os
267
268        #Write test asc file
269        root = 'demtest3'
270
271        filename = root+'.asc'
272        fid = open(filename, 'w')
273        fid.write("""ncols         10
274nrows         10
275xllcorner     2000
276yllcorner     3000
277cellsize      1
278NODATA_value  -9999
279""")
280        #Create linear function
281        ref_points = []
282        ref_elevation = []
283        x0 = 2000
284        y = 3010
285        yvec = range(10)
286        xvec = range(10)
287        #z = range(100)
288        z = num.zeros(100, num.int)     #array default#
289        NODATA_value = -9999
290        count = -1
291        for i in range(10):
292            y = y - 1
293            for j in range(10):
294                x = x0 + xvec[j]
295                ref_points.append ([x,y])
296                count += 1
297                z[count] = (4*i - 3*j)%13
298                if j == 4: z[count] = NODATA_value #column inside clipping region
299                if j == 8: z[count] = NODATA_value #column outside clipping region
300                if i == 6: z[count] = NODATA_value #row on clipping boundary
301                if i == 4 and j == 6: z[count] = NODATA_value #arbitrary point inside clipping region
302                ref_elevation.append( z[count] )
303                fid.write('%f ' %z[count])
304            fid.write('\n')
305
306        fid.close()
307
308        #print 'sending elev', ref_elevation
309
310        #Write prj file with metadata
311        metafilename = root+'.prj'
312        fid = open(metafilename, 'w')
313
314
315        fid.write("""Projection UTM
316Zone 56
317Datum WGS84
318Zunits NO
319Units METERS
320Spheroid WGS84
321Xshift 0.0000000000
322Yshift 10000000.0000000000
323Parameters
324""")
325        fid.close()
326
327        #Convert to NetCDF pts
328        dem2pts(filename, easting_min=2002.0, easting_max=2007.0,
329                northing_min=3003.0, northing_max=3006.0)
330
331        #Check contents
332        #Get NetCDF
333        fid = NetCDFFile(root+'.pts', netcdf_mode_r)
334
335        # Get the variables
336        #print fid.variables.keys()
337        points = fid.variables['points']
338        elevation = fid.variables['elevation']
339
340        #Check values
341        assert fid.xllcorner == 2002.0
342        assert fid.yllcorner == 3003.0
343
344        #create new reference points
345        newz = num.zeros(14, num.int)       #array default#
346        newz[0:2] = ref_elevation[32:34]
347        newz[2:5] = ref_elevation[35:38]
348        newz[5:7] = ref_elevation[42:44]
349        newz[7] = ref_elevation[45]
350        newz[8] = ref_elevation[47]
351        newz[9:11] = ref_elevation[52:54]
352        newz[11:14] = ref_elevation[55:58]
353
354
355
356        ref_elevation = newz
357        ref_points = []
358        new_ref_points = []
359        x0 = 2002
360        y = 3007
361        yvec = range(4)
362        xvec = range(6)
363        for i in range(4):
364            y = y - 1
365            ynew = y - 3003.0
366            for j in range(6):
367                x = x0 + xvec[j]
368                xnew = x - 2002.0
369                if j <> 2 and (i<>1 or j<>4) and i<>3:
370                    ref_points.append([x,y])
371                    new_ref_points.append ([xnew,ynew])
372
373
374        #print points[:],points[:].shape
375        #print new_ref_points, len(new_ref_points)
376
377        assert num.allclose(elevation, ref_elevation)
378        assert num.allclose(points, new_ref_points)
379
380
381        #Cleanup
382        fid.close()
383
384
385        os.remove(root + '.pts')
386        os.remove(root + '.dem')
387        os.remove(root + '.asc')
388        os.remove(root + '.prj')
389
390
391#-------------------------------------------------------------
392
393if __name__ == "__main__":
394    suite = unittest.makeSuite(Test_Dem2Pts,'test')
395    runner = unittest.TextTestRunner()
396    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.