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

"""This scripts recreates fitting problem where values seem to go towards zero
in areas with missing data.
"""
from os.path import join
from anuga.utilities.polygon import read_polygon
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
from anuga.shallow_water.data_manager import dem2pts
from anuga.geospatial_data.geospatial_data import Geospatial_data
from anuga.interface import create_domain_from_regions
from anuga.interface import Reflective_boundary

import project, os
import sys
from anuga.shallow_water.data_manager import start_screen_catcher
from anuga.shallow_water.data_manager import copy_code_files
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
from anuga.shallow_water.data_manager import sww2dem
from os import sep

# Application specific imports
from setup_model import project   # Definition of file names and polygons


#-------------------------------------------------------------------------------
# Setup
#-------------------------------------------------------------------------------

scenario_name = 'fitting_problem_' + project.extent
combined_elevation_basename = join(project.combined_elevation, 'original') + '_' + scenario_name 
combined_elevation = combined_elevation_basename + '.pts'

# Initial bounding polygon for data clipping 
bounding_polygon = read_polygon(join(project.polygons_folder, project.extent)+'.csv')

output_run = join(project.output_folder, project.run_time) + scenario_name

#-------------------------------------------------------------------------------
# Copy scripts to time stamped output directory and capture screen
# output to file. Copy script must be before screen_catcher
#-------------------------------------------------------------------------------
copy_code_files(output_run, __file__, 
                os.path.join(os.path.dirname(project.__file__),
                             project.__name__+'.py'))
start_screen_catcher(output_run, 0, 1)

#-----------------------------------------------------------------------------
# Preparation of topographic data
#
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
# Do for coarse and fine data
# Fine pts file to be clipped to area of interest
#-----------------------------------------------------------------------------

print 'bounding_polygon', bounding_polygon
print 'combined_elevation_basename', combined_elevation_basename
if not os.path.exists(combined_elevation):
 
    # Create Geospatial data from ASCII files
    geospatial_data = {}
    topographies_folder = join(project.topographies_folder, 'original')
    if not project.ascii_grid_filenames == []:
        for filename in project.ascii_grid_filenames:
            absolute_filename = join(topographies_folder, filename)
            convert_dem_from_ascii2netcdf(absolute_filename,
                                          basename_out=absolute_filename,
                                          use_cache=True,
                                          verbose=True)
            dem2pts(absolute_filename, use_cache=True, verbose=True)

            G_grid = Geospatial_data(file_name=absolute_filename+'.pts',
                                                        verbose=True)
            print 'Clip geospatial object'
            geospatial_data[filename] = G_grid.clip(bounding_polygon)

    # Create Geospatial data from TXT files
    if not project.point_filenames == []:
        for filename in project.point_filenames:
            absolute_filename = join(topographies_folder, filename)
            G_points = Geospatial_data(file_name=absolute_filename,
                                                        verbose=True)
            print 'Clip geospatial object'
            geospatial_data[filename] = G_points.clip(bounding_polygon)


    #------------------------------------------------------------------------------
    # Combine, clip and export dataset 
    #------------------------------------------------------------------------------

    print 'Add geospatial objects'
    G = None
    for key in geospatial_data:
        G += geospatial_data[key]
                    
    print 'Export combined DEM file'
    G.export_points_file(combined_elevation_basename + '.pts')

    print 'Do txt version too'
    # Use for comparision in ARC
    G.export_points_file(combined_elevation_basename + '.txt')


#------------------------------------------------------------------------------
# Create ANUGA mesh and fit
#------------------------------------------------------------------------------

tags = {'x': range(len(bounding_polygon))} # All segments tagged the same
domain = create_domain_from_regions(bounding_polygon,
                                    boundary_tags=tags,
                                    maximum_triangle_area=project.extent_maxarea,
                                    interior_regions=project.interior_regions,
                                    mesh_filename=project.meshes[:-4] + '_' + scenario_name + '.msh',
                                    use_cache=True,
                                    verbose=True)
print domain.statistics()


domain.set_name(scenario_name)
domain.set_datadir(output_run) 
domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------

print 'Setup initial conditions'
domain.set_quantity('stage', 0.0)
domain.set_quantity('friction', 0.01)
domain.set_quantity('elevation', 
                    filename=combined_elevation_basename+'.pts',
                    use_cache=True,
                    verbose=True,
                    alpha=0.1)

#------------------------------------------------------------------------------
# Setup boundary conditions 
#------------------------------------------------------------------------------

print 'Set boundary - available tags:', domain.get_boundary_tags()
Br = Reflective_boundary(domain)
domain.set_boundary({'x': Br}) 

#-----------------------------------------------------------------------------
# Take a few steps
#-----------------------------------------------------------------------------


for t in domain.evolve(yieldstep=1, 
                       finaltime=3): 
    print domain.timestepping_statistics()
    print domain.boundary_statistics(tags='x')

print 'Done'


#-----------------------------------------------------------------------------
# Exporting Grid
# Define allowed variable names and associated equations to generate values.
# This would not normally change.
#-----------------------------------------------------------------------------
var_equations = {'stage':     'stage',
                 'momentum':  '(xmomentum**2 + ymomentum**2)**0.5',
                 'depth':     'stage-elevation',
                 'speed':     '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)',
                 'elevation': 'elevation' }

# one or more key strings from var_equations above
var = ['elevation']

######
# Start running the various conversions we require.
######

for which_var in var:
    if which_var not in var_equations:
        print 'Unrecognized variable name: %s' % which_var
        break
    
    name = join(output_run, scenario_name)
    outname = name + '_' + which_var
    quantityname = var_equations[which_var]

    print 'start sww2dem: name=%s' % name
    
    sww2dem(name, basename_out = outname,
                quantity = quantityname,
                timestep = 0,
                cellsize = 10,    
                reduction = max, 
                verbose = True,
                format = 'asc')