from anuga.geospatial_data.geospatial_data import * from os import sep, environ, getenv, getcwd,umask 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 anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files from math import sqrt from time import localtime, strftime, gmtime umask(002) time = strftime('%Y%m%d_%H%M%S',gmtime()) home = getenv('INUNDATIONHOME') + sep +'data'+sep output_dir = home + 'anuga_validation'+sep+'cocos'+sep+'anuga'+sep+'outputs'+sep+time+sep copy_code_files(output_dir,__file__) start_screen_catcher(output_dir,verbose=True) print'testing' 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 grid1_res=370.4 # #grid2_north_max=8713178.04 #grid2_north_min=8625096.92 #grid2_east_max=294652.80 #grid2_east_min=224424.96 grid2_res=74.08 # #grid3_north_max=274295.91 #grid3_north_min=266767.44 #grid3_east_max=8665112.56 #grid3_east_min=8657346.91 grid3_res=14.82 # #grid4_north_max=270900.24 #grid4_north_min=269701.45 #grid4_east_max=8659120.05 #grid4_east_min=8659120.05 grid4_res=2.96 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 = 100 #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,(grid2_res**2/2)*res],[poly3,(grid3_res**2/2)*res],[poly4,(grid4_res**2/2)*res]] interior_polys = [[poly2,(grid2_res**2)*2],[poly3,(grid3_res**2)*res],[poly4,(grid4_res**2)*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) max_tri = (grid1_res**2/2)*res print 'max_tri',max_tri create_mesh_from_regions(poly1, #create_mesh_from_regions(poly3, boundary_tags={'side': [0,1,2,3]}, maximum_triangle_area=(grid1_res**2)*res, # maximum_triangle_area=(grid1_res**2)*1, interior_regions=interior_polys, filename=mesh_file, use_cache=False, 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.1], # boundary_poly=poly3, # mesh_file=mesh_file, # plot_name=None, # split_factor=0.1, # 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) alpha = 0.01) print 'Finished Set quantity' domain.set_name(sww_file) domain.set_datadir(output_dir) domain.set_default_order(1) # 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) 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(output_dir+sww_file+".sww", extra_name_out = 'grid4', quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation timestep = None, reduction = max, # cellsize = grid4_res*sqrt(1), cellsize = grid4_res*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(output_dir+sww_file+".sww", extra_name_out = 'grid3', quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation timestep = None, reduction = max, # cellsize = grid3_res*sqrt(1), cellsize = grid3_res*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(output_dir+sww_file+".sww", extra_name_out = 'grid2', quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation timestep = None, reduction = max, cellsize = grid2_res*sqrt(1), # cellsize = grid2_res*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(output_dir+sww_file+".sww", extra_name_out = 'grid1', quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation timestep = None, reduction = max, # cellsize = grid1_res*sqrt(1), cellsize = grid1_res*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')