source: trunk/anuga_core/source/anuga/file_conversion/test_dem2dem.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: 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(filename, 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(filename, 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.