source: anuga_work/production/broome_2006/project.py @ 3972

Last change on this file since 3972 was 3972, checked in by sexton, 17 years ago

(i) incorporating new supply of interpolated data for Broome (ii) updating report to look at MOST versus ANUGA for Hobart scenario

File size: 6.6 KB
Line 
1# -*- coding: cp1252 -*-
2"""Common filenames and locations for topographic data, meshes and outputs.
3"""
4
5from os import sep, environ, getenv, getcwd
6from os.path import expanduser
7import sys
8from time import localtime, strftime, gmtime
9from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
10#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
11
12if sys.platform == 'win32':
13    home = getenv('INUNDATIONHOME')
14    user = getenv('USERPROFILE')
15
16else:   
17    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
18    user = getenv('LOGNAME')
19    print 'USER:', user
20
21# INUNDATIONHOME is the inundation directory, not the data directory.
22home += sep +'data'
23
24#Making assumptions about the location of scenario data
25state = 'western_australia'
26scenario_dir_name = 'broome_tsunami_scenario_2006'
27
28# onshore data provided by WA DLI
29onshore_name = 'dted2_z51' # original
30
31# AHO + DPI data
32offshore_name1 = 'XY100011610'
33offshore_name2 = 'XY100011611'
34offshore_name3 = 'XY100011613'
35offshore_name4 = 'XY100011614'
36offshore_name5 = 'XY100011616'
37offshore_name6 = 'XY100011617'
38offshore_name7 = 'XY100011618'
39offshore_name8 = 'XY100011621'
40offshore_name9 = 'XY100011623'
41offshore_name10 = 'XY100011745'
42offshore_name11 = 'XY100011746'
43offshore_name12 = 'XY100017530'
44offshore_name13 = 'XY100017532'
45offshore_name14 = 'XY100017538'
46offshore_name15 = 'XY100017540'
47offshore_name16 = 'XYBR66'
48offshore_name17 = 'XYBR70'
49offshore_name18 = 'XYBR80'
50offshore_name19 = 'XYBR88'
51offshore_name20 = 'XYBR93'
52offshore_name21 = 'XYBR0110'
53offshore_name22 = 'XYWADPI'
54
55# developed by NM&I
56coast_name = 'coastline'
57
58# TIN model to fill in data gap
59bathy_interp = 'nearest_neighbour' #'interpolate'
60
61boundary_basename = 'SU-AU' # Mw ?
62
63#swollen/ all data output
64basename = 'source'
65codename = 'project.py'
66
67#Derive subdirectories and filenames
68local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
69meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
70datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
71gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
72polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
73boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
74outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
75outputtimedir = outputdir + local_time + sep
76polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
77
78gauge_filename = gaugedir + 'broome_gauges.csv'
79community_filename = gaugedir + 'CHINS_v2.csv'
80community_broome = gaugedir + 'community_broome.csv'
81codedir = getcwd()+sep                           
82codedirname = codedir + 'project.py'
83meshname = outputtimedir + 'mesh_' + basename
84
85# Necessary if using point datasets, rather than grid
86onshore_dem_name = datadir + onshore_name
87offshore_interp_dem_name = datadir + bathy_interp
88offshore_dem_name1 = datadir + offshore_name1
89offshore_dem_name2 = datadir + offshore_name2
90offshore_dem_name3 = datadir + offshore_name3
91offshore_dem_name4 = datadir + offshore_name4
92offshore_dem_name5 = datadir + offshore_name5
93offshore_dem_name6 = datadir + offshore_name6
94offshore_dem_name7 = datadir + offshore_name7
95offshore_dem_name8 = datadir + offshore_name8
96offshore_dem_name9 = datadir + offshore_name9
97offshore_dem_name10 = datadir + offshore_name10
98offshore_dem_name11 = datadir + offshore_name11
99offshore_dem_name12 = datadir + offshore_name12
100offshore_dem_name13 = datadir + offshore_name13
101offshore_dem_name14 = datadir + offshore_name14
102offshore_dem_name15 = datadir + offshore_name15
103offshore_dem_name16 = datadir + offshore_name16
104offshore_dem_name17 = datadir + offshore_name17
105offshore_dem_name18 = datadir + offshore_name18
106offshore_dem_name19 = datadir + offshore_name19
107offshore_dem_name20 = datadir + offshore_name20
108offshore_dem_name21 = datadir + offshore_name21
109offshore_dem_name22 = datadir + offshore_name22
110
111coast_dem_name = datadir + coast_name
112
113combined_dem_name   = datadir + 'broome_combined_elevation'
114
115###############################
116# Domain definitions
117###############################
118
119# bounding box for clipping MOST/URS output (much bigger than study area)
120##south = degminsec2decimal_degrees(-19,0,0)
121##north = degminsec2decimal_degrees(-17,15,0)
122##west  = degminsec2decimal_degrees(120,0,0)
123##east  = degminsec2decimal_degrees(122,30,0)
124##
125##d0 = [south, west]
126##d1 = [south, east]
127##d2 = [north, east]
128##d3 = [north, west]
129##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
130##refzone = zone
131
132# bounding polygon for study area
133polyAll = read_polygon(polygondir+'extent.csv')
134
135# plot bounding polygon and make sure BC info surrounds it
136#plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
137print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
138
139###################################################################
140# Clipping regions for export to asc and regions for clipping data
141###################################################################
142
143# exporting asc grid
144eastingmin = 412036.43
145eastingmax = 427184.37
146northingmin = 8007984
147northingmax = 8029830
148
149       
150       
151
152
153###############################
154# Interior region definitions
155###############################
156
157# broome digitized polygons
158poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon_update.csv')
159poly_broome2 = read_polygon(polygondir+'Broome_Close2_update.csv')
160poly_broome3 = read_polygon(polygondir+'Broome_Coast_update.csv')
161#poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
162
163plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)
164print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
165print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
166print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
167#print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0
168
169for i in poly_broome3:
170    v = is_inside_polygon(i,poly_broome1, verbose=False)
171    if v == False: print v
172
173def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
174    from anuga.utilities.polygon import polygon_area
175   
176    # TO DO check if any of the regions fall inside one another
177    no_triangles = 0.0
178    area = polygon_area(bounding_poly)
179    for i,j in interior_regions:
180        this_area = polygon_area(i)
181        no_triangles += this_area/j
182        area -= this_area
183        print j, this_area/1000000., area/1000000.
184    no_triangles += area/remainder_res
185    return int(no_triangles/0.7)
Note: See TracBrowser for help on using the repository browser.