# -*- coding: cp1252 -*- """Common filenames and locations for topographic data, meshes and outputs. """ from os import sep, environ, getenv, getcwd from os.path import expanduser import sys 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 if sys.platform == 'win32': home = getenv('INUNDATIONHOME') user = getenv('USERPROFILE') 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 # INUNDATIONHOME is the inundation directory, not the data directory. home += sep +'data' #Making assumptions about the location of scenario data state = 'western_australia' scenario_dir_name = 'broome_tsunami_scenario_2006' # onshore data provided by WA DLI onshore_name = 'dted2_z51' # original # 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' # developed by NM&I coast_name = 'coastline' # TIN model to fill in data gap bathy_interp = 'nearest_neighbour' #'interpolate' boundary_basename = 'SU-AU' # Mw ? #swollen/ all data output basename = 'source' codename = 'project.py' #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 gauge_filename = gaugedir + 'broome_gauges.csv' community_filename = gaugedir + 'CHINS_v2.csv' community_broome = gaugedir + 'community_broome.csv' codedir = getcwd()+sep codedirname = codedir + 'project.py' meshname = outputtimedir + 'mesh_' + basename # Necessary if using point datasets, rather than grid onshore_dem_name = datadir + onshore_name offshore_interp_dem_name = datadir + bathy_interp 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 coast_dem_name = datadir + coast_name combined_dem_name = datadir + 'broome_combined_elevation' ############################### # Domain definitions ############################### # 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') # 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 ################################################################### # exporting asc grid eastingmin = 412036.43 eastingmax = 427184.37 northingmin = 8007984 northingmax = 8029830 ############################### # Interior region definitions ############################### # broome digitized polygons poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon_update.csv') poly_broome2 = read_polygon(polygondir+'Broome_Close2_update.csv') poly_broome3 = read_polygon(polygondir+'Broome_Coast_update2.csv') #poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv') #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 #print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0 for i in poly_broome3: v = is_inside_polygon(i,poly_broome1, verbose=False) if v == False: print v def number_mesh_triangles(interior_regions, bounding_poly, remainder_res): from anuga.utilities.polygon import polygon_area # TO DO check if any of the regions fall inside one another no_triangles = 0.0 area = polygon_area(bounding_poly) for i,j in interior_regions: this_area = polygon_area(i) no_triangles += this_area/j area -= this_area print j, this_area/1000000., area/1000000. no_triangles += area/remainder_res return int(no_triangles/0.7)