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

Last change on this file since 4926 was 4926, checked in by nick, 16 years ago

updates to cocos grid development code

File size: 8.5 KB
Line 
1
2
3from anuga.geospatial_data.geospatial_data import *
4from os import sep, environ, getenv, getcwd,umask
5from anuga.pmesh.mesh_interface import create_mesh_from_regions
6from anuga.shallow_water import Domain
7from anuga.shallow_water import Reflective_boundary
8from anuga.shallow_water.data_manager import export_grid
9from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
10from math import sqrt
11from time import localtime, strftime, gmtime
12
13umask(002)
14time = strftime('%Y%m%d_%H%M%S',gmtime())
15
16home = getenv('INUNDATIONHOME') + sep +'data'+sep
17output_dir = home + 'anuga_validation'+sep+'cocos'+sep+'anuga'+sep+'outputs'+sep+time+sep
18
19copy_code_files(output_dir,__file__)
20
21start_screen_catcher(output_dir,verbose=True)
22print'testing'
23
24data_dir='/d/ehn/4/ehn/tsunami/data/bathymetry/Cocos/bathymetry/'
25data_dir_name = data_dir + 'new_edited_cocos_final'
26
27data_dir_name_small=data_dir_name+'_small'
28data_dir_name_smaller=data_dir_name+'_smaller'
29data_dir_name_small_clip=data_dir_name+'_small_clip'
30data_dir_name_clip=data_dir_name+'_clip'
31#print'create Geospatial data2 objects from coast'
32#G = Geospatial_data(file_name = data_dir_name + '.txt')
33#print'start split'
34#G_small, G_other = G.split(factor=0.1, verbose=True)
35#print'start export'
36#G_small.export_points_file(data_dir_name_smaller + '.txt')
37
38mesh_file=data_dir+'cocos.msh'
39sww_file='test_final'
40
41#grid1_north_max=8783110.00
42#grid1_north_min=8565314.80
43#grid1_east_max=439850.00
44#grid1_east_min=156494.00
45grid1_res=370.4
46#
47#grid2_north_max=8713178.04
48#grid2_north_min=8625096.92
49#grid2_east_max=294652.80
50#grid2_east_min=224424.96
51grid2_res=74.08
52#
53#grid3_north_max=274295.91
54#grid3_north_min=266767.44
55#grid3_east_max=8665112.56
56#grid3_east_min=8657346.91
57grid3_res=14.82
58#
59#grid4_north_max=270900.24
60#grid4_north_min=269701.45
61#grid4_east_max=8659120.05
62#grid4_east_min=8659120.05
63grid4_res=2.96
64
65poly1 = [[439850.00,8565314.80],[439850.00, 8783110.00],
66        [156494.00,8783110.00],[156494.00,8565314.80]]
67#        68450m2
68
69poly2 = [[294652.80,8625096.92],[294652.80,8713178.04],
70         [224424.96,8713178.04],[224424.96,8625096.92]]
71#        2738m2
72
73poly3 = [[274295.91,8657346.91],[274295.91,8665112.56],
74         [266767.44,8665112.56],[266767.44,8657346.91]]
75#        110m2
76       
77poly4 = [[270900.24,8659120.05],[270900.24,8660170.84],
78         [269701.45,8660170.84],[269701.45,8659120.05]]
79#        4.3m2
80
81res = 100
82#G = Geospatial_data(file_name = data_dir_name + '.txt')
83#G1=G.clip(poly1)
84#G1.export_points_file(data_dir_name_clip + '.txt')
85
86
87#interior_polys = [[poly2,(grid2_res**2/2)*res],[poly3,(grid3_res**2/2)*res],[poly4,(grid4_res**2/2)*res]]
88interior_polys = [[poly2,(grid2_res**2)*2],[poly3,(grid3_res**2)*res],[poly4,(grid4_res**2)*res]]
89#interior_polys = [[poly4,4.3*res]]
90
91#create_mesh_from_regions(poly1,
92#                         boundary_tags={'side': [0,1,2,3]},
93#                         maximum_triangle_area=68500*res,
94#                         interior_regions=interior_polys,
95#                         filename=mesh_file,
96#                         use_cache=True,
97#                         verbose=True)
98max_tri = (grid1_res**2/2)*res
99print 'max_tri',max_tri
100create_mesh_from_regions(poly1,
101#create_mesh_from_regions(poly3,
102                         boundary_tags={'side': [0,1,2,3]},
103                         maximum_triangle_area=(grid1_res**2)*res,
104#                         maximum_triangle_area=(grid1_res**2)*1,
105#                         maximum_triangle_area=110*res,
106                         interior_regions=interior_polys,
107                         filename=mesh_file,
108                         use_cache=True,
109                         verbose=True)
110                         
111                         
112#value, alpha = find_optimal_smoothing_parameter(data_file=data_dir_name_clip + '.txt',
113##                                             alpha_list=[0.0001,0.001,0.01,0.1,1.0],
114#                                             alpha_list=[0.1],
115#                                             boundary_poly=poly3,
116#                                             mesh_file=mesh_file,
117#                                             plot_name=None,
118#                                             split_factor=0.1,
119#                                             cache=True,
120#                                             verbose=True)
121                                             
122#print 'Setup computational domain', value, alpha
123
124    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
125    #above don't work
126domain = Domain(mesh_file, use_cache=True, verbose=True)
127
128print 'Start Set quantity'
129
130domain.set_quantity('elevation', 
131                            filename = data_dir_name_clip + '.txt',
132                            use_cache = True,
133                            verbose = True,
134#                            alpha = alpha)
135                            alpha = 0.01)
136print 'Finished Set quantity'
137
138domain.set_name(sww_file)
139domain.set_datadir(output_dir)
140domain.set_default_order(1) # Apply second order scheme
141domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
142domain.set_store_vertices_uniquely(False)
143domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
144#    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
145domain.beta_h = 0
146Br = Reflective_boundary(domain)
147
148
149domain.set_boundary({'side': Br})
150
151for t in domain.evolve(yieldstep = 1,finaltime =1): 
152        domain.write_time()
153        domain.write_boundary_statistics(tags = 'side') 
154
155export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid4',
156                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
157                    timestep = None,
158                    reduction = max,
159#                    cellsize = grid4_res*sqrt(1),
160                    cellsize = grid4_res*sqrt(res),
161                    NODATA_value = -9999,
162                    easting_min = 269701.45,
163                    easting_max = 270900.24,
164                    northing_min = 8659120.05,
165                    northing_max = 8660170.84,
166                    verbose = True,
167                    origin = None,
168                    datum = 'WGS84',
169                    format = 'asc')
170                         
171export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid3',
172                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
173                    timestep = None,
174                    reduction = max,
175#                    cellsize = grid3_res*sqrt(1),
176                    cellsize = grid3_res*sqrt(res),
177                    NODATA_value = -9999,
178                    easting_min = 266767.44,
179                    easting_max = 274295.91,
180                    northing_min = 8657346.91,
181                    northing_max = 8665112.56,
182                    verbose = True,
183                    origin = None,
184                    datum = 'WGS84',
185                    format = 'asc')
186
187export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid2',
188                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
189                    timestep = None,
190                    reduction = max,
191                    cellsize = grid2_res*sqrt(1),
192#                    cellsize = grid2_res*sqrt(res),
193                    NODATA_value = -9999,
194                    easting_min = 224424.96,
195                    easting_max = 294652.80,
196                    northing_min = 8625096.92,
197                    northing_max = 8713178.04,
198                    verbose = True,
199                    origin = None,
200                    datum = 'WGS84',
201                    format = 'asc')
202
203export_grid(output_dir+sww_file+".sww", extra_name_out = 'grid1',
204                    quantities = ['elevation'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
205                    timestep = None,
206                    reduction = max,
207#                    cellsize = grid1_res*sqrt(1),
208                    cellsize = grid1_res*sqrt(res),
209                    NODATA_value = -9999,
210                    easting_min = 156494.00,
211                    easting_max = 439850.00,
212                    northing_min = 8565314.80,
213                    northing_max = 8783110.00,
214                    verbose = True,
215                    origin = None,
216                    datum = 'WGS84',
217                    format = 'asc')
218
219
220
221
Note: See TracBrowser for help on using the repository browser.