1 | """Common filenames and locations for topographic data, meshes and outputs. |
---|
2 | """ |
---|
3 | |
---|
4 | import sys |
---|
5 | from os import sep, environ, getenv, getcwd |
---|
6 | from os.path import expanduser, basename |
---|
7 | |
---|
8 | from anuga.coordinate_transforms.redfearn import\ |
---|
9 | degminsec2decimal_degrees,\ |
---|
10 | convert_from_latlon_to_utm |
---|
11 | |
---|
12 | from time import localtime, strftime, gmtime, ctime |
---|
13 | from anuga.geospatial_data.geospatial_data import * |
---|
14 | from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area |
---|
15 | from anuga.utilities.system_tools import get_user_name |
---|
16 | |
---|
17 | # file and system info |
---|
18 | #--------------------------------- |
---|
19 | codename = 'project.py' |
---|
20 | |
---|
21 | home = getenv('INUNDATIONHOME') #Sandpit's parent dir |
---|
22 | user = get_user_name() |
---|
23 | |
---|
24 | # INUNDATIONHOME is the inundation directory, not the data directory. |
---|
25 | home += sep +'data' |
---|
26 | #---------------------------------- |
---|
27 | # Location and naming of scenario data |
---|
28 | #---------------------------------- |
---|
29 | state = 'western_australia' |
---|
30 | scenario_name = 'dampier_tsunami' |
---|
31 | scenario_datas_name = 'dampier_tsunami_scenario_2006' #name of the directory where the data is stored |
---|
32 | #scenario_datas_name = 'karratha_tsunami_scenario_2005' # Tmp location |
---|
33 | |
---|
34 | #mesh_name = 'elevation50m' |
---|
35 | boundaries_name = 'dampier' |
---|
36 | boundaries_source = 'mag_9_corrected' |
---|
37 | #boundaries_source = 'test' |
---|
38 | |
---|
39 | tide = 2.4 |
---|
40 | #tide = 0.0 |
---|
41 | |
---|
42 | # topography file names |
---|
43 | onshore_name = 'dli_no_islands' |
---|
44 | coast_name = 'DTED_05_Contour' |
---|
45 | islands_name = 'dted_islands' |
---|
46 | offshore_name = 'XY100003902' |
---|
47 | offshore_name1 = 'XY100003903' |
---|
48 | offshore_name2 = 'XY100003951' |
---|
49 | offshore_name3 = 'XY100006321' |
---|
50 | offshore_name4 = 'XY100011756' |
---|
51 | offshore_name5 = 'XY100014243' |
---|
52 | offshore_name6 = 'XY100014244' |
---|
53 | offshore_name7 = 'XY100021081' |
---|
54 | offshore_name8 = 'XY100021082' |
---|
55 | offshore_name9 = 'XY100021083' |
---|
56 | offshore_name10 = 'XY100021085' |
---|
57 | offshore_name11 = 'XY100021086' |
---|
58 | offshore_name12 = 'XY100026309' |
---|
59 | offshore_name13 = 'XY100026338' |
---|
60 | offshore_name14 = 'XYDM83' |
---|
61 | |
---|
62 | offshore_old = 'elevation50m' |
---|
63 | |
---|
64 | combined_name ='dampier_combined_elevation' |
---|
65 | combined_final_name ='dampier_combined_elevation_final' |
---|
66 | |
---|
67 | gauge_name = 'dampier_gauges_up2.csv' |
---|
68 | |
---|
69 | #Derive subdirectories and filenames |
---|
70 | |
---|
71 | meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep |
---|
72 | topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep |
---|
73 | gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep |
---|
74 | polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep |
---|
75 | boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep |
---|
76 | #outputdir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'output'+sep |
---|
77 | tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep |
---|
78 | |
---|
79 | time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir |
---|
80 | gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir |
---|
81 | #cctime = strftime('%Y%m%d_%H%M%S',ctime()) #gets time for new dir |
---|
82 | build_time = time+'_build' |
---|
83 | run_time = time+'_run' |
---|
84 | |
---|
85 | print 'gtime: ', gtime |
---|
86 | output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep |
---|
87 | output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep |
---|
88 | output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep |
---|
89 | output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep |
---|
90 | topographies_time_dir = topographies_dir+build_time+sep |
---|
91 | boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep |
---|
92 | boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep |
---|
93 | meshes_time_dir = meshes_dir+build_time+sep |
---|
94 | |
---|
95 | #ideas |
---|
96 | #boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep |
---|
97 | |
---|
98 | gauge_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep |
---|
99 | gauge_filename = gauge_dir + 'dampier_gauges_up2.csv' |
---|
100 | |
---|
101 | gauges_dir_name = gauges_dir + gauge_name |
---|
102 | |
---|
103 | onshore_dir_name = topographies_dir + onshore_name |
---|
104 | coast_dir_name = topographies_dir + coast_name |
---|
105 | islands_dir_name = topographies_dir + islands_name |
---|
106 | offshore_dir_name = topographies_dir + offshore_name |
---|
107 | offshore_dir_name1 = topographies_dir + offshore_name1 |
---|
108 | offshore_dir_name2 = topographies_dir + offshore_name2 |
---|
109 | offshore_dir_name3 = topographies_dir + offshore_name3 |
---|
110 | offshore_dir_name4 = topographies_dir + offshore_name4 |
---|
111 | offshore_dir_name5 = topographies_dir + offshore_name5 |
---|
112 | offshore_dir_name6 = topographies_dir + offshore_name6 |
---|
113 | offshore_dir_name7 = topographies_dir + offshore_name7 |
---|
114 | offshore_dir_name8 = topographies_dir + offshore_name8 |
---|
115 | offshore_dir_name9 = topographies_dir + offshore_name9 |
---|
116 | offshore_dir_name10 = topographies_dir + offshore_name10 |
---|
117 | offshore_dir_name11 = topographies_dir + offshore_name11 |
---|
118 | offshore_dir_name12 = topographies_dir + offshore_name12 |
---|
119 | offshore_dir_name13 = topographies_dir + offshore_name13 |
---|
120 | offshore_dir_name14 = topographies_dir + offshore_name14 |
---|
121 | |
---|
122 | offshore_dir_name_old = topographies_dir + offshore_old |
---|
123 | |
---|
124 | |
---|
125 | |
---|
126 | #output dir |
---|
127 | combined_dir_name = topographies_dir + combined_name |
---|
128 | combined_time_dir_name = topographies_time_dir + combined_name |
---|
129 | combined_time_dir_final_name = topographies_time_dir + combined_final_name |
---|
130 | |
---|
131 | |
---|
132 | meshes_dir_name = meshes_dir + scenario_name |
---|
133 | meshes_time_dir_name = meshes_time_dir + scenario_name |
---|
134 | #output_build_time_dir_name = output_build_time_dir + scenario_name #Used by post processing |
---|
135 | output_run_time_dir_name = output_run_time_dir + scenario_name #Used by post processing |
---|
136 | boundaries_in_dir_name = boundaries_in_dir + boundaries_name |
---|
137 | boundaries_time_dir_name = boundaries_time_dir + boundaries_name #Used by post processing |
---|
138 | boundaries_dir_name = boundaries_dir + boundaries_name |
---|
139 | |
---|
140 | |
---|
141 | |
---|
142 | # Regions |
---|
143 | |
---|
144 | refzone = 50 |
---|
145 | south = degminsec2decimal_degrees(-20,55,0) |
---|
146 | north = degminsec2decimal_degrees(-20,15,0) |
---|
147 | #north = degminsec2decimal_degrees(-19,15,0) |
---|
148 | west = degminsec2decimal_degrees(116,17,0) |
---|
149 | east = degminsec2decimal_degrees(117,10,0) |
---|
150 | |
---|
151 | #only used to clip boundary condition |
---|
152 | #north_boundary = north + 0.02 |
---|
153 | #south_boundary = south - 0.02 |
---|
154 | #west_boundary = west - 0.02 |
---|
155 | #east_boundary = east + 0.02 |
---|
156 | |
---|
157 | south_boundary = degminsec2decimal_degrees(-21,0,0) |
---|
158 | #north_boundary = degminsec2decimal_degrees(-19,00,0) |
---|
159 | north_boundary = degminsec2decimal_degrees(-20,10,0) |
---|
160 | #west_boundary = degminsec2decimal_degrees(116,0,0) |
---|
161 | #east_boundary = degminsec2decimal_degrees(118,00,0) |
---|
162 | west_boundary = degminsec2decimal_degrees(116,10,0) |
---|
163 | east_boundary = degminsec2decimal_degrees(117,15,0) |
---|
164 | |
---|
165 | |
---|
166 | p0 = [south, degminsec2decimal_degrees(116,32,0)] |
---|
167 | p1 = [south, west] |
---|
168 | p2 = [degminsec2decimal_degrees(-20,23,0), west] |
---|
169 | p3 = [north, degminsec2decimal_degrees(116,45,0)] |
---|
170 | p4 = [north, degminsec2decimal_degrees(117,0,0)] |
---|
171 | p5 = [p2[0], degminsec2decimal_degrees(117,8,0)] |
---|
172 | p6 = [degminsec2decimal_degrees(-20,30,0), east] |
---|
173 | p7 = [degminsec2decimal_degrees(-20,38,0), east] |
---|
174 | p8 = [south, east] |
---|
175 | |
---|
176 | bounding_polygon, zone =\ |
---|
177 | convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8]) |
---|
178 | #bounding_polygon, zone =\ |
---|
179 | # convert_from_latlon_to_utm([p1, p2, p3, p4, p5, p6, p7]) |
---|
180 | refzone = zone |
---|
181 | from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon |
---|
182 | print 'poly area', polygon_area(bounding_polygon)/1000000.0 |
---|
183 | |
---|
184 | #Interior regions |
---|
185 | |
---|
186 | # CIPMA point of interest |
---|
187 | cipma_latitude = -20.588456 |
---|
188 | cipma_longitude = 116.771527 |
---|
189 | |
---|
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 | |
---|
199 | e_min_area = 474000 |
---|
200 | e_max_area = 480000 |
---|
201 | n_min_area = 7719000 |
---|
202 | n_max_area = 7725000 |
---|
203 | |
---|
204 | poly_facility = read_polygon(polygons_dir+'facility.csv') |
---|
205 | |
---|
206 | poly_pipeline = read_polygon(polygons_dir+'pipeline2.csv') |
---|
207 | |
---|
208 | poly_interior = read_polygon(polygons_dir+'interior.csv') |
---|
209 | |
---|
210 | poly_coast = read_polygon(polygons_dir+'coast_final.csv') |
---|
211 | |
---|
212 | clip_poly_e = read_polygon(polygons_dir+'gap_e.csv') |
---|
213 | |
---|
214 | clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv') |
---|
215 | |
---|
216 | clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs') |
---|
217 | |
---|
218 | clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs') |
---|
219 | |
---|
220 | #Interior regions |
---|
221 | karratha_south = degminsec2decimal_degrees(-20,44,0) |
---|
222 | karratha_north = degminsec2decimal_degrees(-20,42,0) |
---|
223 | karratha_west = degminsec2decimal_degrees(116,48,0) |
---|
224 | karratha_east = degminsec2decimal_degrees(116,53,30) |
---|
225 | |
---|
226 | k0 = [karratha_south, karratha_west] |
---|
227 | k1 = [karratha_south, karratha_east] |
---|
228 | k2 = [karratha_north, karratha_east] |
---|
229 | k3 = [karratha_north, karratha_west] |
---|
230 | |
---|
231 | karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3]) |
---|
232 | assert zone == refzone |
---|
233 | |
---|
234 | |
---|
235 | #Interior regions |
---|
236 | dampier_south = degminsec2decimal_degrees(-20,40,0) |
---|
237 | dampier_north = degminsec2decimal_degrees(-20,38,10) |
---|
238 | dampier_west = degminsec2decimal_degrees(116,43,0) |
---|
239 | dampier_east = degminsec2decimal_degrees(116,45,0) |
---|
240 | |
---|
241 | d0 = [dampier_south, dampier_west] |
---|
242 | d1 = [dampier_south, dampier_east] |
---|
243 | d2 = [dampier_north, dampier_east] |
---|
244 | d3 = [dampier_north, dampier_west] |
---|
245 | |
---|
246 | dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3]) |
---|
247 | assert zone == refzone |
---|
248 | |
---|
249 | |
---|
250 | #Interior regions |
---|
251 | refinery_south = degminsec2decimal_degrees(-20,37,50) |
---|
252 | refinery_north = degminsec2decimal_degrees(-20,36,0) |
---|
253 | refinery_west = degminsec2decimal_degrees(116,44,0) |
---|
254 | refinery_east = degminsec2decimal_degrees(116,46,10) |
---|
255 | |
---|
256 | d0 = [refinery_south, refinery_west] |
---|
257 | d1 = [refinery_south, refinery_east] |
---|
258 | d2 = [refinery_north, refinery_east] |
---|
259 | d3 = [refinery_north, refinery_west] |
---|
260 | |
---|
261 | refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3]) |
---|
262 | assert zone == refzone |
---|
263 | |
---|
264 | |
---|
265 | #Interior region around 468899, 7715177: |
---|
266 | #lat (-20, 39, 44.93753), lon (116, 42, 5.09106) |
---|
267 | |
---|
268 | point_south = degminsec2decimal_degrees(-20,39,46) |
---|
269 | point_north = degminsec2decimal_degrees(-20,39,42) |
---|
270 | point_west = degminsec2decimal_degrees(116,42,0) |
---|
271 | point_east = degminsec2decimal_degrees(116,42,10) |
---|
272 | |
---|
273 | d0 = [point_south, point_west] |
---|
274 | d1 = [point_south, point_east] |
---|
275 | d2 = [point_north, point_east] |
---|
276 | d3 = [point_north, point_west] |
---|
277 | |
---|
278 | point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3]) |
---|
279 | assert zone == refzone |
---|
280 | |
---|
281 | |
---|
282 | #Neils areas around interesting points |
---|
283 | neil1_point1 = [degminsec2decimal_degrees(-20,35,34), |
---|
284 | degminsec2decimal_degrees(116,45,18)] |
---|
285 | neil1_point2 = [degminsec2decimal_degrees(-20,36,15), |
---|
286 | degminsec2decimal_degrees(116,46,18)] |
---|
287 | neil1_point3 = [degminsec2decimal_degrees(-20,35,9), |
---|
288 | degminsec2decimal_degrees(116,47,17)] |
---|
289 | neil1_point4 = [degminsec2decimal_degrees(-20,34,26), |
---|
290 | degminsec2decimal_degrees(116,46,17)] |
---|
291 | |
---|
292 | neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1, |
---|
293 | neil1_point2, |
---|
294 | neil1_point3, |
---|
295 | neil1_point4]) |
---|
296 | assert zone == refzone |
---|
297 | |
---|
298 | |
---|
299 | |
---|
300 | |
---|
301 | neil2_point1 = [degminsec2decimal_degrees(-20,39,36), |
---|
302 | degminsec2decimal_degrees(116,41,33)] |
---|
303 | neil2_point2 = [degminsec2decimal_degrees(-20,40,10), |
---|
304 | degminsec2decimal_degrees(116,42,13)] |
---|
305 | neil2_point3 = [degminsec2decimal_degrees(-20,38,39), |
---|
306 | degminsec2decimal_degrees(116,43,49)] |
---|
307 | neil2_point4 = [degminsec2decimal_degrees(-20,38,5), |
---|
308 | degminsec2decimal_degrees(116,43,9)] |
---|
309 | |
---|
310 | neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1, |
---|
311 | neil2_point2, |
---|
312 | neil2_point3, |
---|
313 | neil2_point4]) |
---|
314 | assert zone == refzone |
---|
315 | |
---|
316 | |
---|
317 | |
---|
318 | |
---|
319 | |
---|
320 | #Withnell bay |
---|
321 | wb_point1 = [degminsec2decimal_degrees(-20,35,34), |
---|
322 | degminsec2decimal_degrees(116,45,18)] |
---|
323 | wb_point2 = [degminsec2decimal_degrees(-20,36,15), |
---|
324 | degminsec2decimal_degrees(116,46,18)] |
---|
325 | wb_point3 = [degminsec2decimal_degrees(-20,35,9), |
---|
326 | degminsec2decimal_degrees(116,47,17)] |
---|
327 | wb_point4 = [degminsec2decimal_degrees(-20,34,26), |
---|
328 | degminsec2decimal_degrees(116,46,17)] |
---|
329 | |
---|
330 | wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2, |
---|
331 | wb_point3, wb_point4]) |
---|
332 | assert zone == refzone |
---|
333 | |
---|
334 | |
---|
335 | |
---|
336 | |
---|
337 | |
---|
338 | #Larger Withnell bay |
---|
339 | lwb_point1 = [degminsec2decimal_degrees(-20,35,59), |
---|
340 | degminsec2decimal_degrees(116,42,00)] |
---|
341 | lwb_point2 = [degminsec2decimal_degrees(-20,36,50), |
---|
342 | degminsec2decimal_degrees(116,46,50)] |
---|
343 | lwb_point3 = [degminsec2decimal_degrees(-20,34,00), |
---|
344 | degminsec2decimal_degrees(116,47,39)] |
---|
345 | lwb_point4 = [degminsec2decimal_degrees(-20,33,00), |
---|
346 | degminsec2decimal_degrees(116,42,50)] |
---|
347 | |
---|
348 | lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2, |
---|
349 | lwb_point3, lwb_point4]) |
---|
350 | |
---|
351 | assert zone == refzone |
---|
352 | |
---|
353 | |
---|
354 | |
---|
355 | |
---|