Ignore:
Timestamp:
Sep 24, 2014, 12:41:57 PM (11 years ago)
Author:
davies
Message:

Updates to quantity setting functions + spatialInputUtil
(and various template script modifications)

Location:
trunk/anuga_work/development/gareth/template_scenarios/tsunami
Files:
2 added
19 deleted
2 edited
8 moved

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/run_model.py

    r9341 r9344  
    3636# Routines to setup the domain
    3737
    38 import setup_boundary_conditions
    39 import setup_rainfall
    40 import setup_inlets
    41 import setup_mesh
    42 import setup_initial_conditions
     38from setup import setup_boundary_conditions
     39from setup import setup_rainfall
     40from setup import setup_inlets
     41from setup import setup_mesh
     42from setup import setup_initial_conditions
     43from setup import raster_outputs
     44from setup.prepare_data import PrepareData
    4345
    4446# Routines defined by the user
    … …  
    4850# Functions to make geotifs
    4951
    50 import raster_outputs
    5152
    5253# Record the time so we can report how long the simulation takes
    … …  
    5657# Get key data for the simulation
    5758
    58 from prepare_data import PrepareData
    5959project = PrepareData('ANUGA_setup.xls')
    6060

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup/parse_input_data.py

    r9343 r9344  
    77
    88"""
     9
     10
    911import os
    1012import xlrd
     13import glob
     14
    1115
    1216class ProjectData(object):
    … …  
    2529        file_extension = os.path.splitext(filename)[1]
    2630
    27         if(file_extension in ['.xls', '.xlsx']):
     31        if file_extension in ['.xls', '.xlsx']:
    2832            self.get_data_from_excel(filename)
    2933        else:
    3034            msg = 'File type not supported'
    3135            raise ValueError(msg)
    32    
    3336
    3437    def get_data_from_excel(self, filename):
    … …  
    4346
    4447        # Worksheet names -- define here so they can easily be changed
     48
    4549        project_ws = 'project_settings'
    4650        mesh_ws = 'mesh'
    … …  
    5761        # #####################################################################
    5862
    59         self.scenario = data_source.get_var(
    60             project_ws, 'scenario', [0, 1], post_process = str)
     63        self.scenario = data_source.get_var(project_ws, 'scenario', [0, 1],
     64                                            post_process=str)
    6165
    6266        self.output_basedir = data_source.get_var(
    63             project_ws, 'output_base_directory', [0, 1], post_process= str)
     67            project_ws, 'output_base_directory', [
     68                0, 1], post_process=str)
    6469
    6570        self.yieldstep = data_source.get_var(project_ws, 'yieldstep', [0, 1],
    66                                            post_process=float)
     71                                             post_process=float)
    6772
    6873        self.finaltime = data_source.get_var(project_ws, 'finaltime', [0, 1],
    69                                            post_process=float)
     74                                             post_process=float)
    7075
    7176        self.projection_information = data_source.get_var(
    72             project_ws, 'projection_information', [0, 1])
     77            project_ws, 'projection_information', [
     78                0, 1])
    7379
    7480        # Could be unicode string or float -- if float, we need int
    … …  
    7682        if isinstance(self.projection_information, float):
    7783            self.projection_information = int(self.projection_information)
    78         if type(self.projection_information) is unicode:
     84        if isinstance(self.projection_information, unicode):
    7985            self.projection_information = str(self.projection_information)
    8086
    81         self.flow_algorithm = data_source.get_var(project_ws, 'flow_algorithm',
    82                                                 [0, 1], post_process = str)
     87        self.flow_algorithm = data_source.get_var(
     88            project_ws, 'flow_algorithm', [
     89                0, 1], post_process=str)
    8390
    8491        # Coerce this to a logical variable
    … …  
    8693        self.use_local_extrapolation_and_flux_updating = \
    8794            data_source.get_var(project_ws,
    88                               'use_local_extrapolation_and_flux_updating',
    89                               offset=[0, 1], post_process=bool)
     95                                'use_local_extrapolation_and_flux_updating',
     96                                offset=[0, 1], post_process=bool)
    9097
    9198        self.output_tif_cellsize = data_source.get_var(
    9299            project_ws,
    93100            'output_tif_cellsize',
    94             offset=[0, 1],
     101            offset=[
     102                0,
     103                1],
    95104            post_process=float)
    96105
    … …  
    98107            project_ws,
    99108            'output_tif_bounding_polygon',
    100             offset=[0, 1],
     109            offset=[
     110                0,
     111                1],
    101112            post_process=empty_string_to_none)
    102113
    … …  
    104115            project_ws,
    105116            'max_quantity_update_frequency',
    106             offset=[0, 1],
     117            offset=[
     118                0,
     119                1],
    107120            post_process=int)
    108121
    … …  
    110123            project_ws,
    111124            'max_quantity_collection_starttime',
    112             offset=[0, 1],
    113             post_process = float)
     125            offset=[
     126                0,
     127                1],
     128            post_process=float)
    114129
    115130        self.store_vertices_uniquely = data_source.get_var(
    116131            project_ws,
    117132            'store_vertices_uniquely',
    118             offset=[0, 1],
     133            offset=[
     134                0,
     135                1],
    119136            post_process=bool)
    120137
    … …  
    122139            project_ws,
    123140            'spatial_text_output_dir',
    124             offset = [0, 1],
     141            offset=[
     142                0,
     143                1],
    125144            post_process=str)
    126145
    … …  
    132151
    133152        self.bounding_polygon_and_tags_file = data_source.get_var(
    134             mesh_ws, 'bounding_polygon', [1, 1], post_process=str)
     153            mesh_ws, 'bounding_polygon', [
     154                1, 1], post_process=str)
    135155
    136156        self.default_res = data_source.get_var(mesh_ws, 'bounding_polygon',
    137                                              [2, 1], post_process = str)
    138 
    139         # Loop to read the interior regions data
    140         # Finish when we hit an empty cell
     157                                               [2, 1], post_process=str)
    141158
    142159        self.interior_regions_data = data_source.get_paired_list(
    143             mesh_ws, 'interior_regions', [1, 1],
    144             post_process = string_or_float)
     160            mesh_ws, 'interior_regions', [
     161                1, 1], post_process=string_or_float)
    145162
    146163        # breaklines, riverwalls, point movement threshold
    147 
    148         self.breakline_files = data_source.get_list(mesh_ws, 'breakline_files',
    149                                                   [0, 1], post_process = str)
    150 
    151         self.riverwall_csv_files = data_source.get_list(
    152             mesh_ws, 'riverwall_csv_files', [0, 1], post_process = str)
     164        # We allow wildcard type names as input files
     165
     166        breakline_files = data_source.get_list(
     167            mesh_ws, 'breakline_files', [
     168                0, 1], post_process=str)
     169
     170        riverwall_csv_files = data_source.get_list(
     171            mesh_ws, 'riverwall_csv_files', [
     172                0, 1], post_process=str)
     173
     174        # Allow wildcard names as input for breaklines
     175
     176        self.breakline_files = []
     177        if len(breakline_files) > 0:
     178            self.breakline_files = [self.breakline_files + glob.glob(bl_file)
     179                                    for bl_file in breakline_files]
     180
     181        # Allow wildcards for breaklines
     182
     183        self.riverwall_csv_files = []
     184        if len(riverwall_csv_files) > 0:
     185            self.riverwall_csv_files = [
     186                self.riverwall_csv_files +
     187                glob.glob(rw_file) for rw_file in riverwall_csv_files]
    153188
    154189        self.break_line_intersect_point_movement_threshold = \
    155190            data_source.get_var(mesh_ws, 'breakline_intersection_threshold',
    156                               [0, 1], post_process = float)
     191                                [0, 1], post_process=float)
    157192
    158193        self.pt_areas = data_source.get_var(
    159             mesh_ws, 'region_areas_file', [0, 1], empty_string_to_none)
     194            mesh_ws, 'region_areas_file', [
     195                0, 1], empty_string_to_none)
    160196
    161197        # Check we are only using interior_region OR breakline methods
     198
    162199        use_ir = self.interior_regions_data != []
    163         use_bl = (
    164                   (self.breakline_files != []) or\
    165                   ((self.riverwall_csv_files != []) or\
    166                    (self.pt_areas is not None))
    167                  )
    168 
    169         if (use_ir and use_bl):
    170             msg = ' In worksheet ' + mesh_ws +\
    171                  'Cannot have both *interior_regions_data* and *breaklines ' \
     200        use_bl = self.breakline_files != [] or self.riverwall_csv_files != [] \
     201            or self.pt_areas is not None
     202
     203        if use_ir and use_bl:
     204            msg = ' In worksheet ' + mesh_ws \
     205                + 'Cannot have both *interior_regions_data* and *breaklines ' \
    172206                + '/ riverwalls / region point areas* being non-empty.'
    173207            raise ValueError(msg)
    … …  
    180214
    181215        self.elevation_data = data_source.get_paired_list(
    182             init_ws, 'Elevation', [1, 1], post_process = string_or_float)
     216            init_ws, 'Elevation', [
     217                1, 1], post_process=string_or_float)
    183218
    184219        self.friction_data = data_source.get_paired_list(
    185             init_ws, 'Friction', [1, 1], post_process = string_or_float)
     220            init_ws, 'Friction', [
     221                1, 1], post_process=string_or_float)
    186222
    187223        self.stage_data = data_source.get_paired_list(
    188             init_ws, 'Stage', [1, 1], post_process = string_or_float)
     224            init_ws, 'Stage', [
     225                1, 1], post_process=string_or_float)
    189226
    190227        self.xmomentum_data = data_source.get_paired_list(
    191             init_ws, 'X-velocity_times_depth', [1, 1],
    192             post_process = string_or_float)
     228            init_ws, 'X-velocity_times_depth',
     229            [1, 1], post_process=string_or_float)
    193230
    194231        self.ymomentum_data = data_source.get_paired_list(
    195             init_ws, 'Y-velocity_times_depth', [1, 1], 
    196             post_process = string_or_float)
     232            init_ws, 'Y-velocity_times_depth', [1, 1],
     233            post_process=string_or_float)
    197234
    198235        self.elevation_additions = data_source.get_paired_list(
    199             init_add_ws, 'Elevation', [1, 1],
    200             post_process = string_or_float)
     236            init_add_ws, 'Elevation', [
     237                1, 1], post_process=string_or_float)
    201238
    202239        self.friction_additions = data_source.get_paired_list(
    203             init_add_ws, 'Friction', [1, 1],
    204             post_process = string_or_float)
     240            init_add_ws, 'Friction', [
     241                1, 1], post_process=string_or_float)
    205242
    206243        self.stage_additions = data_source.get_paired_list(
    207             init_add_ws, 'Stage', [1, 1],
    208             post_process = string_or_float)
     244            init_add_ws, 'Stage', [
     245                1, 1], post_process=string_or_float)
    209246
    210247        self.xmomentum_additions = data_source.get_paired_list(
    211248            init_add_ws, 'X-velocity_times_depth', [1, 1],
    212             post_process = string_or_float)
     249            post_process=string_or_float)
    213250
    214251        self.ymomentum_additions = data_source.get_paired_list(
    215252            init_add_ws, 'Y-velocity_times_depth', [1, 1],
    216             post_process = string_or_float)
     253            post_process=string_or_float)
    217254
    218255        # #####################################################################
    … …  
    223260
    224261        # This can fail if we don't coerce to str
     262
    225263        self.boundary_tags_attribute_name = data_source.get_var(
    226             boundary_ws, 'boundary_tags_attribute_name', [0, 1], str)
     264            boundary_ws, 'boundary_tags_attribute_name', [
     265                0, 1], str)
    227266
    228267        self.boundary_data = data_source.get_subtable(
    229             boundary_ws, 'Boundary_Condition_Tag', [1, 0],
    230             post_process = string_or_float)
    231 
     268            boundary_ws, 'Boundary_Condition_Tag', [
     269                1, 0], post_process=string_or_float)
    232270
    233271        # #####################################################################
    … …  
    237275        # #####################################################################
    238276
    239         self.inlet_data = data_source.get_subtable(
    240             inlet_ws, 'Inlet_name', [1, 0])
     277        self.inlet_data = data_source.get_subtable(inlet_ws, 'Inlet_name', [1,
     278                                                                            0])
    241279
    242280        # #####################################################################
    … …  
    246284        # #####################################################################
    247285
    248         self.rain_data = data_source.get_subtable(
    249             rain_ws, 'Rainfall', [1, 1])
    250 
    251         # #####################################################################
    252 
    253 
    254 
    255        
    256 
    257 # #############################################################################
     286        self.rain_data = data_source.get_subtable(rain_ws, 'Rainfall', [1, 1])
     287
     288        # #####################################################################
     289
     290##############################################################################
    258291#
    259292# END OF CLASS
    260293#
    261 # #############################################################################
     294##############################################################################
     295
    262296
    263297class AnugaXls(object):
    … …  
    307341
    308342        # Loop through the rows and look for a match in the first column
     343
    309344        for i in range(len(sheet_data)):
    310345            if flag in sheet_data[i][0]:
    … …  
    320355
    321356        if matching_row + offset[0] >= len(sheet_data):
    322             return ""
     357            return ''
    323358
    324359        # If the col offset exceeds the data length, return a blank
     360
    325361        if offset[1] < len(sheet_data[i]):
    326362            output = sheet_data[i + offset[0]][offset[1]]
    … …  
    338374        flag,
    339375        offset=[0, 0],
    340         post_process = None,
     376        post_process=None,
    341377    ):
    342378        """Find the cell in the first column of worksheet_name that (partially)
    … …  
    366402            else:
    367403                cell_x = offset[0] + 1
    368                 next_row2 = self.get_var(worksheet_name, flag,
    369                                          [cell_x, cell_y],
    370                                         post_process)
     404                next_row2 = self.get_var(
     405                    worksheet_name, flag, [
     406                        cell_x, cell_y], post_process)
    371407
    372408            paired_list.append([next_row1, next_row2])
    … …  
    398434            cell_x = offset[0]
    399435            cell_y = counter + offset[1]
    400             next_row1 = self.get_var(worksheet_name, flag,
    401                                      [cell_x, cell_y],
     436            next_row1 = self.get_var(worksheet_name, flag, [cell_x, cell_y],
    402437                                     post_process)
    403438            if next_row1 == '':
    … …  
    415450        flag,
    416451        offset=[0, 0],
    417         post_process=None
     452        post_process=None,
    418453    ):
    419454        """Find the cell in the first column of worksheet_name that (partially)
    420            matches 'flag', and read data from a nearby 'sub-table' in the
    421            spreadsheet as a list of lists.
    422 
    423            The 'sub-table' starts at a location
    424            offset fom the matching cell by offset = [nrow, ncol],
    425            and finishes when a blank row of the sub-table is found
    426 
    427            Rows can have different numbers of non-empty columns, but each row
    428            ends when its first empty cell is found
     455        matches 'flag', and read data from a nearby 'sub-table' in the
     456        spreadsheet as a list of lists.
     457
     458        The 'sub-table' starts at a location
     459        offset fom the matching cell by offset = [nrow, ncol],
     460        and finishes when a blank row of the sub-table is found
     461
     462        Rows can have different numbers of non-empty columns, but each row
     463        ends when its first empty cell is found
    429464        """
    430        
    431         out_list = [ ]
     465
     466        out_list = []
    432467
    433468        counter = -1
    434         while True: 
     469        while True:
    435470            counter = counter + 1
    436471            cell_x = offset[0] + counter
    437472            cell_y = offset[1]
    438             next_row = self.get_list(worksheet_name, flag,
    439                                      [cell_x, cell_y],
     473            next_row = self.get_list(worksheet_name, flag, [cell_x, cell_y],
    440474                                     post_process)
    441475            if next_row == []:
    … …  
    446480        return out_list
    447481
    448 # #############################################################################   
     482
     483##############################################################################
    449484#
    450 # END OF CLASS
    451 #           
    452 # #############################################################################
    453 
     485# END OF CLASS
     486#
     487##############################################################################
    454488
    455489# Functions below here
    … …  
    478512def empty_string_to_none(string):
    479513    """Given a string, return None if string == "", and the string value
    480        (with type str) otherwise
     514    (with type str) otherwise
    481515    """
    482516
    … …  
    489523def string_or_float(value):
    490524    """ If value is unicode, return a value with type string.
    491         If value is a number, return it with type float
    492         If value is anything else, return value
     525    If value is a number, return it with type float
     526    If value is anything else, return value
    493527    """
    494     if type(value) == unicode:
     528
     529    if isinstance(value, unicode):
    495530        return str(value)
    496531    elif isinstance(value, (int, long, float, complex)):
    … …  
    498533    else:
    499534        return value
    500 

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup/prepare_data.py

    r9343 r9344  
    2121# Local modules
    2222import user_functions
    23 from parse_input_data import ProjectData
     23from setup.parse_input_data import ProjectData
     24
    2425
    2526class PrepareData(ProjectData):
    … …  
    2728    def __init__(self, filename):
    2829        """Parse the input data then process it for ANUGA
    29              
    30         """
    31         # Get the 'raw' data       
     30
     31        """
     32        # Get the 'raw' data
    3233        ProjectData.__init__(self, filename)
    3334
    … …  
    7071
    7172        # Hack to override resolution
    72         # region_point_areas=[ region_point_areas[i][0:2]+[150*150*0.5] for i in \
     73        # region_point_areas=\
     74        # [ region_point_areas[i][0:2]+[150*150*0.5] for i in \
    7375        #                      range(len(region_point_areas))]
    7476
    … …  
    7678
    7779        self.interior_regions = [[su.read_polygon(ir[0]), ir[1]] for ir in
    78                                   self.interior_regions_data]
    79 
    80         # Deal with intersections in the bounding polygon / breaklines / riverwalls
     80                                 self.interior_regions_data]
     81
     82        # Deal with intersections in the bounding polygon / breaklines /
     83        # riverwalls
    8184
    8285        (self.bounding_polygon, self.breaklines, self.riverwalls) = \
    8386            su.add_intersections_to_domain_features(
    84                 self.bounding_polygon, 
    85                 self.breaklines, 
     87                self.bounding_polygon,
     88                self.breaklines,
    8689                self.riverwalls,
    87                 point_movement_threshold=\
    88                     self.break_line_intersect_point_movement_threshold,
     90                point_movement_threshold=self.break_line_intersect_point_movement_threshold,
    8991                verbose=True)
    9092
    … …  
    113115        #
    114116
    115         if type(self.projection_information) is int:
     117        if isinstance(self.projection_information, int):
    116118
    117119            # projection_information describes a UTM zone
    … …  
    126128                    + str(self.projection_information) \
    127129                    + ' +datum=WGS84 +units=m +no_defs'
    128         elif type(self.projection_information) is str:
     130        elif isinstance(self.projection_information, str):
    129131            self.proj4string = self.projection_information
    130132        else:
    … …  
    146148
    147149        self.output_dir = self.output_basedir + 'RUN_' + str(runtime) +\
    148                           '_' + self.scenario
     150            '_' + self.scenario
    149151        self.partition_basedir = 'PARTITIONS/'
    150152        self.partition_dir = self.partition_basedir + 'Mesh_' +\
    151                              str(self.mesh_id_hash)
     153            str(self.mesh_id_hash)
    152154        self.meshname = self.output_dir + '/mesh.tsh'
    153 
    154155
    155156    def setup_directory_structure(self):
    … …  
    215216                for i in range(len(self.breaklines)):
    216217                    sav_lines(self.breaklines.values()[i],
    217                               filename=self.spatial_text_output_dir +\
    218                                          '/breakline_' + str(i) + '.txt')
     218                              filename=self.spatial_text_output_dir +
     219                              '/breakline_' + str(i) + '.txt')
    219220            if self.riverwalls is not {}:
    220221                for i in range(len(self.riverwalls)):
    221222                    sav_lines(self.riverwalls.values()[i],
    222                               filename=self.spatial_text_output_dir +\
    223                                  '/riverwall_' + str(i) + '.txt')
    224 
    225             sav_lines(self.bounding_polygon, 
     223                              filename=self.spatial_text_output_dir +
     224                              '/riverwall_' + str(i) + '.txt')
     225
     226            sav_lines(self.bounding_polygon,
    226227                      filename=self.spatial_text_output_dir
    227228                      + '/boundingPoly.txt')
    … …  
    235236            # directory
    236237
    237             shutil.copytree(self.spatial_text_output_dir, 
    238                             self.output_dir +\
    239                              '/' + self.spatial_text_output_dir)
     238            shutil.copytree(self.spatial_text_output_dir,
     239                            self.output_dir +
     240                            '/' + self.spatial_text_output_dir)
    240241        else:
    241242            pass

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/setup/setup_boundary_conditions.py

    r9343 r9344  
    5151                skip_header=1)
    5252            # Make start time = 0
    53             boundary_data[:,0]-=start_time
     53            boundary_data[:, 0] -= start_time
    5454            # Make interpolation function
    5555            stage_time_fun = scipy.interpolate.interp1d(
    56                 boundary_data[:,0],
    57                 boundary_data[:,1])
    58             # Make boundary condition 
     56                boundary_data[:, 0],
     57                boundary_data[:, 1])
     58            # Make boundary condition
    5959            boundary_function = \
    6060                asb.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(
    6161                    domain,
    6262                    stage_time_fun)
    63            
     63
    6464            boundary_tags_and_conditions[boundary_tag] = \
    6565                boundary_function
    … …  
    6868            msg = 'Boundary condition type ' + boundary_condition_type +\
    6969                  ' for tag ' + boundary_tag + ' is not implemented'
    70 
    71                        
    7270
    7371    # Check that all tags in boundary_info have been set
    … …  
    8078            raise Exception(msg)
    8179
    82     # Set the boundary 
     80    # Set the boundary
    8381
    8482    domain.set_boundary(boundary_tags_and_conditions)

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami3/scripts/user_functions.py

    r9343 r9344  
    1515"""
    1616
    17 
     17import os
    1818import anuga
    1919from anuga_parallel import myid, numprocs, finalize, barrier
     20import ogr
    2021
    2122
    … …  
    6566    """
    6667
     68    if not os.path.exists(shapefile_name):
     69        msg = 'Cannot find file ' + shapefile_name
     70        raise ValueError(msg)
     71
    6772    # Step 1: Read the data
    68 
    69     import ogr
    70 
    7173    data_source = ogr.Open(shapefile_name)
    7274
    … …  
    9799    # by dropping the last point
    98100    if bounding_polygon[0] == bounding_polygon[-1]:
    99         bounding_polygon = [bounding_polygon[k] for k in range(len(bounding_polygon)-1)]
     101        lbp = len(bounding_polygon)-1
     102        bounding_polygon = [bounding_polygon[k] for k in range(lbp)]
    100103
    101104    for i in range(boundary_segnum - 1):
    … …  
    103106        # starts at the end point of bounding_polygon
    104107
    105         found_match = False # Flag for error if no matching segments found
     108        found_match = False  # Flag for error if no matching segments found
    106109
    107110        for j in range(boundary_segnum):
    … …  
    126129                              'polygon are not identical'
    127130                        raise ValueError(msg)
    128                    
     131
    129132                    # Remove the last point
    130133                    bounding_polygon = [bounding_polygon[k] for k in
    … …  
    147150        # Check that a 'match' was found for every 'i' iteration
    148151        if not found_match:
    149             raise Exception, 'Did not find a match'
     152            raise Exception('Did not find a match')
    150153
    151154    return (bounding_polygon, boundary_tags)
Note: See TracChangeset for help on using the changeset viewer.