source: trunk/anuga_core/source/anuga/file_conversion/test_dem2dem.py @ 7778

Last change on this file since 7778 was 7778, checked in by hudson, 14 years ago

Cleaned up unit tests to use API.

File size: 7.4 KB
Line 
1import unittest
2import copy
3import os
4import numpy as num
5
6from anuga.config import netcdf_float
7
8from anuga.coordinate_transforms.geo_reference import Geo_reference
9from anuga.geometry.polygon import is_inside_polygon
10from anuga.abstract_2d_finite_volumes.util import file_function
11from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
12from anuga.config import g
13
14from boundaries import Reflective_boundary, \
15            Field_boundary, Transmissive_momentum_set_stage_boundary, \
16            Transmissive_stage_zero_momentum_boundary
17from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
18     import Transmissive_boundary, Dirichlet_boundary, \
19            Time_boundary, File_boundary, AWI_boundary
20
21from anuga.file.sww import get_mesh_and_quantities_from_file
22           
23from shallow_water_domain import Domain
24
25from anuga.abstract_2d_finite_volumes.mesh_factory \
26        import rectangular_cross
27from sww_interrogate import get_maximum_inundation_elevation, \
28            get_maximum_inundation_location, get_maximum_inundation_data, \
29            get_flow_through_cross_section, get_energy_through_cross_section
30           
31           
32from dem2dem import dem2dem
33               
34
35class Test_Dem2Dem(unittest.TestCase):
36    def test_decimate_dem(self):
37        """Test decimation of dem file
38        """
39
40        import os
41        from Scientific.IO.NetCDF import NetCDFFile
42
43        #Write test dem file
44        root = 'decdemtest'
45
46        filename = root + '.dem'
47        fid = NetCDFFile(filename, netcdf_mode_w)
48
49        fid.institution = 'Geoscience Australia'
50        fid.description = 'NetCDF DEM format for compact and portable ' +\
51                          'storage of spatial point data'
52
53        nrows = 15
54        ncols = 18
55
56        fid.ncols = ncols
57        fid.nrows = nrows
58        fid.xllcorner = 2000.5
59        fid.yllcorner = 3000.5
60        fid.cellsize = 25
61        fid.NODATA_value = -9999
62
63        fid.zone = 56
64        fid.false_easting = 0.0
65        fid.false_northing = 0.0
66        fid.projection = 'UTM'
67        fid.datum = 'WGS84'
68        fid.units = 'METERS'
69
70        fid.createDimension('number_of_points', nrows*ncols)
71
72        fid.createVariable('elevation', netcdf_float, ('number_of_points',))
73
74        elevation = fid.variables['elevation']
75
76        elevation[:] = (num.arange(nrows*ncols))
77
78        fid.close()
79
80        #generate the elevation values expected in the decimated file
81        ref_elevation = [(  0+  1+  2+ 18+ 19+ 20+ 36+ 37+ 38) / 9.0,
82                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
83                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
84                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
85                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
86                         ( 76+ 77+ 78+ 94+ 95+ 96+112+113+114) / 9.0,
87                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
88                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
89                         (144+145+146+162+163+164+180+181+182) / 9.0,
90                         (148+149+150+166+167+168+184+185+186) / 9.0,
91                         (152+153+154+170+171+172+188+189+190) / 9.0,
92                         (156+157+158+174+175+176+192+193+194) / 9.0,
93                         (216+217+218+234+235+236+252+253+254) / 9.0,
94                         (220+221+222+238+239+240+256+257+258) / 9.0,
95                         (224+225+226+242+243+244+260+261+262) / 9.0,
96                         (228+229+230+246+247+248+264+265+266) / 9.0]
97
98        # generate a stencil for computing the decimated values
99        stencil = num.ones((3,3), num.float) / 9.0
100
101        dem2dem(root, stencil=stencil, cellsize_new=100)
102
103        # Open decimated NetCDF file
104        fid = NetCDFFile(root + '_100.dem', netcdf_mode_r)
105
106        # Get decimated elevation
107        elevation = fid.variables['elevation']
108
109        # Check values
110        assert num.allclose(elevation, ref_elevation)
111
112        # Cleanup
113        fid.close()
114
115        os.remove(root + '.dem')
116        os.remove(root + '_100.dem')
117
118    def test_decimate_dem_NODATA(self):
119        """Test decimation of dem file that includes NODATA values
120        """
121
122        import os
123        from Scientific.IO.NetCDF import NetCDFFile
124
125        # Write test dem file
126        root = 'decdemtest'
127
128        filename = root + '.dem'
129        fid = NetCDFFile(filename, netcdf_mode_w)
130
131        fid.institution = 'Geoscience Australia'
132        fid.description = 'NetCDF DEM format for compact and portable ' +\
133                          'storage of spatial point data'
134
135        nrows = 15
136        ncols = 18
137        NODATA_value = -9999
138
139        fid.ncols = ncols
140        fid.nrows = nrows
141        fid.xllcorner = 2000.5
142        fid.yllcorner = 3000.5
143        fid.cellsize = 25
144        fid.NODATA_value = NODATA_value
145
146        fid.zone = 56
147        fid.false_easting = 0.0
148        fid.false_northing = 0.0
149        fid.projection = 'UTM'
150        fid.datum = 'WGS84'
151        fid.units = 'METERS'
152
153        fid.createDimension('number_of_points', nrows*ncols)
154
155        fid.createVariable('elevation', netcdf_float, ('number_of_points',))
156
157        elevation = fid.variables['elevation']
158
159        # Generate initial elevation values
160        elevation_tmp = (num.arange(nrows*ncols))
161
162        # Add some NODATA values
163        elevation_tmp[0]   = NODATA_value
164        elevation_tmp[95]  = NODATA_value
165        elevation_tmp[188] = NODATA_value
166        elevation_tmp[189] = NODATA_value
167        elevation_tmp[190] = NODATA_value
168        elevation_tmp[209] = NODATA_value
169        elevation_tmp[252] = NODATA_value
170
171        elevation[:] = elevation_tmp
172
173        fid.close()
174
175        # Generate the elevation values expected in the decimated file
176        ref_elevation = [NODATA_value,
177                         (  4+  5+  6+ 22+ 23+ 24+ 40+ 41+ 42) / 9.0,
178                         (  8+  9+ 10+ 26+ 27+ 28+ 44+ 45+ 46) / 9.0,
179                         ( 12+ 13+ 14+ 30+ 31+ 32+ 48+ 49+ 50) / 9.0,
180                         ( 72+ 73+ 74+ 90+ 91+ 92+108+109+110) / 9.0,
181                         NODATA_value,
182                         ( 80+ 81+ 82+ 98+ 99+100+116+117+118) / 9.0,
183                         ( 84+ 85+ 86+102+103+104+120+121+122) / 9.0,
184                         (144+145+146+162+163+164+180+181+182) / 9.0,
185                         (148+149+150+166+167+168+184+185+186) / 9.0,
186                         NODATA_value,
187                         (156+157+158+174+175+176+192+193+194) / 9.0,
188                         NODATA_value,
189                         (220+221+222+238+239+240+256+257+258) / 9.0,
190                         (224+225+226+242+243+244+260+261+262) / 9.0,
191                         (228+229+230+246+247+248+264+265+266) / 9.0]
192
193        # Generate a stencil for computing the decimated values
194        stencil = num.ones((3,3), num.float) / 9.0
195
196        dem2dem(root, stencil=stencil, cellsize_new=100)
197
198        # Open decimated NetCDF file
199        fid = NetCDFFile(root + '_100.dem', netcdf_mode_r)
200
201        # Get decimated elevation
202        elevation = fid.variables['elevation']
203
204        # Check values
205        assert num.allclose(elevation, ref_elevation)
206
207        # Cleanup
208        fid.close()
209
210        os.remove(root + '.dem')
211        os.remove(root + '_100.dem')     
212
213#################################################################################
214
215if __name__ == "__main__":
216    suite = unittest.makeSuite(Test_Dem2Dem, 'test')
217    runner = unittest.TextTestRunner(verbosity=1)
218    runner.run(suite)
219       
Note: See TracBrowser for help on using the repository browser.