1 | """This scripts recreates fitting problem where values seem to go towards zero |
---|
2 | in areas with missing data. |
---|
3 | """ |
---|
4 | from os.path import join |
---|
5 | from anuga.utilities.polygon import read_polygon |
---|
6 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf |
---|
7 | from anuga.shallow_water.data_manager import dem2pts |
---|
8 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
9 | from anuga.interface import create_domain_from_regions |
---|
10 | from anuga.interface import Reflective_boundary |
---|
11 | |
---|
12 | import project, os |
---|
13 | import sys |
---|
14 | from anuga.shallow_water.data_manager import start_screen_catcher |
---|
15 | from anuga.shallow_water.data_manager import copy_code_files |
---|
16 | from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts |
---|
17 | from anuga.shallow_water.data_manager import sww2dem |
---|
18 | from os import sep |
---|
19 | |
---|
20 | # Application specific imports |
---|
21 | from setup_model import project # Definition of file names and polygons |
---|
22 | |
---|
23 | |
---|
24 | #------------------------------------------------------------------------------- |
---|
25 | # Setup |
---|
26 | #------------------------------------------------------------------------------- |
---|
27 | |
---|
28 | scenario_name = 'fitting_problem_' + project.extent |
---|
29 | combined_elevation_basename = join(project.combined_elevation, 'original') + '_' + scenario_name |
---|
30 | combined_elevation = combined_elevation_basename + '.pts' |
---|
31 | |
---|
32 | # Initial bounding polygon for data clipping |
---|
33 | bounding_polygon = read_polygon(join(project.polygons_folder, project.extent)+'.csv') |
---|
34 | |
---|
35 | output_run = join(project.output_folder, project.run_time) + scenario_name |
---|
36 | |
---|
37 | #------------------------------------------------------------------------------- |
---|
38 | # Copy scripts to time stamped output directory and capture screen |
---|
39 | # output to file. Copy script must be before screen_catcher |
---|
40 | #------------------------------------------------------------------------------- |
---|
41 | copy_code_files(output_run, __file__, |
---|
42 | os.path.join(os.path.dirname(project.__file__), |
---|
43 | project.__name__+'.py')) |
---|
44 | start_screen_catcher(output_run, 0, 1) |
---|
45 | |
---|
46 | #----------------------------------------------------------------------------- |
---|
47 | # Preparation of topographic data |
---|
48 | # |
---|
49 | # Convert ASC 2 DEM 2 PTS using source data and store result in source data |
---|
50 | # Do for coarse and fine data |
---|
51 | # Fine pts file to be clipped to area of interest |
---|
52 | #----------------------------------------------------------------------------- |
---|
53 | |
---|
54 | print 'bounding_polygon', bounding_polygon |
---|
55 | print 'combined_elevation_basename', combined_elevation_basename |
---|
56 | if not os.path.exists(combined_elevation): |
---|
57 | |
---|
58 | # Create Geospatial data from ASCII files |
---|
59 | geospatial_data = {} |
---|
60 | topographies_folder = join(project.topographies_folder, 'original') |
---|
61 | if not project.ascii_grid_filenames == []: |
---|
62 | for filename in project.ascii_grid_filenames: |
---|
63 | absolute_filename = join(topographies_folder, filename) |
---|
64 | convert_dem_from_ascii2netcdf(absolute_filename, |
---|
65 | basename_out=absolute_filename, |
---|
66 | use_cache=True, |
---|
67 | verbose=True) |
---|
68 | dem2pts(absolute_filename, use_cache=True, verbose=True) |
---|
69 | |
---|
70 | G_grid = Geospatial_data(file_name=absolute_filename+'.pts', |
---|
71 | verbose=True) |
---|
72 | print 'Clip geospatial object' |
---|
73 | geospatial_data[filename] = G_grid.clip(bounding_polygon) |
---|
74 | |
---|
75 | # Create Geospatial data from TXT files |
---|
76 | if not project.point_filenames == []: |
---|
77 | for filename in project.point_filenames: |
---|
78 | absolute_filename = join(topographies_folder, filename) |
---|
79 | G_points = Geospatial_data(file_name=absolute_filename, |
---|
80 | verbose=True) |
---|
81 | print 'Clip geospatial object' |
---|
82 | geospatial_data[filename] = G_points.clip(bounding_polygon) |
---|
83 | |
---|
84 | |
---|
85 | #------------------------------------------------------------------------------ |
---|
86 | # Combine, clip and export dataset |
---|
87 | #------------------------------------------------------------------------------ |
---|
88 | |
---|
89 | print 'Add geospatial objects' |
---|
90 | G = None |
---|
91 | for key in geospatial_data: |
---|
92 | G += geospatial_data[key] |
---|
93 | |
---|
94 | print 'Export combined DEM file' |
---|
95 | G.export_points_file(combined_elevation_basename + '.pts') |
---|
96 | |
---|
97 | print 'Do txt version too' |
---|
98 | # Use for comparision in ARC |
---|
99 | G.export_points_file(combined_elevation_basename + '.txt') |
---|
100 | |
---|
101 | |
---|
102 | #------------------------------------------------------------------------------ |
---|
103 | # Create ANUGA mesh and fit |
---|
104 | #------------------------------------------------------------------------------ |
---|
105 | |
---|
106 | tags = {'x': range(len(bounding_polygon))} # All segments tagged the same |
---|
107 | domain = create_domain_from_regions(bounding_polygon, |
---|
108 | boundary_tags=tags, |
---|
109 | maximum_triangle_area=project.extent_maxarea, |
---|
110 | interior_regions=project.interior_regions, |
---|
111 | mesh_filename=project.meshes[:-4] + '_' + scenario_name + '.msh', |
---|
112 | use_cache=True, |
---|
113 | verbose=True) |
---|
114 | print domain.statistics() |
---|
115 | |
---|
116 | |
---|
117 | domain.set_name(scenario_name) |
---|
118 | domain.set_datadir(output_run) |
---|
119 | domain.set_minimum_storable_height(0.01) # Don't store depth less than 1cm |
---|
120 | |
---|
121 | #------------------------------------------------------------------------------ |
---|
122 | # Setup initial conditions |
---|
123 | #------------------------------------------------------------------------------ |
---|
124 | |
---|
125 | print 'Setup initial conditions' |
---|
126 | domain.set_quantity('stage', 0.0) |
---|
127 | domain.set_quantity('friction', 0.01) |
---|
128 | domain.set_quantity('elevation', |
---|
129 | filename=combined_elevation_basename+'.pts', |
---|
130 | use_cache=True, |
---|
131 | verbose=True, |
---|
132 | alpha=0.1) |
---|
133 | |
---|
134 | #------------------------------------------------------------------------------ |
---|
135 | # Setup boundary conditions |
---|
136 | #------------------------------------------------------------------------------ |
---|
137 | |
---|
138 | print 'Set boundary - available tags:', domain.get_boundary_tags() |
---|
139 | Br = Reflective_boundary(domain) |
---|
140 | domain.set_boundary({'x': Br}) |
---|
141 | |
---|
142 | #----------------------------------------------------------------------------- |
---|
143 | # Take a few steps |
---|
144 | #----------------------------------------------------------------------------- |
---|
145 | |
---|
146 | |
---|
147 | for t in domain.evolve(yieldstep=1, |
---|
148 | finaltime=3): |
---|
149 | print domain.timestepping_statistics() |
---|
150 | print domain.boundary_statistics(tags='x') |
---|
151 | |
---|
152 | print 'Done' |
---|
153 | |
---|
154 | |
---|
155 | #----------------------------------------------------------------------------- |
---|
156 | # Exporting Grid |
---|
157 | # Define allowed variable names and associated equations to generate values. |
---|
158 | # This would not normally change. |
---|
159 | #----------------------------------------------------------------------------- |
---|
160 | var_equations = {'stage': 'stage', |
---|
161 | 'momentum': '(xmomentum**2 + ymomentum**2)**0.5', |
---|
162 | 'depth': 'stage-elevation', |
---|
163 | 'speed': '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)', |
---|
164 | 'elevation': 'elevation' } |
---|
165 | |
---|
166 | # one or more key strings from var_equations above |
---|
167 | var = ['elevation'] |
---|
168 | |
---|
169 | ###### |
---|
170 | # Start running the various conversions we require. |
---|
171 | ###### |
---|
172 | |
---|
173 | for which_var in var: |
---|
174 | if which_var not in var_equations: |
---|
175 | print 'Unrecognized variable name: %s' % which_var |
---|
176 | break |
---|
177 | |
---|
178 | name = join(output_run, scenario_name) |
---|
179 | outname = name + '_' + which_var |
---|
180 | quantityname = var_equations[which_var] |
---|
181 | |
---|
182 | print 'start sww2dem: name=%s' % name |
---|
183 | |
---|
184 | sww2dem(name, basename_out = outname, |
---|
185 | quantity = quantityname, |
---|
186 | timestep = 0, |
---|
187 | cellsize = 10, |
---|
188 | reduction = max, |
---|
189 | verbose = True, |
---|
190 | format = 'asc') |
---|
191 | |
---|
192 | |
---|
193 | |
---|