Changeset 9322


Ignore:
Timestamp:
Sep 4, 2014, 1:36:40 PM (11 years ago)
Author:
davies
Message:

Template scripts: Allowing earthquake sources / add_quantities + various
tif related changes

Location:
trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py

    r9321 r9322  
    4040yieldStep=20.0
    4141
    42 
     42## Spatial Information
     43#
     44# Give UTM Zone as an integer (negative for Southern Hemisphere)
     45utm_zone = 54
    4346
    4447#####################################################################
    … …  
    141144# pixle size (m) for output tifs
    142145output_tif_cellsize = 250.
    143 # EPSG CODE (projection information) for outputs tifs
    144 EPSG_CODE = 32647
     146# Polygon for output tif [either a filename, or if None the bounding_polygon is used]
     147output_tif_bounding_polygon = None
     148
    145149
    146150## Storage of max quantities
    … …  
    214218            [bounding_polygon, interior_regions, riverWalls, \
    215219             breakLines, RegionPtAreas, default_res, BoundaryInfo]) ).hexdigest()
     220
     221# Fix the output tif bounding polygon
     222if output_tif_bounding_polygon is None:
     223    output_tif_bounding_polygon = bounding_polygon
     224else:
     225    output_tif_bounding_polygon = su.read_polygon(output_tif_bounding_polygon)
    216226
    217227# Set up directories etc

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/raster_outputs.py

    r9321 r9322  
    1010import numpy
    1111import gdal
     12from pyproj import Proj
    1213from anuga.utilities import plot_utils as util
    1314from anuga.utilities import spatialInputUtil as su
    … …  
    1920#
    2021
    21 def make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir):
     22def make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, proj4string , tif_output_dir):
    2223    """
    2324        Call gdalwarp to get the resampled elevation data
    … …  
    4243    yres=CellSize
    4344
    44     allSRS='EPSG:'+str(EPSG_CODE)
     45    #allSRS='EPSG:'+str(EPSG_CODE)
     46    allSRS = proj4string
    4547    newDEM=tif_output_dir+'LIDAR_resampled.tif'
    4648    newMask=tif_output_dir+'Mask_resampled.tif'
    … …  
    6668
    6769    return
     70
     71##############################################################################
    6872 
    6973def gdal_calc_command(stage, depth, elevation, mask, depth_threshold):
    … …  
    8993def make_me_some_tifs(swwFile,
    9094                    bounding_polygon,
    91                     EPSG_CODE,
     95                    utm_zone,
    9296                    myTimeStep='collected_max',
    9397                    tif_output_subdir='/TIFS/',
    … …  
    104108    - swwFile -- Full path name of sww to read outputs from
    105109    - bounding_polygon -- ANUGA's bounding polygon, or another polygon to clip the rasters to
    106     - EPSG_CODE -- EPSG projection info as an integer like 3123 [which is over Manila, Philippines]
     110    - utm_zone -- UTM zone as an integer [make negative for Southern Hemisphere]
    107111    - myTimeStep -- option from util.Make_Geotif
    108112       use 'max' to plot the maxima
    … …  
    125129    """
    126130
     131    # Convert utm_zone to proj4string
     132    p=Proj(proj='utm', south=(utm_zone<0.), zone=abs(utm_zone), ellps='WGS84')
     133    proj4string = p.srs
     134
    127135    tif_output_dir=os.path.dirname(swwFile)+tif_output_subdir
    128136    # Make the geotifs
    … …  
    144152            tmpArr=maxQs[:,[0, 1, i+2]]
    145153            util.Make_Geotif(tmpArr, output_quantities=[quant+'_MAX'],
    146                              CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True,
     154                             CellSize=CellSize, proj4string = proj4string , verbose=True,
    147155                             bounding_polygon=bounding_polygon, output_dir=tif_output_dir)
    148156       
    … …  
    150158        util.Make_Geotif(swwFile, myTimeStep=[0],
    151159                         output_quantities=['elevation', 'friction'],
    152                          CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True,
     160                         CellSize=CellSize, proj4string = proj4string , verbose=True,
    153161                         bounding_polygon=bounding_polygon,
    154162                         output_dir=tif_output_dir)
    … …  
    156164        util.Make_Geotif(swwFile, myTimeStep=myTimeStep,
    157165                         output_quantities=['depth', 'stage', 'elevation', 'velocity', 'depthIntegratedVelocity', 'friction'],
    158                          CellSize=CellSize, EPSG_CODE=EPSG_CODE, verbose=True,
     166                         CellSize=CellSize, proj4string = proj4string , verbose=True,
    159167                         bounding_polygon=bounding_polygon,
    160168                         output_dir=tif_output_dir)
    … …  
    169177   
    170178    # Make the resampled elevation data
    171     make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, EPSG_CODE, tif_output_dir)
     179    make_resampled_elevation(elevation_raster, rasterExtent, CellSize, clipPolygon, clipPolygonLayer, proj4string, tif_output_dir)
    172180   
    173181    elevation=glob.glob(tif_output_dir+'LIDAR_resampled*')[0]

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

    r9321 r9322  
    6262#
    6363####################################################################################
    64 print 'Making rainfall conditions'
    65 setup_rainfall.setup_rainfall(domain, project)
    66 print 'Making inlets '
    67 setup_inlets.setup_inlets(domain, project)
     64# print 'Making rainfall conditions'
     65# setup_rainfall.setup_rainfall(domain, project)
     66# print 'Making inlets '
     67# setup_inlets.setup_inlets(domain, project)
    6868
    6969#####################################################################################
    … …  
    138138        raster_outputs.make_me_some_tifs(swwFile = './'+project.scenario+'.sww',
    139139                                    bounding_polygon=project.bounding_polygon,
    140                                     EPSG_CODE = project.EPSG_CODE,
     140                                    utm_zone = project.utm_zone,
    141141                                    CellSize= project.output_tif_cellsize)
    142142    except:

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

    r9320 r9322  
    4949        THE USER MUST MODIFY THIS TO SETUP THE BOUNDARIES AS DESIRED
    5050
    51         1) First make sure the right boundary conditions are created
     51        1) First make sure the required boundary conditions are created
    5252           (see examples below, modify or add new ones as required)
    5353

  • trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_initial_conditions.py

    r9320 r9322  
    6363    #START EDITING HERE
    6464    #@@@@@@@@@@@@@@@@@#
     65
     66
     67    ##############################################################################
     68    #
     69    # PRELIMINARY DEFINITIONS
     70    #
     71    ##############################################################################
     72
     73    # If the earthquake slip is given then convert to deformation here
     74    # So that we can add it to both elevation AND stage
     75    make_earthquake_deformation_from_source_model=False
     76
     77    if make_earthquake_deformation_from_source_model:
     78        source_model='../anuga/source/XXXXX.csv'
     79        # Make an earthquake source model
     80        eq_function = okada_tsunami.earthquake_source(source=None,
     81                                                  filename=source_model,
     82                                                  domain=domain,
     83                                                  lon_lat_degrees=True,
     84                                                  lon_before_lat=False,
     85                                                  utm_zone=project.utm_zone,
     86                                                  verbose=False)
     87    else:
     88        # Make a function that is zero everywhere
     89        # If we pass this to add_quantity, nothing should happen
     90        def eq_function(x,y):
     91            return x*0.
    6592
    6693
    … …  
    80107    elevation_default = numpy.genfromtxt('../anuga/topographies/andaman_3s_grid_Clip2_Project.txt', delimiter=",",skip_header=1)
    81108
    82 
     109    # Set_quantity data
    83110    elevation_Data=[
    84111                    ['../anuga/polygons/patong_10m.txt', elevation_patong_10m],
    … …  
    88115                   ]
    89116
     117    # Define a function for add_quantity in the same way
     118    # [ ['All', 0.] ] will do nothing
     119    elevation_additions = [
     120                            ['All', eq_function]
     121                          ]
     122
    90123   
    91124    ##############################################################################
    … …  
    96129   
    97130    default_friction=0.025
     131    # Set_quantity data
    98132    friction_Data =[
    99133                    ['All', default_friction]
    100134                   ]
     135   
     136    # Add_quantity data
     137    friction_additions = [
     138                          ['All', 0.]
     139                         ]
    101140   
    102141    ###############################################################################
    … …  
    106145    ###############################################################################
    107146
     147    # Set_quantity data
    108148    stage_Data = [
    109149                  ['All', project.tide]
    110150                 ]
    111151
     152   
     153    # Add_quantity data
     154    stage_additions = [
     155                       ['All', eq_function]
     156                      ]
    112157    ###############################################################################
    113158    #
    … …  
    117162    ###############################################################################
    118163   
     164    # Set_quantity data
    119165    xmomentum_Data = [
    120166                      ['All', 0.]
    121167                     ]
    122 
     168   
     169    # Add_quantity data
     170    xmomentum_additions = [
     171                            ['All', 0.]
     172                          ]
     173
     174    # Set_quantity data
    123175    ymomentum_Data = [
    124176                      ['All', 0.]
    125177                     ]
     178   
     179    # Add_quantity data
     180    ymomentum_additions = [
     181                            ['All', 0.]
     182                          ]
    126183
    127184    #@@@@@@@@@@@@@@@@@#
    … …  
    140197    elevFun=qs.composite_quantity_setting_function(elevation_Data, domain)
    141198    domain.set_quantity('elevation', elevFun, location='vertices')
     199    elevAddFun=qs.composite_quantity_setting_function(elevation_additions, domain)
     200    domain.add_quantity('elevation', elevAddFun)
    142201   
    143202    frictionFun=qs.composite_quantity_setting_function(friction_Data, domain)
    144203    domain.set_quantity('friction',  frictionFun,  location='centroids')
     204    frictionAddFun=qs.composite_quantity_setting_function(friction_additions, domain)
     205    domain.add_quantity('friction', frictionAddFun)
    145206   
    146207    stageFun=qs.composite_quantity_setting_function(stage_Data, domain)
    147208    domain.set_quantity('stage',  stageFun,  location='centroids')
     209    stageAddFun=qs.composite_quantity_setting_function(stage_additions, domain)
     210    domain.add_quantity('stage', stageAddFun)
    148211
    149212    xmomentumFun=qs.composite_quantity_setting_function(xmomentum_Data, domain)
    150213    domain.set_quantity('xmomentum',  xmomentumFun,  location='centroids')
     214    xmomentumAddFun=qs.composite_quantity_setting_function(xmomentum_additions, domain)
     215    domain.add_quantity('xmomentum', xmomentumAddFun)
    151216   
    152217    ymomentumFun=qs.composite_quantity_setting_function(ymomentum_Data, domain)
    153218    domain.set_quantity('ymomentum',  ymomentumFun,  location='centroids')
     219    ymomentumAddFun=qs.composite_quantity_setting_function(ymomentum_additions, domain)
     220    domain.add_quantity('ymomentum', ymomentumAddFun)
    154221
    155222    ################################################################################
Note: See TracChangeset for help on using the changeset viewer.