Changeset 7647
- Timestamp:
- Mar 2, 2010, 2:31:30 PM (15 years ago)
- Location:
- anuga_work/production
- Files:
-
- 1 added
- 1 deleted
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/production/bunbury_storm_surge_2009/build_elevation.py
r7613 r7647 84 84 keylist = project.ascii_grid_filenames + project.point_filenames 85 85 86 G1 = geospatial_data['swwa_10m '].clip_outside(project.elevation_clip_box)87 G2 = geospatial_data['bunbury_5m '].clip(project.elevation_clip_box)86 G1 = geospatial_data['swwa_10m_IC'].clip_outside(project.elevation_clip_box) 87 G2 = geospatial_data['bunbury_5m_IC'].clip(project.elevation_clip_box) 88 88 #~ G2 = geospatial_data['bunbury_nth_a'].clip(project.elevation_clip_box) 89 89 #~ G3 = geospatial_data['bunbury_nth_b'].clip(project.elevation_clip_box) … … 94 94 G5 = geospatial_data['Busselton_NavyFinal_Clip_ss.txt'] 95 95 G6 = geospatial_data['Leschenault_TIN.txt'] 96 G7 = geospatial_data['DPI5U1A02_01a_edited.txt'] 97 G8 = geospatial_data['DPI5U1A02_01b_edited.txt'] 98 G9 = geospatial_data['DPI5U1A02_01c_edited.txt'] 99 G10 = geospatial_data['DPI5U1A02_01d_edited.txt'] 100 G11 = geospatial_data['DPI5U1A02_01e_edited.txt'] 96 G7 = geospatial_data['Leschenault_Inlet_TIN.txt'] 97 G8 = geospatial_data['DPI5U1A02_01a_edited.txt'] 98 G9 = geospatial_data['DPI5U1A02_01b_edited.txt'] 99 G10 = geospatial_data['DPI5U1A02_01c_edited.txt'] 100 G11 = geospatial_data['DPI5U1A02_01d_edited.txt'] 101 G12 = geospatial_data['DPI5U1A02_01e_edited.txt'] 101 102 102 103 #~ G_list = [G1, G2, G3, G4, G5, G6, G7, G8, G9, G10] … … 113 114 #~ print 'G10', G10.attributes 114 115 115 G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 + G11 116 G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 + G9 + G10 + G11 + G12 116 117 117 118 # G = None … … 145 146 G10.export_points_file(join(project.topographies_folder, keylist[9] +'_export.txt')) 146 147 G11.export_points_file(join(project.topographies_folder, keylist[10] +'_export.txt')) 148 G12.export_points_file(join(project.topographies_folder, keylist[11] +'_export.txt')) -
anuga_work/production/bunbury_storm_surge_2009/get_timeseries.py
r7613 r7647 19 19 directory = project.output_folder 20 20 21 time_dir1 = '20090508_150215_run_final_0_51469_lfountai' 22 time_dir2 = '20090511_161526_run_final_1.1_51469_kvanputt' 23 time_dir3 = '20090511_165539_run_final_0_50863_lfountai' 24 time_dir4 = '20090521_220101_run_final_1.1_50863_kvanputt' 25 time_dir5 = '20090522_164526_run_final_0_51392_lfountai' 26 time_dir6 = '20090522_164640_run_final_1.1_51392_lfountai' 27 time_dir7 = '20090518_154710_run_final_0_50994_lfountai' 28 time_dir8 = '20090519_160510_run_final_1.1_50994_lfountai' 29 time_dir9 = '20090522_164948_run_final_0_51423_lfountai' 30 time_dir10 = '20090522_165600_run_final_1.1_51423_lfountai' 21 time_dir1 = '20100103_102505_run_storm_surge_final_0_alby_coarse_lfountai' 31 22 32 time_dirs = [time_dir1, time_dir2, time_dir3, time_dir4, time_dir5, 33 time_dir6, time_dir7, time_dir8, time_dir9, time_dir10]23 24 time_dirs = [time_dir1] 34 25 35 26 for time_dir in time_dirs: -
anuga_work/production/bunbury_storm_surge_2009/project.py
r7613 r7647 46 46 # storm_surge_final - as 'final' but with a longer yieldstep (12 mins) 47 47 48 #------------------------------------------------------------------------------49 #Floodgate parameters50 #------------------------------------------------------------------------------51 52 #changeable parameters specific to the floodgates53 54 floodgate_height = 4.4 #the relief of the floodgates above adjacent topography55 floodgate_thickness = 1 #the thickness of the floodgates in meters56 floodgate_LHS_x = 373546.02 #the x coordinate of the left hand side of the floodgates57 floodgate_LHS_y = 6312332.62 #the y coordinate of the left hand side of the floodgates58 floodgate_RHS_x = 373553.97 #the x coordinate of the right hand side of the floodgates59 floodgate_RHS_y = 6312331.97 #the y coordinate of the right hand side of the floodgates60 61 floodgate_start = 86400 #the time in seconds at which the floodgates begin to close62 floodgate_duration = 10 #the number of seconds that it takes the floodgates to close63 64 #polygon around the floodgates used to define an interior region when building mesh65 floodgate_boundary = [[floodgate_LHS_x,(floodgate_LHS_y - floodgate_thickness)],66 [floodgate_RHS_x,(floodgate_RHS_y - floodgate_thickness)],67 [floodgate_RHS_x,(floodgate_RHS_y + floodgate_thickness)],68 [floodgate_LHS_x,(floodgate_LHS_y + floodgate_thickness)]]69 floodgate_resolution = 1.0 #max triangle area to use when building mesh around floodgates70 48 71 49 #------------------------------------------------------------------------------- … … 91 69 # Used in build_elevation.py 92 70 # Format for ascii grids, as produced in ArcGIS + a projection file 93 ascii_grid_filenames = ['swwa_10m', # LiDAR mosaiced and resampled to 10m 94 'bunbury_5m'] # Bunbury LiDar grid resampled to 5 m - 1m data caused memory errors 95 #~ 'bunbury_nth_a', # 1m LiDAR from WA Dot - nothern quarter (file split as too large for dem2pts) 96 #~ 'bunbury_nth_b', # 1m LiDAR from WA Dot - 2nd northmost quarter (file split as too large for dem2pts) 97 #~ 'bunbury_sth_a', # 1m LiDAR from WA Dot - 2nd southmost quarter (file split as too large for dem2pts) 98 #~ 'bunbury_sth_b'] # 1m LiDAR from WA Dot - southern quarter (file split as too large for dem2pts) 99 # 'bunbury_aoi'] # 1m LiDAR from WA DoT - clipped to area_of_interest.csv 71 ascii_grid_filenames = ['swwa_10m_IC', # LiDAR mosaiced and resampled to 10m, clipped to IC to remove values in water 72 'bunbury_5m_IC'] # Bunbury LiDar grid resampled to 5 m (1m data caused memory errors) 73 # and clipped to Initial Conditions (to remove values in water) 74 # 'bunbury_nth_a', # 1m LiDAR from WA Dot - nothern quarter (file split as too large for dem2pts) 75 # 'bunbury_nth_b', # 1m LiDAR from WA Dot - 2nd northmost quarter (file split as too large for dem2pts) 76 # 'bunbury_sth_a', # 1m LiDAR from WA Dot - 2nd southmost quarter (file split as too large for dem2pts) 77 # 'bunbury_sth_b'] # 1m LiDAR from WA Dot - southern quarter (file split as too large for dem2pts) 78 # 'bunbury_aoi'] # 1m LiDAR from WA DoT - clipped to area_of_interest.csv 100 79 # 'bunbury_nth', # 1m LiDAR from WA Dot - nothern half (file split as too large for dem2pts) 101 80 # 'bunbury_sth'] # 1m LiDAR data from WA DoT - southern half … … 103 82 104 83 # Format for point is x,y,elevation (with header) 105 point_filenames = ['DPI.txt', # Bathymetry data from DPI 106 'Busselton_Chart_Clip_ss.txt', # Clipped from Busselton_Chart - see Busselton Tsunami Scenario 2009 107 'Busselton_NavyFinal_Clip_ss.txt', # Clipped from Busselton_NavyFinal - see Busselton Tsunami Scenario 2009 108 'Leschenault_TIN.txt', # TIN created over the Leschenault Estuary 109 'DPI5U1A02_01a_edited.txt', 'DPI5U1A02_01b_edited.txt', 110 'DPI5U1A02_01c_edited.txt','DPI5U1A02_01d_edited.txt', 111 'DPI5U1A02_01e_edited.txt'] # Bathymetric LiDAR from DPI - split into manageable pieces 84 point_filenames = ['DPI.txt', # Bathymetry data from DPI 85 'Busselton_Chart_Clip_ss.txt', # Clipped from Busselton_Chart - see Busselton Tsunami Scenario 2009 86 'Busselton_NavyFinal_Clip_ss.txt', # Clipped from Busselton_NavyFinal - see Busselton Tsunami Scenario 2009 87 'Leschenault_TIN.txt', # TIN created over the Leschenault Estuary 88 'Leschenault_inlet_TIN.txt' # TIN created over the Leschenault Inlet 89 'DPI5U1A02_01a_edited.txt', # Bathymetric LiDAR from DPI - split into manageable pieces 90 'DPI5U1A02_01b_edited.txt', 91 'DPI5U1A02_01c_edited.txt', 92 'DPI5U1A02_01d_edited.txt', 93 'DPI5U1A02_01e_edited.txt'] 112 94 113 95 # BOUNDING POLYGON - for data clipping and estimate of triangles in mesh … … 122 104 interior_regions_data = [['intermediate.csv', 2500], 123 105 ['area_of_interest.csv', 100], 124 [floodgate_boundary, floodgate_resolution]] 106 ['storm_gate_area.csv', 1], 107 ['floodgates.csv', floodgate_resolution]] 125 108 126 109 # LAND - used to set the initial stage/water to be offcoast only -
anuga_work/production/bunbury_storm_surge_2009/run_model.py
r7613 r7647 159 159 'ocean': Bf}) 160 160 161 #-----------------------------------------------------------------------------162 # Setup dynamic topography163 #-----------------------------------------------------------------------------164 from project import floodgate_height, floodgate_thickness,floodgate_LHS_x,floodgate_LHS_y,floodgate_RHS_x,floodgate_RHS_y,floodgate_start,floodgate_duration165 166 floodgate_end = floodgate_start + floodgate_duration #the time in seconds at which the floodgates are fully closed167 168 def floodgate(x,y):169 z = 0*x170 N = len(x)171 for i in range(N):172 ymin = ((floodgate_RHS_y - floodgate_LHS_y)/(floodgate_RHS_x - floodgate_LHS_x))*(x[i]-floodgate_LHS_x)+ floodgate_LHS_y - 0.5*floodgate_thickness173 ymax = ((floodgate_RHS_y - floodgate_LHS_y)/(floodgate_RHS_x - floodgate_LHS_x))*(x[i]-floodgate_LHS_x)+ floodgate_LHS_y + 0.5*floodgate_thickness174 xmin = floodgate_LHS_x175 xmax = floodgate_RHS_x176 if ymin < y[i] < ymax and xmin < x[i] < xmax:177 z[i] += (floodgate_height/((floodgate_duration)/project.yieldstep))178 179 return z180 181 161 #------------------------------------------------------------------------------- 182 162 # Evolve system through time 183 163 #------------------------------------------------------------------------------- 184 185 raising_floodgates = False186 lowering_floodgates = False187 164 188 165 t0 = time.time() … … 193 170 print domain.boundary_statistics(tags='ocean') 194 171 195 #raise floodgate at a given time interval (to make floodgate lower again repeat this code below196 #for a later time interval and change rising to false and falling to true197 if floodgate_start < t < floodgate_end:198 rising = True199 falling = False200 201 if rising:202 domain.add_quantity('elevation',floodgate)203 204 if falling:205 domain.add_quantity('elevation', lambda x,y: -1*floodgate(x,y))206 172 207 173 print 'Simulation took %.2f seconds' % (time.time()-t0) -
anuga_work/production/mandurah_storm_surge_2009/build_elevation.py
r7627 r7647 87 87 G1 = geospatial_data['m_peel_aoi'].clip(project.elevation_clip_box) 88 88 G2 = geospatial_data['ph10m_ss'].clip_outside(project.elevation_clip_box) 89 G3 = geospatial_data['MA-46893-SNDS_AHD.csv'] 90 G4 = geospatial_data['MS0205HY_AHD.csv'] 91 G5 = geospatial_data['MS0404_AHD.csv'] 92 G6 = geospatial_data['YU0403HY_AHD.csv'] 93 G7 = geospatial_data['original_data_ss.csv'] 89 G3 = geospatial_data['DpiU1A03_ss.csv'] 90 G4 = geospatial_data['MA-46893-SNDS_AHD.csv'] 91 G5 = geospatial_data['MS0205HY_AHD.csv'] 92 G6 = geospatial_data['MS0404_AHD.csv'] 93 G7 = geospatial_data['YU0403HY_AHD.csv'] 94 G8 = geospatial_data['original_data_ss.csv'] 94 95 95 G = G1 + G2 + G3 + G4 + G5 + G6 + G7 96 G = G1 + G2 + G3 + G4 + G5 + G6 + G7 + G8 96 97 97 98 print 'Export combined DEM file' -
anuga_work/production/mandurah_storm_surge_2009/gems_comparison/build_elevation_GEMS.py
r7628 r7647 82 82 print 'project_GEMS.elevation_clip_box', project_GEMS.elevation_clip_box 83 83 84 G1 = geospatial_data[' Man_25m'].clip(project_GEMS.elevation_clip_box)84 G1 = geospatial_data['man_25m_prj'].clip(project_GEMS.elevation_clip_box) 85 85 G2 = geospatial_data['man10m_ss'].clip_outside(project_GEMS.elevation_clip_box) 86 86 -
anuga_work/production/mandurah_storm_surge_2009/gems_comparison/project_GEMS.py
r7628 r7647 62 62 #------------------------------------------------------------------------------- 63 63 64 output_comment = [setup, tide, event, ' _GEMS_compare'] # event_number will have to64 output_comment = [setup, tide, event, 'GEMS_compare'] # event_number will have to 65 65 # change to something relevent 66 66 # for storm surge … … 77 77 ## 'm_peel_10m.asc' 78 78 ## 'm_harvey_10m.asc'] 79 ascii_grid_filenames = [' Man_25m', # this is the latest 25m DEM from GEMS79 ascii_grid_filenames = ['man_25m_prj', # this is the latest 25m DEM from GEMS 80 80 'man10m_ss'] # this is to fill in areas not covered by the 81 81 # GEMS grid and is derived from the 10m resampled LiDAR from DoT, -
anuga_work/production/mandurah_storm_surge_2009/gems_comparison/run_model_GEMS.py
r7628 r7647 110 110 print domain.statistics() 111 111 112 domain.s tarttime(project_GEMS.starttime)112 domain.set_starttime(project_GEMS.starttime) 113 113 domain.set_name(project_GEMS.scenario_name) 114 114 domain.set_datadir(project_GEMS.output_run) -
anuga_work/production/mandurah_storm_surge_2009/get_timeseries.py
r7614 r7647 19 19 directory = project.output_folder 20 20 21 time_dir1 = '20090508_150215_run_final_0_51469_lfountai' 22 time_dir2 = '20090511_161526_run_final_1.1_51469_kvanputt' 23 time_dir3 = '20090511_165539_run_final_0_50863_lfountai' 24 time_dir4 = '20090521_220101_run_final_1.1_50863_kvanputt' 25 time_dir5 = '20090522_164526_run_final_0_51392_lfountai' 26 time_dir6 = '20090522_164640_run_final_1.1_51392_lfountai' 27 time_dir7 = '20090518_154710_run_final_0_50994_lfountai' 28 time_dir8 = '20090519_160510_run_final_1.1_50994_lfountai' 29 time_dir9 = '20090522_164948_run_final_0_51423_lfountai' 30 time_dir10 = '20090522_165600_run_final_1.1_51423_lfountai' 21 time_dir1 = '20091229_100831_run_storm_surge_final_0_alby_waves_lfountai' 31 22 32 time_dirs = [time_dir1, time_dir2, time_dir3, time_dir4, time_dir5, 33 time_dir6, time_dir7, time_dir8, time_dir9, time_dir10] 23 time_dirs = [time_dir1] 34 24 35 25 for time_dir in time_dirs: -
anuga_work/production/mandurah_storm_surge_2009/project.py
r7627 r7647 86 86 # Format for point is x,y,elevation (with header) 87 87 # Don't use these for comparison with GEMS inundation 88 point_filenames = ['MA-46893-SNDS_AHD.csv', # These files contain inf-fill bathymetry for 88 point_filenames = ['DpiU1A03_ss.csv', # Bathymetric LiDAR data from DPI, clipped to bounding polygon 89 'MA-46893-SNDS_AHD.csv', # These files contain inf-fill bathymetry for 89 90 'MS0205HY_AHD.csv', # the Peel & Harvey estuaries, as well as around 90 91 'MS0404_AHD.csv', # the canals -
anuga_work/production/mandurah_storm_surge_2009/run_model.py
r7614 r7647 110 110 print domain.statistics() 111 111 112 domain.set_starttime(project.starttime) 112 113 domain.set_name(project.scenario_name) 113 114 domain.set_datadir(project.output_run)
Note: See TracChangeset
for help on using the changeset viewer.