source: anuga_work/production/perth_2006/project.py @ 3908

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

adding script for locating indigeneous communities for FESA nominated areas

File size: 6.4 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 = 'perth_tsunami_scenario_2006'
27
28# onshore data provided by WA DLI
29onshore_name = '' # 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
58boundary_basename = 'SU-AU' # Mw ?
59
60#swollen/ all data output
61basename = 'source'
62codename = 'project.py'
63
64#Derive subdirectories and filenames
65local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
66meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
67datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
68gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
69polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
70boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
71outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
72outputtimedir = outputdir + local_time + sep
73polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
74
75gauge_filename = gaugedir + 'perth_gauges.csv'
76community_filename = gaugedir + 'CHINS_v2.csv'
77community_perth = gaugedir + 'community_perth.csv'
78codedir = getcwd()+sep                           
79codedirname = codedir + 'project.py'
80meshname = outputtimedir + 'mesh_' + basename
81
82# Necessary if using point datasets, rather than grid
83onshore_dem_name = datadir + onshore_name
84offshore_dem_name1 = datadir + offshore_name1
85offshore_dem_name2 = datadir + offshore_name2
86offshore_dem_name3 = datadir + offshore_name3
87offshore_dem_name4 = datadir + offshore_name4
88offshore_dem_name5 = datadir + offshore_name5
89offshore_dem_name6 = datadir + offshore_name6
90offshore_dem_name7 = datadir + offshore_name7
91offshore_dem_name8 = datadir + offshore_name8
92offshore_dem_name9 = datadir + offshore_name9
93offshore_dem_name10 = datadir + offshore_name10
94offshore_dem_name11 = datadir + offshore_name11
95offshore_dem_name12 = datadir + offshore_name12
96offshore_dem_name13 = datadir + offshore_name13
97offshore_dem_name14 = datadir + offshore_name14
98offshore_dem_name15 = datadir + offshore_name15
99offshore_dem_name16 = datadir + offshore_name16
100offshore_dem_name17 = datadir + offshore_name17
101offshore_dem_name18 = datadir + offshore_name18
102offshore_dem_name19 = datadir + offshore_name19
103offshore_dem_name20 = datadir + offshore_name20
104offshore_dem_name21 = datadir + offshore_name21
105offshore_dem_name22 = datadir + offshore_name22
106
107coast_dem_name = datadir + coast_name
108
109combined_dem_name   = datadir + 'perth_combined_elevation'
110
111###############################
112# Domain definitions
113###############################
114
115# bounding box for clipping MOST/URS output (much bigger than study area)
116##south = degminsec2decimal_degrees(-19,0,0)
117##north = degminsec2decimal_degrees(-17,15,0)
118##west  = degminsec2decimal_degrees(120,0,0)
119##east  = degminsec2decimal_degrees(122,30,0)
120##
121##d0 = [south, west]
122##d1 = [south, east]
123##d2 = [north, east]
124##d3 = [north, west]
125##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
126##refzone = zone
127
128# bounding polygon for study area
129polyAll = read_polygon(polygondir+'extent.csv')
130
131# plot bounding polygon and make sure BC info surrounds it
132#plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
133print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
134
135###################################################################
136# Clipping regions for export to asc and regions for clipping data
137###################################################################
138
139# exporting asc grid
140eastingmin = 406215.87
141eastingmax = 440208.78
142northingmin = 7983427.73
143northingmax = 8032834.52
144
145###############################
146# Interior region definitions
147###############################
148
149# perth digitized polygons
150poly_perth1 = read_polygon(polygondir+'perth_Local_Polygon_update.csv')
151poly_perth2 = read_polygon(polygondir+'perth_Close2_update.csv')
152poly_perth3 = read_polygon(polygondir+'perth_Coast_update.csv')
153#poly_perth4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
154
155plot_polygons([polyAll,poly_perth1,poly_perth2,poly_perth3],'boundingpoly2',verbose=False)
156print 'Area of local polygon', polygon_area(poly_perth1)/1000000.0
157print 'Area of close polygon', polygon_area(poly_perth2)/1000000.0
158print 'Area of coastal polygon', polygon_area(poly_perth3)/1000000.0
159#print 'Area of cable beach polygon', polygon_area(poly_perth4)/1000000.0
160
161for i in poly_perth3:
162    v = is_inside_polygon(i,poly_perth1, verbose=False)
163    if v == False: print v
164
165def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
166    from anuga.utilities.polygon import polygon_area
167   
168    # TO DO check if any of the regions fall inside one another
169    no_triangles = 0.0
170    area = polygon_area(bounding_poly)
171    for i,j in interior_regions:
172        this_area = polygon_area(i)
173        no_triangles += this_area/j
174        area -= this_area
175        print j, this_area/1000000., area/1000000.
176    no_triangles += area/remainder_res
177    return int(no_triangles/0.7)
Note: See TracBrowser for help on using the repository browser.