1 | import unittest |
---|
2 | import copy |
---|
3 | import os |
---|
4 | import numpy as num |
---|
5 | |
---|
6 | from anuga.config import netcdf_float |
---|
7 | |
---|
8 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
9 | from anuga.geometry.polygon import is_inside_polygon |
---|
10 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
11 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
12 | from anuga.config import g |
---|
13 | |
---|
14 | from boundaries import Reflective_boundary, \ |
---|
15 | Field_boundary, Transmissive_momentum_set_stage_boundary, \ |
---|
16 | Transmissive_stage_zero_momentum_boundary |
---|
17 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
18 | import Transmissive_boundary, Dirichlet_boundary, \ |
---|
19 | Time_boundary, File_boundary, AWI_boundary |
---|
20 | |
---|
21 | from anuga.file.sww import get_mesh_and_quantities_from_file |
---|
22 | |
---|
23 | from shallow_water_domain import Domain |
---|
24 | |
---|
25 | from anuga.abstract_2d_finite_volumes.mesh_factory \ |
---|
26 | import rectangular_cross |
---|
27 | from 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 | |
---|
32 | from dem2dem import dem2dem |
---|
33 | |
---|
34 | |
---|
35 | class 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 | |
---|
215 | if __name__ == "__main__": |
---|
216 | suite = unittest.makeSuite(Test_Dem2Dem, 'test') |
---|
217 | runner = unittest.TextTestRunner(verbosity=1) |
---|
218 | runner.run(suite) |
---|
219 | |
---|