Changeset 4094

Timestamp:

Dec 19, 2006, 3:30:10 PM (17 years ago)

Author:

nick

Message:

update busselton

File:

• anuga_work/production/busselton_2006/project.py (modified) (4 diffs)

anuga_work/production/busselton_2006/project.py

@@ -8 +8 @@
 from time import localtime, strftime, gmtime
 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
+from anuga.utilities.system_tools import get_user_name
 
-if sys.platform == 'win32':
-    home = getenv('INUNDATIONHOME')
-    user = getenv('USERPROFILE')
+# file and system info
+#---------------------------------
+codename = 'project.py'
 
-else:    
-    home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
-    user = getenv('LOGNAME')
-    print 'USER:', user
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+user = get_user_name()
 
 # INUNDATIONHOME is the inundation directory, not the data directory.
 home += sep +'data'
 
+#time stuff
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
+build_time = time+'_build'
+run_time = time+'_run'
+print 'gtime: ', gtime
+
+tide = 0.6
+
 #Making assumptions about the location of scenario data
 state = 'western_australia'
-scenario_dir_name = 'busselton_tsunami_scenario_2006'
+scenario_name = 'busselton'
+scenario = 'busselton_tsunami_scenario_2006'
 
 # onshore data provided by WA DLI
-onshore_name =  ''# original
+onshore_name = 'DLI_orthophoto_DEM' # original
+#islands
 
 # AHO + DPI data
-offshore_name1 = 'XY100011610'
-offshore_name2 = 'XY100011611'
-offshore_name3 = 'XY100011613'
-offshore_name4 = 'XY100011614'
-offshore_name5 = 'XY100011616'
-offshore_name6 = 'XY100011617'
-offshore_name7 = 'XY100011618'
-offshore_name8 = 'XY100011621'
-offshore_name9 = 'XY100011623'
-offshore_name10 = 'XY100011745'
-offshore_name11 = 'XY100011746'
-offshore_name12 = 'XY100017530'
-offshore_name13 = 'XY100017532'
-offshore_name14 = 'XY100017538'
-offshore_name15 = 'XY100017540'
-offshore_name16 = 'XYBR66'
-offshore_name17 = 'XYBR70'
-offshore_name18 = 'XYBR80'
-offshore_name19 = 'XYBR88'
-offshore_name20 = 'XYBR93'
-offshore_name21 = 'XYBR0110'
-offshore_name22 = 'XYWADPI'
+coast_name = '100k_coastline'
+offshore_name = 'Bathymetry'
 
-# developed by NM&I
-coast_name = 'coastline'
+#final topo name
+combined_name ='busselton_combined_elevation'
+combined_smaller_name = 'busselton_combined_elevation_smaller'
 
-boundary_basename = 'SU-AU' # Mw ?
 
-#swollen/ all data output
-basename = 'source'
-codename = 'project.py'
+topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
+topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
+topographies_time_dir = topographies_dir+build_time+sep
 
-#Derive subdirectories and filenames
-local_time = strftime('%Y%m%d_%H%M%S',gmtime()) 
-meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
-datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
-gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
-polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
-boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
-outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
-outputtimedir = outputdir + local_time + sep
-polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+#input topo file location
+onshore_in_dir_name = topographies_in_dir + onshore_name
+#island_in_dir_name = topographies_in_dir + island_name
 
-gauge_filename = gaugedir + 'busselton_gauges.csv'
-community_filename = gaugedir + 'CHINS_v2.csv'
-community_broome = gaugedir + 'community_busselton.csv'
-codedir = getcwd()+sep                            
-codedirname = codedir + 'project.py'
-meshname = outputtimedir + 'mesh_' + basename
+coast_in_dir_name = topographies_in_dir + coast_name
+offshore_in_dir_name = topographies_in_dir + offshore_name
 
-# Necessary if using point datasets, rather than grid
-onshore_dem_name = datadir + onshore_name
-offshore_dem_name1 = datadir + offshore_name1
-offshore_dem_name2 = datadir + offshore_name2
-offshore_dem_name3 = datadir + offshore_name3
-offshore_dem_name4 = datadir + offshore_name4
-offshore_dem_name5 = datadir + offshore_name5
-offshore_dem_name6 = datadir + offshore_name6
-offshore_dem_name7 = datadir + offshore_name7
-offshore_dem_name8 = datadir + offshore_name8
-offshore_dem_name9 = datadir + offshore_name9
-offshore_dem_name10 = datadir + offshore_name10
-offshore_dem_name11 = datadir + offshore_name11
-offshore_dem_name12 = datadir + offshore_name12
-offshore_dem_name13 = datadir + offshore_name13
-offshore_dem_name14 = datadir + offshore_name14
-offshore_dem_name15 = datadir + offshore_name15
-offshore_dem_name16 = datadir + offshore_name16
-offshore_dem_name17 = datadir + offshore_name17
-offshore_dem_name18 = datadir + offshore_name18
-offshore_dem_name19 = datadir + offshore_name19
-offshore_dem_name20 = datadir + offshore_name20
-offshore_dem_name21 = datadir + offshore_name21
-offshore_dem_name22 = datadir + offshore_name22
+onshore_dir_name = topographies_dir + onshore_name
+#island_dir_name = topographies_dir + island_name
 
