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

Last change on this file since 7814 was 7814, checked in by hudson, 13 years ago

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