source: anuga_work/production/hobart_2009/fitting_problem.py @ 7855

Last change on this file since 7855 was 7082, checked in by kristy, 16 years ago
File size: 7.6 KB
Line 
1"""This scripts recreates fitting problem where values seem to go towards zero
2in areas with missing data.
3"""
4from os.path import join
5from anuga.utilities.polygon import read_polygon
6from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
7from anuga.shallow_water.data_manager import dem2pts
8from anuga.geospatial_data.geospatial_data import Geospatial_data
9from anuga.interface import create_domain_from_regions
10from anuga.interface import Reflective_boundary
11
12import project, os
13import sys
14from anuga.shallow_water.data_manager import start_screen_catcher
15from anuga.shallow_water.data_manager import copy_code_files
16from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
17from anuga.shallow_water.data_manager import sww2dem
18from os import sep
19
20# Application specific imports
21from setup_model import project   # Definition of file names and polygons
22
23
24#-------------------------------------------------------------------------------
25# Setup
26#-------------------------------------------------------------------------------
27
28scenario_name = 'fitting_problem_' + project.extent
29combined_elevation_basename = join(project.combined_elevation, 'original') + '_' + scenario_name
30combined_elevation = combined_elevation_basename + '.pts'
31
32# Initial bounding polygon for data clipping
33bounding_polygon = read_polygon(join(project.polygons_folder, project.extent)+'.csv')
34
35output_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#-------------------------------------------------------------------------------
41copy_code_files(output_run, __file__, 
42                os.path.join(os.path.dirname(project.__file__),
43                             project.__name__+'.py'))
44start_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
54print 'bounding_polygon', bounding_polygon
55print 'combined_elevation_basename', combined_elevation_basename
56if 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
106tags = {'x': range(len(bounding_polygon))} # All segments tagged the same
107domain = 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)
114print domain.statistics()
115
116
117domain.set_name(scenario_name)
118domain.set_datadir(output_run) 
119domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
120
121#------------------------------------------------------------------------------
122# Setup initial conditions
123#------------------------------------------------------------------------------
124
125print 'Setup initial conditions'
126domain.set_quantity('stage', 0.0)
127domain.set_quantity('friction', 0.01)
128domain.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
138print 'Set boundary - available tags:', domain.get_boundary_tags()
139Br = Reflective_boundary(domain)
140domain.set_boundary({'x': Br}) 
141
142#-----------------------------------------------------------------------------
143# Take a few steps
144#-----------------------------------------------------------------------------
145
146
147for t in domain.evolve(yieldstep=1, 
148                       finaltime=3): 
149    print domain.timestepping_statistics()
150    print domain.boundary_statistics(tags='x')
151
152print '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#-----------------------------------------------------------------------------
160var_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
167var = ['elevation']
168
169######
170# Start running the various conversions we require.
171######
172
173for 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
Note: See TracBrowser for help on using the repository browser.