Changeset 6261 for anuga_work


Ignore:
Timestamp:
Feb 3, 2009, 11:51:59 AM (15 years ago)
Author:
ole
Message:

Looped elevation data conversions

Location:
anuga_work/production/busselton/standardised_version
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/busselton/standardised_version/build_busselton.py

    r6254 r6261  
    1414
    1515# Standard modules
    16 from os import sep
    17 from os.path import dirname, basename
     16from os.path import dirname
    1817from os import mkdir, access, F_OK
    19 from shutil import copy
    20 import time
    21 import sys
    2218
    2319# Related major packages
    … …  
    3531#------------------------------------------------------------------------------
    3632
    37 copy_code_files(project.output_build_time_dir,__file__,
    38                dirname(project.__file__)+sep+ project.__name__+'.py' )
     33#copy_code_files(project.output_build_time_dir,__file__,
     34#               dirname(project.__file__)+os.sep+ project.__name__+'.py' )
     35#start_screen_catcher(project.output_build_time_dir)
    3936
    40 start_screen_catcher(project.output_build_time_dir)
    41 
    42 print 'USER:    ', project.user
    4337
    4438#-------------------------------------------------------------------------------
    … …  
    5044#-------------------------------------------------------------------------------
    5145
    52 # FIXME(Ole): Here's the Shark Bay way.
     46print 'project.poly_all', project.poly_all
     47print 'project.combined_dir_name', project.combined_dir_name
    5348
    54 ## FIXME (Ole): Clip data by interior regions if possible
     49# input elevation directory filenames
     50#onshore_in_dir_name = project.onshore_in_dir_name
     51#coast_in_dir_name = project.coast_in_dir_name
     52#coast_in_dir_name1 = project.coast_in_dir_name1
     53#offshore_in_dir_name = project.offshore_in_dir_name
     54#offshore_in_dir_name1 = project.offshore_in_dir_name1
     55#offshore_in_dir_name2 = project.offshore_in_dir_name2
     56#offshore_in_dir_name3 = project.offshore_in_dir_name3
     57#offshore_in_dir_name4 = project.offshore_in_dir_name4
     58#offshore_in_dir_name5 = project.offshore_in_dir_name5
     59
     60# output elevation directory filenames
     61#onshore_dir_name = project.onshore_dir_name
     62#coast_dir_name = project.coast_dir_name
     63#coast_dir_name1 = project.coast_dir_name1
     64#offshore_dir_name = project.offshore_dir_name
     65#offshore_dir_name1 = project.offshore_dir_name1
     66#offshore_dir_name2 = project.offshore_dir_name2
     67#offshore_dir_name3 = project.offshore_dir_name3
     68#offshore_dir_name4 = project.offshore_dir_name4
     69#offshore_dir_name5 = project.offshore_dir_name5
     70
     71
     72# Create Geospatial data from ASCII files
     73geospatial_data = {}
     74for filename in project.ascii_grid_filenames:
     75    absolute_filename = project.topographies_in_dir + filename
     76    convert_dem_from_ascii2netcdf(absolute_filename,
     77                                  basename_out=absolute_filename,
     78                                  use_cache=True,
     79                                  verbose=True)
     80    dem2pts(absolute_filename, use_cache=True, verbose=True)
     81
     82    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename+'.pts',
     83                                                verbose=True)
     84
     85# Create Geospatial data from TXT files
     86for filename in project.point_filenames:
     87    absolute_filename = project.topographies_in_dir + filename
     88    geospatial_data[filename] = Geospatial_data(file_name=absolute_filename,
     89                                                verbose=True)
     90
     91
     92#print "creates DEMs from asc data"
     93#print onshore_in_dir_name, onshore_dir_name
     94#print offshore_in_dir_name5, offshore_dir_name5
    5595#
    56 #print 'creating geospatial data objects from asc data (via dem and pts formats)'
    57 #for filename in project.ascii_grid_filenames:
    58 #    convert_dem_from_ascii2netcdf(filename,
    59 #                                  basename_out=filename,
    60 #                                  use_cache=True, verbose=True)
    61 #    dem2pts(filename, use_cache=True, verbose=True)#
     96#convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
    6297#
    63 #    geospatial_data += Geospatial_data(file_name=filename + '.pts',
    64 #                                       verbose=True)
     98#convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)#
    6599##
    66100#
    67 #print 'creating geospatial data objects from txt data'
    68 #for filename in project.point_filenames:
    69 #    geospatial_data += Geospatial_data(file_name=filename + '.txt',
    70 #                                       verbose=True)
    71 #
    72 #
    73 #print 'clip combined geospatial object by bounding polygon'
    74 #G = geospatial_data.clip(project.bounding_polygon)
    75 #
    76 #
    77 #print 'export combined geospatial data'
    78 #if access(project.topographies_dir, F_OK) == 0:
    79 #    mkdir (project.topographies_dir)
    80 #G.export_points_file(project.combined_dir_name + '.pts')
    81 
    82 
    83 
    84 print "project.poly_all", project.poly_all
    85 print "project.combined_dir_name", project.combined_dir_name
    86 
    87 # input elevation directory filenames
    88 onshore_in_dir_name = project.onshore_in_dir_name
    89 coast_in_dir_name = project.coast_in_dir_name
    90 coast_in_dir_name1 = project.coast_in_dir_name1
    91 offshore_in_dir_name = project.offshore_in_dir_name
    92 offshore_in_dir_name1 = project.offshore_in_dir_name1
    93 offshore_in_dir_name2 = project.offshore_in_dir_name2
    94 offshore_in_dir_name3 = project.offshore_in_dir_name3
    95 offshore_in_dir_name4 = project.offshore_in_dir_name4
    96 offshore_in_dir_name5 = project.offshore_in_dir_name5
    97 
    98 # output elevation directory filenames
    99 onshore_dir_name = project.onshore_dir_name
    100 coast_dir_name = project.coast_dir_name
    101 coast_dir_name1 = project.coast_dir_name1
    102 offshore_dir_name = project.offshore_dir_name
    103 offshore_dir_name1 = project.offshore_dir_name1
    104 offshore_dir_name2 = project.offshore_dir_name2
    105 offshore_dir_name3 = project.offshore_dir_name3
    106 offshore_dir_name4 = project.offshore_dir_name4
    107 offshore_dir_name5 = project.offshore_dir_name5
    108 # creates DEM from asc data
    109 print "creates DEMs from asc data"
    110 convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
    111 convert_dem_from_ascii2netcdf(offshore_in_dir_name5, basename_out=offshore_dir_name5, use_cache=True, verbose=True)
    112 
    113101# creates pts file for onshore DEM
    114 print "creates pts file for onshore DEM"
    115 dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
    116 dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
     102#print "creates pts file for onshore DEM"
     103#dem2pts(onshore_dir_name ,use_cache=True,verbose=True)
     104#dem2pts(offshore_dir_name5 ,use_cache=True,verbose=True)
    117105
    118106# create onshore pts files
    119 print'create Geospatial data1 objects from topographies'
    120 G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
     107#print'create Geospatial data1 objects from topographies'
     108#G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
    121109
    122110# create coastal and offshore pts files
    123 print'create Geospatial data2 objects from topographies'
    124 G2 = Geospatial_data(file_name = coast_in_dir_name)
    125 print'create Geospatial data3 objects from topographies'
    126 G3 = Geospatial_data(file_name = coast_in_dir_name1)
    127 print'create Geospatial data4 objects from topographies'
    128 G_off = Geospatial_data(file_name = offshore_in_dir_name)
    129 print'create Geospatial data5 objects from topographies'
    130 G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
    131 print'create Geospatial data6 objects from topographies'
    132 G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
    133 print'create Geospatial data7 objects from topographies'
    134 G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
    135 print'create Geospatial data8 objects from topographies'
    136 G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
    137 print'create Geospatial data9 objects from topographies'
    138 G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
     111#print'create Geospatial data2 objects from topographies'
     112#G2 = Geospatial_data(file_name = coast_in_dir_name)
     113#print'create Geospatial data3 objects from topographies'
     114#G3 = Geospatial_data(file_name = coast_in_dir_name1)
     115#print'create Geospatial data4 objects from topographies'
     116#G_off = Geospatial_data(file_name = offshore_in_dir_name)
     117#print'create Geospatial data5 objects from topographies'
     118#G_off1 = Geospatial_data(file_name = offshore_in_dir_name1)
     119#print'create Geospatial data6 objects from topographies'
     120#G_off2 = Geospatial_data(file_name = offshore_in_dir_name2)
     121#print'create Geospatial data7 objects from topographies'
     122#G_off3 = Geospatial_data(file_name = offshore_in_dir_name3)
     123#print'create Geospatial data8 objects from topographies'
     124#G_off4 = Geospatial_data(file_name = offshore_in_dir_name4)
     125#print'create Geospatial data9 objects from topographies'
     126#G_off5 = Geospatial_data(file_name = offshore_dir_name5 + '.pts')
    139127
    140128
    … …  
    143131#-------------------------------------------------------------------------------
    144132
    145 print'add all geospatial objects'
    146 G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
     133#print'add all geospatial objects'
     134#G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
    147135
    148 print'clip combined geospatial object by bounding polygon'
     136print 'Add geospatial objects except', project.offshore_name5
     137print geospatial_data.keys()
     138G = None
     139for key in geospatial_data:
     140    if key != project.offshore_name5:
     141        G += geospatial_data[key]
     142       
     143print 'Clip (outside) combined geospatial object except %s by area of interest.'
    149144G_clip = G.clip_outside(project.poly_aoi1)
    150 G_all = G_clip + G_off5
     145
     146print 'Clip combined geospatial object by bounding polygon'
     147G_all = G_clip + geospatial_data[project.offshore_name5]   # G_off5
    151148G_clipped = G_all.clip(project.poly_all)
    152149
    153 print'export combined DEM file'
    154 if access(project.topographies_dir,F_OK) == 0:
    155     mkdir (project.topographies_dir)
     150print 'Export combined DEM file'
     151#if access(project.topographies_dir,F_OK) == 0:
     152#    mkdir (project.topographies_dir)
    156153G_clipped.export_points_file(project.combined_dir_name + '.pts')
    157 G_clipped.export_points_file(project.combined_dir_name + '.txt') #Use for comparision in ARC
     154print 'Do txt version too'
     155G_clipped.export_points_file(project.combined_dir_name + '.txt') # Use for comparision in ARC
    158156

  • anuga_work/production/busselton/standardised_version/project.py

    r6256 r6261  
    66#------------------------------------------------------------------------------
    77
     8import os
    89from os import sep, environ, getenv, getcwd
    910from os.path import expanduser
    … …  
    5253finaltime=80000         # final time for simulation
    5354
    54 setup='final'  # Final can be replaced with trial or basic.
     55setup='trial'  # Final can be replaced with trial or basic.
    5556               # Either will result in a coarser mesh that will allow a
    5657               # faster, but less accurate, simulation.
    … …  
    5859if setup =='trial':
    5960    print'trial'
    60     res_factor=10
    61     time_thinning=48
     61    res_factor=100
     62    time_thinning=96
    6263    yieldstep=240
    6364if setup =='basic':
    … …  
    102103#------------------------------------------------------------------------------
    103104
    104 # FIXME(Ole): Use lists of elevation data files as in Shark Bay study:
    105 ## elevation data filenames
    106 #ascii_grid_filenames = ['10m_dem_without_survey', '50m_dem_without_10m_dem',
    107 #                        'bathysteeppt', 'bathyleft', 'bathyright']
    108 #point_filenames = ['field_survey_north', 'field_survey_south',
    109 #                   'clipped_bathymetry_final', 'coast_points_final']
    110 
    111 
    112105
    113106# elevation data used in build_busselton.py
    … …  
    140133landward = 'landward_bounding_polygon.csv'
    141134
     135
    142136#------------------------------------------------------------------------------
    143137# Output Elevation Data
    … …  
    164158# Location of input and output data
    165159#------------------------------------------------------------------------------
    166 # where the input data sits
    167 onshore_in_dir_name = topographies_in_dir + onshore_name #topo
    168 coast_in_dir_name = topographies_in_dir + coast_name #coastline
    169 coast_in_dir_name1 = topographies_in_dir + coast_name1 #beach survey
    170 offshore_in_dir_name = topographies_in_dir + offshore_name #bathymetry
    171 offshore_in_dir_name1 = topographies_in_dir + offshore_name1 #bathymetry Charts
    172 offshore_in_dir_name2 = topographies_in_dir + offshore_name2 #Digitised Fairsheet
    173 offshore_in_dir_name3 = topographies_in_dir + offshore_name3 #250m
    174 offshore_in_dir_name4 = topographies_in_dir + offshore_name4 #Bunbury DPI
    175 offshore_in_dir_name5 = topographies_in_dir + offshore_name5 #Busselton Topo
    176 
    177 # where the output data sits
    178 onshore_dir_name = topographies_dir + onshore_name
    179 
    180 coast_dir_name = topographies_dir + coast_name
    181 coast_dir_name1 = topographies_dir + coast_name1
    182 
    183 offshore_dir_name = topographies_dir + offshore_name
    184 offshore_dir_name1 = topographies_dir + offshore_name1
    185 offshore_dir_name2 = topographies_dir + offshore_name2
    186 offshore_dir_name3 = topographies_dir + offshore_name3
    187 offshore_dir_name4 = topographies_dir + offshore_name4
    188 offshore_dir_name5 = topographies_dir + offshore_name5
    189 
    190 # where the combined elevation file sits
     160
     161ascii_grid_filenames = [onshore_name,   # Topo
     162                        offshore_name5] # Busselton Topo
     163
     164point_filenames = [coast_name,     # Coastline
     165                   coast_name1,    # Beach survey
     166                   offshore_name,  # Bathymetry
     167                   offshore_name1, # Bathymetry Charts
     168                   offshore_name2, # Digitised Fairsheet
     169                   offshore_name3, # 250m
     170                   offshore_name4] # Bunbury DPI
     171
     172
     173# Where the combined elevation file sits
    191174combined_dir_name = topographies_dir + combined_name
    192175combined_smaller_name_dir = topographies_dir + combined_smaller_name
    193176
    194 # where the mesh sits (this is created during the run_busselton.py)
     177# Where the mesh sits (this is created during the run_busselton.py)
    195178meshes_dir_name = meshes_dir + scenario_name+'.msh'
    196179
    197 # where the boundary ordering files sit (this is used within build_boundary.py)
     180# Where the boundary ordering files sit (this is used within build_boundary.py)
    198181order_filename_dir = boundaries_dir + order_filename
    199182
    200 # where the landward points of boundary extent sit (this is used within run_busselton.py)
     183# Where the landward points of boundary extent sit (this is used within run_busselton.py)
    201184landward_dir = boundaries_dir + landward
    202185
    203 # where the event sts files sits (this is created during the build_boundary.py)
     186# Where the event sts files sits (this is created during the build_boundary.py)
    204187boundaries_dir_event = boundaries_dir + str(event_number) + sep
    205188boundaries_dir_mux = muxhome
    206 
    207 # where the directory of the output filename sits
     189urs_boundary_name = os.path.join(boundaries_dir_event,
     190                                 scenario_name)
     191
     192
     193# Where the directory of the output filename sits
    208194output_build_time_dir = output_dir+build_time+dir_comment+sep   #used for build_busselton.py
    209195output_run_time_dir = output_dir+run_time+dir_comment+sep       #used for run_busselton.py
    210196output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
    211197
    212 #w here the directory of the gauges sit
     198# Where the directory of the gauges sit
    213199gauges_dir_name = gauges_dir + gauge_name       #used for get_timeseries.py
    214200building_in_dir_name = gauges_dir + building + '.csv'    #used for run_building_inundation.py
    … …  
    218204#------------------------------------------------------------------------------
    219205
    220 #Land, to set the initial stage/water to be offcoast only
     206# Land, to set the initial stage/water to be offcoast only
    221207poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
    222208
    223 #Land, to set the initial stage/water to be offcoast only
     209# Land, to set the initial stage/water to be offcoast only
    224210poly_marina = read_polygon(polygons_dir+'initial_condition_marina.csv')
    225211
    … …  
    262248
    263249   
    264 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
     250trigs_min = number_mesh_triangles(interior_regions, poly_all,
     251                                  res_poly_all)
    265252print 'min estimated number of triangles', trigs_min
    266253   
    … …  
    270257#------------------------------------------------------------------------------
    271258
    272 # exporting asc grid for Busselton
     259# ASCII export grid for Busselton
    273260xminBusselton = 340000
    274261xmaxBusselton = 352000
    … …  
    276263ymaxBusselton = 6280000
    277264
    278 # exporting asc grid for Bunbury
     265# ASCII export grid for Bunbury
    279266xminBunbury = 369000
    280267xmaxBunbury = 381000

  • anuga_work/production/busselton/standardised_version/run_busselton.py

    r6259 r6261  
    2525
    2626# Standard modules
    27 from os import sep
    2827import os
    29 from os.path import dirname, basename
    30 from os import mkdir, access, F_OK
    31 from shutil import copy
     28from os.path import dirname
    3229import time
    33 import sys
    3430
    3531# Related major packages
    … …  
    4844import project  # Definition of file names and polygons
    4945
    50    
     46
    5147#-----------------------------------------------------------------------
    5248# Copy scripts to time stamped output directory and capture screen
    … …  
    5652#copy script must be before screen_catcher
    5753#copy_code_files(project.output_run_time_dir, __file__,
    58 #         dirname(project.__file__)+sep+ project.__name__+'.py' )
     54#         dirname(project.__file__)+os.sep+ project.__name__+'.py' )
    5955#start_screen_catcher(project.output_run_time_dir, myid, numprocs)
    6056
    … …  
    9591print domain.statistics()
    9692
    97 
    9893domain.set_name(project.scenario_name)
    9994domain.set_datadir(project.output_run_time_dir)
    … …  
    111106                      geo_reference=domain.geo_reference)
    112107domain.set_quantity('stage', IC, use_cache=True, verbose=True)
    113 
    114108domain.set_quantity('friction', project.friction)
    115 
    116109domain.set_quantity('elevation',
    117110                    filename=project.combined_dir_name+'.pts',
    … …  
    124117# Setup boundary conditions
    125118#-------------------------------------------------------------------------
    126 
    127119print 'Set boundary - available tags:', domain.get_boundary_tags()
    128 
    129 boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
    130120
    131121Br = Reflective_boundary(domain)
    132122Bd = Dirichlet_boundary([project.tide,0,0])
    133 
    134 
    135 Bf = Field_boundary(boundary_urs_out+'.sts',
     123Bf = Field_boundary(project.urs_boundary_name+'.sts',
    136124                    domain, mean_stage=project.tide,
    137125                    time_thinning=1,
    … …  
    140128                    use_cache=True,
    141129                    verbose=True)
    142 
    143130
    144131domain.set_boundary({'back': Br,
Note: See TracChangeset for help on using the changeset viewer.