[3627] | 1 | """Common filenames and locations for topographic data, meshes and outputs. |
---|
| 2 | """ |
---|
| 3 | |
---|
[4147] | 4 | from os import sep, environ, getenv, getcwd |
---|
| 5 | from os.path import expanduser |
---|
[3802] | 6 | import sys |
---|
[4147] | 7 | from time import localtime, strftime, gmtime |
---|
| 8 | from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon, number_mesh_triangles |
---|
| 9 | from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm |
---|
[3937] | 10 | from anuga.utilities.system_tools import get_user_name |
---|
[3627] | 11 | |
---|
[3828] | 12 | # file and system info |
---|
| 13 | #--------------------------------- |
---|
[3627] | 14 | codename = 'project.py' |
---|
| 15 | |
---|
| 16 | home = getenv('INUNDATIONHOME') #Sandpit's parent dir |
---|
[3937] | 17 | user = get_user_name() |
---|
[3627] | 18 | |
---|
[3669] | 19 | # INUNDATIONHOME is the inundation directory, not the data directory. |
---|
| 20 | home += sep +'data' |
---|
[3627] | 21 | |
---|
[4147] | 22 | #time stuff |
---|
| 23 | time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir |
---|
| 24 | gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir |
---|
| 25 | build_time = time+'_build' |
---|
| 26 | run_time = time+'_run' |
---|
| 27 | print 'gtime: ', gtime |
---|
[3828] | 28 | |
---|
[4147] | 29 | tide = 0.6 |
---|
[3871] | 30 | |
---|
[4147] | 31 | #Making assumptions about the location of scenario data |
---|
| 32 | state = 'western_australia' |
---|
| 33 | scenario_name = 'dampier' |
---|
| 34 | scenario = 'dampier_tsunami_scenario_2006' |
---|
[3828] | 35 | |
---|
[4147] | 36 | # onshore data provided by WA DLI |
---|
| 37 | onshore_name = 'dampier_dli_ext' # original |
---|
| 38 | #island |
---|
| 39 | island_name = 'rott_dli_ext' # original |
---|
| 40 | island_name1 = 'gard_dli_ext' |
---|
| 41 | island_name2 = 'carnac_island_dted' |
---|
| 42 | island_name3 = 'penguin_dted' |
---|
[3940] | 43 | |
---|
[4147] | 44 | # AHO + DPI data + colin French coastline |
---|
| 45 | coast_name = 'waterline' |
---|
| 46 | offshore_name = 'dampier_bathymetry' |
---|
| 47 | offshore1_name = 'missing_fairsheets' |
---|
| 48 | |
---|
| 49 | #final topo name |
---|
[3828] | 50 | combined_name ='dampier_combined_elevation' |
---|
[4147] | 51 | combined_smaller_name = 'dampier_combined_elevation_smaller' |
---|
[3828] | 52 | |
---|
[4147] | 53 | topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep |
---|
| 54 | topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep |
---|
[3851] | 55 | topographies_time_dir = topographies_dir+build_time+sep |
---|
[3627] | 56 | |
---|
[4147] | 57 | # input topo file location |
---|
| 58 | onshore_in_dir_name = topographies_in_dir + onshore_name |
---|
| 59 | island_in_dir_name = topographies_in_dir + island_name |
---|
| 60 | island_in_dir_name1 = topographies_in_dir + island_name1 |
---|
| 61 | island_in_dir_name2 = topographies_in_dir + island_name2 |
---|
| 62 | island_in_dir_name3 = topographies_in_dir + island_name3 |
---|
[3828] | 63 | |
---|
[4147] | 64 | coast_in_dir_name = topographies_in_dir + coast_name |
---|
| 65 | offshore_in_dir_name = topographies_in_dir + offshore_name |
---|
| 66 | offshore1_in_dir_name = topographies_in_dir + offshore1_name |
---|
[3828] | 67 | |
---|
[4147] | 68 | onshore_dir_name = topographies_dir + onshore_name |
---|
| 69 | island_dir_name = topographies_dir + island_name |
---|
| 70 | island_dir_name1 = topographies_dir + island_name1 |
---|
| 71 | island_dir_name2 = topographies_dir + island_name2 |
---|
| 72 | island_dir_name3 = topographies_dir + island_name3 |
---|
[3627] | 73 | |
---|
[3828] | 74 | coast_dir_name = topographies_dir + coast_name |
---|
| 75 | offshore_dir_name = topographies_dir + offshore_name |
---|
[3627] | 76 | |
---|
[4147] | 77 | #final topo files |
---|
| 78 | combined_dir_name = topographies_dir + combined_name |
---|
| 79 | combined_time_dir_name = topographies_time_dir + combined_name |
---|
| 80 | combined_smaller_name_dir = topographies_dir + combined_smaller_name |
---|
| 81 | #combined_time_dir_final_name = topographies_time_dir + combined_final_name |
---|
[3940] | 82 | |
---|
[4147] | 83 | meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep |
---|
| 84 | meshes_dir_name = meshes_dir + scenario_name |
---|
[3940] | 85 | |
---|
[4147] | 86 | polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep |
---|
| 87 | tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep |
---|
[3940] | 88 | |
---|
[4147] | 89 | boundaries_source = '????' |
---|
| 90 | #boundaries locations |
---|
| 91 | boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep |
---|
| 92 | boundaries_in_dir_name = boundaries_in_dir + scenario_name |
---|
| 93 | boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep |
---|
| 94 | boundaries_dir_name = boundaries_dir + scenario_name |
---|
| 95 | #boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep |
---|
| 96 | #boundaries_time_dir_name = boundaries_time_dir + boundaries_name #Used by post processing |
---|
[3627] | 97 | |
---|
[4147] | 98 | #output locations |
---|
| 99 | output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep |
---|
| 100 | output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep |
---|
| 101 | output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep |
---|
[3851] | 102 | output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing |
---|
[3828] | 103 | |
---|
[4147] | 104 | #gauges |
---|
| 105 | gauge_name = 'dampier.csv' |
---|
| 106 | gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep |
---|
| 107 | gauges_dir_name = gauges_dir + gauge_name |
---|
[3851] | 108 | |
---|
| 109 | |
---|
[4147] | 110 | ############################### |
---|
| 111 | # Domain definitions |
---|
| 112 | ############################### |
---|
[3627] | 113 | |
---|
| 114 | refzone = 50 |
---|
| 115 | south = degminsec2decimal_degrees(-20,55,0) |
---|
[3835] | 116 | north = degminsec2decimal_degrees(-20,15,0) |
---|
[3627] | 117 | west = degminsec2decimal_degrees(116,17,0) |
---|
[3836] | 118 | east = degminsec2decimal_degrees(117,10,0) |
---|
[3627] | 119 | |
---|
| 120 | p0 = [south, degminsec2decimal_degrees(116,32,0)] |
---|
| 121 | p1 = [south, west] |
---|
| 122 | p2 = [degminsec2decimal_degrees(-20,23,0), west] |
---|
| 123 | p3 = [north, degminsec2decimal_degrees(116,45,0)] |
---|
| 124 | p4 = [north, degminsec2decimal_degrees(117,0,0)] |
---|
| 125 | p5 = [p2[0], degminsec2decimal_degrees(117,8,0)] |
---|
| 126 | p6 = [degminsec2decimal_degrees(-20,30,0), east] |
---|
| 127 | p7 = [degminsec2decimal_degrees(-20,38,0), east] |
---|
| 128 | p8 = [south, east] |
---|
| 129 | |
---|
[4147] | 130 | poly_all, zone = convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8]) |
---|
[3627] | 131 | refzone = zone |
---|
[4147] | 132 | print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0 |
---|
[3627] | 133 | |
---|
[4147] | 134 | res_poly_all = 100000 |
---|
[3744] | 135 | |
---|
[4147] | 136 | ############################### |
---|
| 137 | # Interior region definitions |
---|
| 138 | ############################### |
---|
| 139 | |
---|
| 140 | poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20_pts.csv') |
---|
| 141 | res_pos20_neg20 = 20000 |
---|
| 142 | |
---|
| 143 | poly_dampier = read_polygon(polygons_dir+'dampier_pts.csv') |
---|
| 144 | res_dampier = 500 |
---|
| 145 | |
---|
| 146 | poly_karratha = read_polygon(polygons_dir+'karratha_pts.csv') |
---|
| 147 | res_karratha = 500 |
---|
| 148 | |
---|
| 149 | poly_delambre = read_polygon(polygons_dir+'delambre_pts.csv') |
---|
| 150 | res_delambre = 1000 |
---|
| 151 | |
---|
| 152 | poly_mainisland = read_polygon(polygons_dir+'mainisland_pts.csv') |
---|
| 153 | res_mainisland = 1000 |
---|
| 154 | |
---|
| 155 | poly_NWislands = read_polygon(polygons_dir+'NWislands_pts.csv') |
---|
| 156 | res_NWislands = 1000 |
---|
| 157 | |
---|
| 158 | plot_polygons([poly_pos20_neg20,poly_dampier,poly_karratha,poly_delambre,polylmainisland, |
---|
| 159 | polyNWislands,poly_all],output_run_time_dir + 'poly_pic') |
---|
| 160 | |
---|
| 161 | interior_regions = [[poly_pos20_neg20,res_pos20_neg20],[poly_dampier,res_dampier], |
---|
| 162 | [poly_karratha,res_karratha],[poly_delambre,res_delambre], |
---|
| 163 | [poly_mainisland,res_mainisland],[poly_NWislands,res_NWislands]] |
---|
| 164 | |
---|
| 165 | trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all) |
---|
| 166 | |
---|
| 167 | print 'min number triangles', trigs_min |
---|
| 168 | |
---|
| 169 | ################################################################### |
---|
| 170 | # Clipping regions for export to asc and regions for clipping data |
---|
| 171 | ################################################################### |
---|
| 172 | |
---|
| 173 | # exporting asc grid - Dampier |
---|
| 174 | e_min_area = 474000 |
---|
| 175 | e_max_area = 480000 |
---|
| 176 | n_min_area = 7719000 |
---|
| 177 | n_max_area = 7725000 |
---|
| 178 | |
---|
| 179 | # exporting asc grid - Karratha |
---|
| 180 | e_min_area = |
---|
| 181 | e_max_area = |
---|
| 182 | n_min_area = |
---|
| 183 | n_max_area = |
---|
| 184 | |
---|
| 185 | """ |
---|
| 186 | # used in the CIPMA 2006 scenario |
---|
[3744] | 187 | # CIPMA point of interest |
---|
| 188 | cipma_latitude = -20.588456 |
---|
| 189 | cipma_longitude = 116.771527 |
---|
| 190 | |
---|
| 191 | k0 = [cipma_latitude-0.02, cipma_longitude-0.02] |
---|
| 192 | k1 = [cipma_latitude-0.02, cipma_longitude+0.02] |
---|
| 193 | k2 = [cipma_latitude+0.02, cipma_longitude+0.02] |
---|
| 194 | k3 = [cipma_latitude+0.02, cipma_longitude-0.02] |
---|
| 195 | |
---|
| 196 | cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3]) |
---|
| 197 | assert zone == refzone |
---|
| 198 | |
---|
[3905] | 199 | poly_facility = read_polygon(polygons_dir+'facility.csv') |
---|
[3997] | 200 | poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv') |
---|
| 201 | poly_interior = read_polygon(polygons_dir+'interior.csv') |
---|
| 202 | poly_coast = read_polygon(polygons_dir+'coast_final.csv') |
---|
[3940] | 203 | clip_poly_e = read_polygon(polygons_dir+'gap_e.csv') |
---|
| 204 | clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv') |
---|
| 205 | clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs') |
---|
| 206 | clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs') |
---|
| 207 | |
---|
[4147] | 208 | # exporting asc grid |
---|
| 209 | e_min_area = 474000 |
---|
| 210 | e_max_area = 480000 |
---|
| 211 | n_min_area = 7719000 |
---|
| 212 | n_max_area = 7725000 |
---|
| 213 | |
---|
| 214 | # used in the original 2005 scenario |
---|
[3744] | 215 | #Interior regions |
---|
[3627] | 216 | karratha_south = degminsec2decimal_degrees(-20,44,0) |
---|
| 217 | karratha_north = degminsec2decimal_degrees(-20,42,0) |
---|
| 218 | karratha_west = degminsec2decimal_degrees(116,48,0) |
---|
| 219 | karratha_east = degminsec2decimal_degrees(116,53,30) |
---|
| 220 | |
---|
| 221 | k0 = [karratha_south, karratha_west] |
---|
| 222 | k1 = [karratha_south, karratha_east] |
---|
| 223 | k2 = [karratha_north, karratha_east] |
---|
| 224 | k3 = [karratha_north, karratha_west] |
---|
| 225 | |
---|
[3744] | 226 | karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3]) |
---|
[3716] | 227 | assert zone == refzone |
---|
[3627] | 228 | |
---|
| 229 | #Interior regions |
---|
| 230 | dampier_south = degminsec2decimal_degrees(-20,40,0) |
---|
| 231 | dampier_north = degminsec2decimal_degrees(-20,38,10) |
---|
| 232 | dampier_west = degminsec2decimal_degrees(116,43,0) |
---|
| 233 | dampier_east = degminsec2decimal_degrees(116,45,0) |
---|
| 234 | |
---|
| 235 | d0 = [dampier_south, dampier_west] |
---|
| 236 | d1 = [dampier_south, dampier_east] |
---|
| 237 | d2 = [dampier_north, dampier_east] |
---|
| 238 | d3 = [dampier_north, dampier_west] |
---|
| 239 | |
---|
[3744] | 240 | dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3]) |
---|
[3716] | 241 | assert zone == refzone |
---|
[3627] | 242 | |
---|
| 243 | #Interior regions |
---|
| 244 | refinery_south = degminsec2decimal_degrees(-20,37,50) |
---|
| 245 | refinery_north = degminsec2decimal_degrees(-20,36,0) |
---|
| 246 | refinery_west = degminsec2decimal_degrees(116,44,0) |
---|
| 247 | refinery_east = degminsec2decimal_degrees(116,46,10) |
---|
| 248 | |
---|
| 249 | d0 = [refinery_south, refinery_west] |
---|
| 250 | d1 = [refinery_south, refinery_east] |
---|
| 251 | d2 = [refinery_north, refinery_east] |
---|
| 252 | d3 = [refinery_north, refinery_west] |
---|
| 253 | |
---|
[3744] | 254 | refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3]) |
---|
[3716] | 255 | assert zone == refzone |
---|
[3627] | 256 | |
---|
| 257 | #Interior region around 468899, 7715177: |
---|
| 258 | #lat (-20, 39, 44.93753), lon (116, 42, 5.09106) |
---|
| 259 | |
---|
| 260 | point_south = degminsec2decimal_degrees(-20,39,46) |
---|
| 261 | point_north = degminsec2decimal_degrees(-20,39,42) |
---|
| 262 | point_west = degminsec2decimal_degrees(116,42,0) |
---|
| 263 | point_east = degminsec2decimal_degrees(116,42,10) |
---|
| 264 | |
---|
| 265 | d0 = [point_south, point_west] |
---|
| 266 | d1 = [point_south, point_east] |
---|
| 267 | d2 = [point_north, point_east] |
---|
| 268 | d3 = [point_north, point_west] |
---|
| 269 | |
---|
[3744] | 270 | point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3]) |
---|
[3716] | 271 | assert zone == refzone |
---|
[3627] | 272 | |
---|
| 273 | #Neils areas around interesting points |
---|
| 274 | neil1_point1 = [degminsec2decimal_degrees(-20,35,34), |
---|
| 275 | degminsec2decimal_degrees(116,45,18)] |
---|
| 276 | neil1_point2 = [degminsec2decimal_degrees(-20,36,15), |
---|
| 277 | degminsec2decimal_degrees(116,46,18)] |
---|
| 278 | neil1_point3 = [degminsec2decimal_degrees(-20,35,9), |
---|
| 279 | degminsec2decimal_degrees(116,47,17)] |
---|
| 280 | neil1_point4 = [degminsec2decimal_degrees(-20,34,26), |
---|
| 281 | degminsec2decimal_degrees(116,46,17)] |
---|
| 282 | |
---|
[3744] | 283 | neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1, |
---|
| 284 | neil1_point2, |
---|
| 285 | neil1_point3, |
---|
| 286 | neil1_point4]) |
---|
[3716] | 287 | assert zone == refzone |
---|
[3627] | 288 | |
---|
| 289 | neil2_point1 = [degminsec2decimal_degrees(-20,39,36), |
---|
| 290 | degminsec2decimal_degrees(116,41,33)] |
---|
| 291 | neil2_point2 = [degminsec2decimal_degrees(-20,40,10), |
---|
| 292 | degminsec2decimal_degrees(116,42,13)] |
---|
| 293 | neil2_point3 = [degminsec2decimal_degrees(-20,38,39), |
---|
| 294 | degminsec2decimal_degrees(116,43,49)] |
---|
| 295 | neil2_point4 = [degminsec2decimal_degrees(-20,38,5), |
---|
| 296 | degminsec2decimal_degrees(116,43,9)] |
---|
| 297 | |
---|
[3744] | 298 | neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1, |
---|
| 299 | neil2_point2, |
---|
| 300 | neil2_point3, |
---|
| 301 | neil2_point4]) |
---|
[3716] | 302 | assert zone == refzone |
---|
[3627] | 303 | |
---|
| 304 | #Withnell bay |
---|
| 305 | wb_point1 = [degminsec2decimal_degrees(-20,35,34), |
---|
| 306 | degminsec2decimal_degrees(116,45,18)] |
---|
| 307 | wb_point2 = [degminsec2decimal_degrees(-20,36,15), |
---|
| 308 | degminsec2decimal_degrees(116,46,18)] |
---|
| 309 | wb_point3 = [degminsec2decimal_degrees(-20,35,9), |
---|
| 310 | degminsec2decimal_degrees(116,47,17)] |
---|
| 311 | wb_point4 = [degminsec2decimal_degrees(-20,34,26), |
---|
| 312 | degminsec2decimal_degrees(116,46,17)] |
---|
| 313 | |
---|
[3744] | 314 | wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2, |
---|
| 315 | wb_point3, wb_point4]) |
---|
[3716] | 316 | assert zone == refzone |
---|
[3627] | 317 | |
---|
| 318 | #Larger Withnell bay |
---|
| 319 | lwb_point1 = [degminsec2decimal_degrees(-20,35,59), |
---|
| 320 | degminsec2decimal_degrees(116,42,00)] |
---|
| 321 | lwb_point2 = [degminsec2decimal_degrees(-20,36,50), |
---|
| 322 | degminsec2decimal_degrees(116,46,50)] |
---|
| 323 | lwb_point3 = [degminsec2decimal_degrees(-20,34,00), |
---|
| 324 | degminsec2decimal_degrees(116,47,39)] |
---|
| 325 | lwb_point4 = [degminsec2decimal_degrees(-20,33,00), |
---|
| 326 | degminsec2decimal_degrees(116,42,50)] |
---|
| 327 | |
---|
[3744] | 328 | lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2, |
---|
| 329 | lwb_point3, lwb_point4]) |
---|
[3716] | 330 | |
---|
| 331 | assert zone == refzone |
---|
[4147] | 332 | """ |
---|