Changeset 4077


Ignore:
Timestamp:
Dec 14, 2006, 11:41:51 AM (17 years ago)
Author:
nick
Message:

start perth scenario

Location:
anuga_work/production/perth_2006
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/perth_2006/project.py

    r3908 r4077  
    1010#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
    1111
    12 if sys.platform == 'win32':
    13     home = getenv('INUNDATIONHOME')
    14     user = getenv('USERPROFILE')
     12# file and system info
     13#---------------------------------
     14codename = 'project.py'
    1515
    16 else:   
    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
     16home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
     17user = get_user_name()
    2018
    2119# INUNDATIONHOME is the inundation directory, not the data directory.
     
    2422#Making assumptions about the location of scenario data
    2523state = 'western_australia'
     24scenario_name = 'perth'
    2625scenario_dir_name = 'perth_tsunami_scenario_2006'
    2726
    2827# onshore data provided by WA DLI
    29 onshore_name = '' # original
     28onshore_name = 'perth_dli_ext' # original
     29#island
     30island_name = 'rott_dli_ext' # original
     31island_name1 = 'gard_dli_ext'
     32island_name2 = 'carnac_island_dted'
    3033
    3134# AHO + DPI data
    32 offshore_name1 = 'XY100011610'
    33 offshore_name2 = 'XY100011611'
    34 offshore_name3 = 'XY100011613'
    35 offshore_name4 = 'XY100011614'
    36 offshore_name5 = 'XY100011616'
    37 offshore_name6 = 'XY100011617'
    38 offshore_name7 = 'XY100011618'
    39 offshore_name8 = 'XY100011621'
    40 offshore_name9 = 'XY100011623'
    41 offshore_name10 = 'XY100011745'
    42 offshore_name11 = 'XY100011746'
    43 offshore_name12 = 'XY100017530'
    44 offshore_name13 = 'XY100017532'
    45 offshore_name14 = 'XY100017538'
    46 offshore_name15 = 'XY100017540'
    47 offshore_name16 = 'XYBR66'
    48 offshore_name17 = 'XYBR70'
    49 offshore_name18 = 'XYBR80'
    50 offshore_name19 = 'XYBR88'
    51 offshore_name20 = 'XYBR93'
    52 offshore_name21 = 'XYBR0110'
    53 offshore_name22 = 'XYWADPI'
     35coast_name = 'waterline'
     36offshore_name = 'perth_bathymetry'
    5437
    55 # developed by NM&I
    56 coast_name = 'coastline'
     38#final topo name
     39combined_name ='perth_combined_elevation'
    5740
    58 boundary_basename = 'SU-AU' # Mw ?
     41topographies_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'topographies'+sep
     42topographies_time_dir = topographies_dir+build_time+sep
    5943
    60 #swollen/ all data output
    61 basename = 'source'
    62 codename = 'project.py'
     44#input topo file location
     45onshore_dir_name = topographies_dir + onshore_name
     46island_dir_name = topographies_dir + island_name
     47island_dir_name1 = topographies_dir + island_name1
     48island_dir_name2 = topographies_dir + island_name2
     49
     50coast_dir_name = topographies_dir + coast_name
     51offshore_dir_name = topographies_dir + offshore_name
     52
     53#final topo files
     54combined_dir_name = topographies_dir + combined_name
     55combined_time_dir_name = topographies_time_dir + combined_name
     56combined_time_dir_final_name = topographies_time_dir + combined_final_name
     57
    6358
    6459#Derive subdirectories and filenames
    65 local_time = strftime('%Y%m%d_%H%M%S',gmtime())
    66 meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
    67 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
    68 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
    69 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
    70 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
    71 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
    72 outputtimedir = outputdir + local_time + sep
    73 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
    7460
    75 gauge_filename = gaugedir + 'perth_gauges.csv'
    76 community_filename = gaugedir + 'CHINS_v2.csv'
    77 community_perth = gaugedir + 'community_perth.csv'
    78 codedir = getcwd()+sep                           
    79 codedirname = codedir + 'project.py'
    80 meshname = outputtimedir + 'mesh_' + basename
     61meshes_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'meshes'+sep
     62meshes_dir_name = meshes_dir + scenario_name
    8163
    82 # Necessary if using point datasets, rather than grid
    83 onshore_dem_name = datadir + onshore_name
    84 offshore_dem_name1 = datadir + offshore_name1
    85 offshore_dem_name2 = datadir + offshore_name2
    86 offshore_dem_name3 = datadir + offshore_name3
    87 offshore_dem_name4 = datadir + offshore_name4
    88 offshore_dem_name5 = datadir + offshore_name5
    89 offshore_dem_name6 = datadir + offshore_name6
    90 offshore_dem_name7 = datadir + offshore_name7
    91 offshore_dem_name8 = datadir + offshore_name8
    92 offshore_dem_name9 = datadir + offshore_name9
    93 offshore_dem_name10 = datadir + offshore_name10
    94 offshore_dem_name11 = datadir + offshore_name11
    95 offshore_dem_name12 = datadir + offshore_name12
    96 offshore_dem_name13 = datadir + offshore_name13
    97 offshore_dem_name14 = datadir + offshore_name14
    98 offshore_dem_name15 = datadir + offshore_name15
    99 offshore_dem_name16 = datadir + offshore_name16
    100 offshore_dem_name17 = datadir + offshore_name17
    101 offshore_dem_name18 = datadir + offshore_name18
    102 offshore_dem_name19 = datadir + offshore_name19
    103 offshore_dem_name20 = datadir + offshore_name20
    104 offshore_dem_name21 = datadir + offshore_name21
    105 offshore_dem_name22 = datadir + offshore_name22
     64polygons_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'polygons'+sep
     65tide_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'tide_data'+sep
    10666
    107 coast_dem_name = datadir + coast_name
     67#boundaries locations
     68boundaries_in_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
     69boundaries_in_dir_name = boundaries_in_dir + boundaries_name
     70boundaries_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep
     71boundaries_dir_name = boundaries_dir + boundaries_name
     72#boundaries_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
     73#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
     74#ideas
     75#boundaries_time_dir = boundaries_in_dir+'urs'+sep+boundaries_source+sep
    10876
    109 combined_dem_name   = datadir + 'perth_combined_elevation'
     77#time stuff
     78time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     79gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
     80build_time = time+'_build'
     81run_time = time+'_run'
     82print 'gtime: ', gtime
     83
     84#output locations
     85output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
     86output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep
     87output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep
     88output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
     89
     90#gauges
     91gauge_name = 'perth.csv'
     92gauges_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'gauges'+sep
     93gauges_dir_name = gauges_dir + gauge_name
     94
     95
     96
     97
     98
     99
     100
    110101
    111102###############################
    112103# Domain definitions
    113104###############################
     105from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    114106
    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
     107bounding_polygon = read_polygon(polygons_dir+'bounding_poly.csv')
    127108
    128 # bounding polygon for study area
    129 polyAll = read_polygon(polygondir+'extent.csv')
     109refzone = 50
    130110
    131 # plot bounding polygon and make sure BC info surrounds it
    132 #plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
    133 print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
     111print "bounding_polygon", bounding_polygon
     112print 'poly area', polygon_area(bounding_polygon)/1000000.0
     113
     114
     115#Interior regions
     116
     117#cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
     118#assert zone == refzone
     119
     120poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20pts.csv')
     121
     122poly_cbd = read_polygon(polygons_dir+'cbd_pts.csv')
     123
     124poly_penguin = read_polygon(polygons_dir+'penguin_pts.csv')
     125assert zone == refzone
     126
    134127
    135128###################################################################
    136129# Clipping regions for export to asc and regions for clipping data
    137130###################################################################
    138 
     131'''
    139132# exporting asc grid
    140133eastingmin = 406215.87
     
    158151print 'Area of coastal polygon', polygon_area(poly_perth3)/1000000.0
    159152#print 'Area of cable beach polygon', polygon_area(poly_perth4)/1000000.0
     153'''
    160154
    161 for i in poly_perth3:
    162     v = is_inside_polygon(i,poly_perth1, verbose=False)
    163     if v == False: print v
    164 
    165 def 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 TracChangeset for help on using the changeset viewer.