-coast_dem_name = datadir + coast_name
+coast_dir_name = topographies_dir + coast_name
+offshore_dir_name = topographies_dir + offshore_name
 
-combined_dem_name   = datadir + 'broome_combined_elevation'
+#final topo files
+combined_dir_name = topographies_dir + combined_name
+combined_time_dir_name = topographies_time_dir + combined_name
+combined_smaller_dir_name = topographies_dir + combined_smaller_name
+
+meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
+meshes_dir_name = meshes_dir + scenario_name
+
+polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep
+tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep
+
+boundaries_source = '????'
+#boundaries locations
+boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
+boundaries_in_dir_name = boundaries_in_dir + scenario_name
+boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
+boundaries_dir_name = boundaries_dir + scenario_name
+#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
+#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
+
+#output locations
+output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
+output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
+output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
+
+#gauges
+gauge_name = '???.csv'
+gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+gauges_dir_name = gauges_dir + gauge_name
+
 
 ###############################
@@ -113 +102 @@
 
 # bounding box for clipping MOST/URS output (much bigger than study area)
-##south = degminsec2decimal_degrees(-19,0,0)
-##north = degminsec2decimal_degrees(-17,15,0)
-##west  = degminsec2decimal_degrees(120,0,0)
-##east  = degminsec2decimal_degrees(122,30,0) 
-##
-##d0 = [south, west]
-##d1 = [south, east]
-##d2 = [north, east]
-##d3 = [north, west]
-##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
-##refzone = zone
 
 # bounding polygon for study area
-polyAll = read_polygon(polygondir+'extent.csv')
+poly_all = read_polygon(polygons_dir+'poly_all.csv')
+res_poly_all = 100000
 
-# plot bounding polygon and make sure BC info surrounds it
-#plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
-print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
 
+'''
 ###################################################################
 # Clipping regions for export to asc and regions for clipping data
@@ -147 +124 @@
 
 # broome digitized polygons
-poly_busselton1 = read_polygon(polygondir+'buss_Local_Polygon_update.csv')
-poly_busselton2 = read_polygon(polygondir+'buss_Close2_update.csv')
-poly_busselton3 = read_polygon(polygondir+'buss_Coast_update.csv')
+poly_busselton1 = read_polygon(polygons_dir+'buss_Local_Polygon_update.csv')
+res_busselton1 = 500000
+poly_busselton2 = read_polygon(polygons_dir+'buss_Close2_update.csv')
+res_busselton2 = 500000
+poly_busselton3 = read_polygon(polygons_dir+'buss_Coast_update.csv')
+res_busselton3 = 500000
 
-plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)
-print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
-print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
-print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
+#plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)
 
-for i in poly_broome3:
-    v = is_inside_polygon(i,poly_broome1, verbose=False)
-    if v == False: print v
+interior_regions = [[poly_busselton1,res_busselton1],[poly_busselton2,res_busselton2],
+                    [poly_busselton3,res_busselton3]]
+#trigs_all = polygon_area(poly_all)/res_poly_all
+#trigs_bus1 = polygon_area(poly_busselton1)/res_busselton1
+#trigs_bus2 = polygon_area(poly_busselton2)/res_busselton2
+#trigs_bus3 = polygon_area(poly_busselton3)/res_busselton3
+#trigs_min = trigs_bound + trigs_pos + trigs_cbd + trigs_penguin
+
+#print 'Area of bounding poly', trigs_all
+#print 'Area of busselton1', trigs_bus1
+#print 'Area of busselton2', trigs_bus2
+#print 'Area of busselton3', trigs_bus3
+#print 'min number triangles', trigs_min
+
+trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+
+print 'min number triangles', trigs_min
 
 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
@@ -173 +164 @@
     no_triangles += area/remainder_res
     return int(no_triangles/0.7)
+'''
+