source: anuga_work/development/cocos/grid_data.py @ 4832

from anuga.geospatial_data.geospatial_data import *
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water.data_manager import export_grid
from math import sqrt

data_dir='/d/ehn/4/ehn/tsunami/data/bathymetry/Cocos/bathymetry/'
data_dir_name = data_dir + 'new_edited_cocos_final'

data_dir_name_small=data_dir_name+'_small'
data_dir_name_smaller=data_dir_name+'_smaller'
data_dir_name_small_clip=data_dir_name+'_small_clip'
data_dir_name_clip=data_dir_name+'_clip'
#print'create Geospatial data2 objects from coast'
#G = Geospatial_data(file_name = data_dir_name + '.txt')
#print'start split'
#G_small, G_other = G.split(factor=0.1, verbose=True)
#print'start export'
#G_small.export_points_file(data_dir_name_smaller + '.txt')

mesh_file=data_dir+'cocos.msh'
sww_file='test_final'

#grid1_north_max=8783110.00
#grid1_north_min=8565314.80
#grid1_east_max=439850.00
#grid1_east_min=156494.00
#
#grid2_north_max=8713178.04
#grid2_north_min=8625096.92
#grid2_east_max=294652.80
#grid2_east_min=224424.96
#
#grid3_north_max=274295.91
#grid3_north_min=266767.44
#grid3_east_max=8665112.56
#grid3_east_min=8657346.91
#
#grid4_north_max=270900.24
#grid4_north_min=269701.45
#grid4_east_max=8659120.05
#grid4_east_min=8659120.05

poly1 = [[439850.00,8565314.80],[439850.00, 8783110.00],
        [156494.00,8783110.00],[156494.00,8565314.80]]
#        68450m2

poly2 = [[294652.80,8625096.92],[294652.80,8713178.04],
         [224424.96,8713178.04],[224424.96,8625096.92]]
#        2738m2

poly3 = [[274295.91,8657346.91],[274295.91,8665112.56],
         [266767.44,8665112.56],[266767.44,8657346.91]]
#        110m2
        
poly4 = [[270900.24,8659120.05],[270900.24,8660170.84],
         [269701.45,8660170.84],[269701.45,8659120.05]]
#        4.3m2

res = 1
#G = Geospatial_data(file_name = data_dir_name + '.txt')
#G1=G.clip(poly1)
#G1.export_points_file(data_dir_name_clip + '.txt')


interior_polys = [[poly2,2740*res],[poly3,110*res],[poly4,4.3*res]]
#interior_polys = [[poly4,4.3*res]]

#create_mesh_from_regions(poly1,
#                         boundary_tags={'side': [0,1,2,3]},
#                         maximum_triangle_area=68500*res,
#                         interior_regions=interior_polys,
#                         filename=mesh_file,
#                         use_cache=True,
#                         verbose=True)
create_mesh_from_regions(poly1,
#create_mesh_from_regions(poly3,
                         boundary_tags={'side': [0,1,2,3]},
                         maximum_triangle_area=68500*res,
#                         maximum_triangle_area=110*res,
                         interior_regions=interior_polys,
                         filename=mesh_file,
                         use_cache=True,
                         verbose=True)
                         
                         
value, alpha = find_optimal_smoothing_parameter(data_file=data_dir_name_clip + '.txt', 
                                             alpha_list=[0.0001,0.001,0.01,0.1,1.0],
#                                             alpha_list=[0.01],
                                             boundary_poly=poly3,
                                             mesh_file=mesh_file,
                                             plot_name=None,
                                             cache=True,
                                             verbose=True)
                                             
print 'Setup computational domain', value, alpha

    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
    #above don't work
domain = Domain(mesh_file, use_cache=True, verbose=True)

print 'Start Set quantity'

domain.set_quantity('elevation', 
                            filename = data_dir_name_clip + '.txt',
                            use_cache = True,
                            verbose = True,
                            alpha = alpha)
print 'Finished Set quantity'

domain.set_name(sww_file)
domain.set_datadir(data_dir)
domain.set_default_order(2) # Apply second order scheme
domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
domain.set_store_vertices_uniquely(False)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
#    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
domain.beta_h = 0
Br = Reflective_boundary(domain)


domain.set_boundary({'side': Br})

for t in domain.evolve(yieldstep = 1,finaltime =1): 
        domain.write_time()
        domain.write_boundary_statistics(tags = 'side')  

export_grid(data_dir+sww_file+".sww", extra_name_out = 'grid4',
                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
                    timestep = None,
                    reduction = max,
                    cellsize = 2.96*sqrt(res),
                    NODATA_value = -9999,
                    easting_min = 269701.45,
                    easting_max = 270900.24,
                    northing_min = 8659120.05,
                    northing_max = 8660170.84,
                    verbose = True,
                    origin = None,
                    datum = 'WGS84',
                    format = 'asc')
                         
export_grid(data_dir+sww_file+".sww", extra_name_out = 'grid3',
                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
                    timestep = None,
                    reduction = max,
                    cellsize = 14.82*sqrt(res),
                    NODATA_value = -9999,
                    easting_min = 266767.44,
                    easting_max = 274295.91,
                    northing_min = 8657346.91,
                    northing_max = 8665112.56,
                    verbose = True,
                    origin = None,
                    datum = 'WGS84',
                    format = 'asc')

export_grid(data_dir+sww_file+".sww", extra_name_out = 'grid2',
                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
                    timestep = None,
                    reduction = max,
                    cellsize = 74.08*sqrt(res),
                    NODATA_value = -9999,
                    easting_min = 224424.96,
                    easting_max = 294652.80,
                    northing_min = 8625096.92,
                    northing_max = 8713178.04,
                    verbose = True,
                    origin = None,
                    datum = 'WGS84',
                    format = 'asc')

export_grid(data_dir+sww_file+".sww", extra_name_out = 'grid1',
                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
                    timestep = None,
                    reduction = max,
                    cellsize = 370.4*sqrt(res),
                    NODATA_value = -9999,
                    easting_min = 156494.00,
                    easting_max = 439850.00,
                    northing_min = 8565314.80,
                    northing_max = 8783110.00,
                    verbose = True,
                    origin = None,
                    datum = 'WGS84',
                    format = 'asc')