Changeset 6842 for anuga_work


Ignore:
Timestamp:
Apr 20, 2009, 4:45:56 PM (16 years ago)
Author:
myall
Message:

using ARC script to turn ascii stage files into rasters, and dividing by max wave height; done for GC clockwise to perth. not checked thoroughly.

Location:
anuga_work/production/australia_ph2
Files:
4 added
7 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/australia_ph2/albany/Arc_asc2raster_GDA94z50.py

    r6484 r6842  
    2020gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
    2121gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    22 gp.overwriteoutput = 0
     22gp.overwriteoutput = 1
    2323
    24 scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\ceduna\\"
     24scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\albany\\"
    2525output_dir="anuga\\outputs\\"
    2626
    27 time_dir1 = '20090225_153629_run_final_0_64469_mhingee'
     27time_dir1 = '20090408_152157_run_final_0_27319_1330_Tb__mhingee'
     28time_dir2 = '20090408_201701_run_final_0_64344_1330_Tb__mhingee'
     29time_dir3 = '20090409_010309_run_final_0_46697_1330_Tb__mhingee'
    2830
    29 time_dirs = [time_dir1] #, time_dir2] #4, time_dir5] #, time_dir3, time_dir4, time_dir5, time_dir6]
    30  
    31 for time_dir in time_dirs:
     31events = [[time_dir1,0.395764],[time_dir2,0.390969],[time_dir3,0.392536]]
    3232
     33for event in events:
     34##for time_dir in time_dirs:
     35    time_dir = event[0]
     36    max_wave = event[1]
     37    print time_dir
     38    print max_wave
    3339    # Local variables...
    34     folder = scenario_dir + output_dir + time_dir +'\\'
     40    folder = scenario_dir + output_dir +  time_dir +'\\'
    3541    raster_gbd = folder + 'raster.gdb'
    36     land = scenario_dir + "map_work\\Busselton.gdb\\Internal_polygons\\initial_conditions_extend_Di"
    37     ocean = scenario_dir + "map_work\\Busselton.gdb\\input_boundaries\\Ocean"
     42##    land = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_condition"
     43##    ocean = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_conditions_ocean"
    3844
    39    
    4045    print 'Process: Create File GDB'
    4146    gp.CreateFileGDB_management(folder, "raster")
    … …  
    4348    gp.Workspace = raster_gbd
    4449
    45     print time_dir
     50    print gp.Workspace
    4651   
    4752    #replication dictionary
    48     replicate = (('ceduna', ''),('_', ''),
     53    replicate = (('albany', ''),('_', ''),('max','_M'),
     54                 ('CBD', 'CDB'),('All',''),
    4955                 ('depth','_depth'),('speed', '_speed'),
    5056                 ('elevation', '_ele_'), ('stage','_stage'))
    5157
    5258    generate_filename = []
    53     input_ascii = glob.glob(folder + '*.asc')
     59    input_ascii = glob.glob(folder + '*stage.asc')
    5460
    5561    for infile in input_ascii:
    … …  
    5763        for (key, rep) in replicate:
    5864            output_DEM = output_DEM.replace(key,rep)
    59         output_DEM = output_DEM[:12]
     65        output_DEM = output_DEM[:10]
    6066        if output_DEM in generate_filename:
    6167            print 'Output_DEM filename (%s) already in use' % output_DEM
    … …  
    6874        gp.ASCIIToRaster_conversion(infile, output_DEM, "FLOAT")
    6975
    70         print 'Process: Define Projection'
    71         gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_53',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    72                                                    ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',135.0],PARAMETER['Scale_Factor',0.9996]"
     76        print 'Process: Define Projection'
     77        # GDA_1994_MGA_Zone_54
     78        gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_50',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
     79                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',117.0],PARAMETER['Scale_Factor',0.9996]"
    7380                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    74 ##        output_extract = output_DEM + '_E'
     81##        output_extract = output_DEM + 'E'
    7582##        print 'Output Extract ',output_extract
    7683##        print 'Process: Extract by Mask'
    7784##        gp.ExtractByMask_sa(output_DEM, land, output_extract)
     85       
     86        div_file = output_DEM.replace('stage','div_stage')
     87        if not div_file == output_DEM:
     88            print 'divide', output_DEM,' by', max_wave,' and call',div_file
     89            gp.Divide_sa(output_DEM,max_wave,div_file)
    7890
    79 

  • anuga_work/production/australia_ph2/ceduna/Arc_asc2raster_GDA94z50.py

    r6474 r6842  
    2020gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
    2121gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    22 gp.overwriteoutput = 0
     22gp.overwriteoutput = 1
    2323
    2424scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\ceduna\\"
    2525output_dir="anuga\\outputs\\"
    2626
    27 time_dir1 = '20090225_153629_run_final_0_64469_mhingee'
     27time_dir1 = '20090408_182701_run_final_0_27347_1884_Tb__kvanputt'
     28time_dir2 = '20090409_095708_run_final_0_64448_1884_Tb__kvanputt'
     29time_dir3 = '20090410_000802_run_final_0_58331_1884_Tb__kvanputt'
    2830
    29 time_dirs = [time_dir1] #, time_dir2] #4, time_dir5] #, time_dir3, time_dir4, time_dir5, time_dir6]
    30  
    31 for time_dir in time_dirs:
     31events = [[time_dir1,0.314774],[time_dir2,0.27307],[time_dir3,0.315495]]
    3232
     33for event in events:
     34##for time_dir in time_dirs:
     35    time_dir = event[0]
     36    max_wave = event[1]
     37    print time_dir
     38    print max_wave
    3339    # Local variables...
    34     folder = scenario_dir + output_dir + time_dir +'\\'
     40    folder = scenario_dir + output_dir +  time_dir +'\\'
    3541    raster_gbd = folder + 'raster.gdb'
    36     land = scenario_dir + "map_work\\Busselton.gdb\\Internal_polygons\\initial_conditions_extend_Di"
    37     ocean = scenario_dir + "map_work\\Busselton.gdb\\input_boundaries\\Ocean"
     42##    land = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_condition"
     43##    ocean = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_conditions_ocean"
    3844
    39    
    4045    print 'Process: Create File GDB'
    4146    gp.CreateFileGDB_management(folder, "raster")
    … …  
    4348    gp.Workspace = raster_gbd
    4449
    45     print time_dir
     50    print gp.Workspace
    4651   
    4752    #replication dictionary
    48     replicate = (('ceduna', ''),('_', ''),
     53    replicate = (('ceduna', ''),('_', ''),('max','_M'),
     54                 ('CBD', 'CDB'),('All',''),
    4955                 ('depth','_depth'),('speed', '_speed'),
    5056                 ('elevation', '_ele_'), ('stage','_stage'))
    5157
    5258    generate_filename = []
    53     input_ascii = glob.glob(folder + '*.asc')
     59    input_ascii = glob.glob(folder + '*max.asc')
    5460
    5561    for infile in input_ascii:
    … …  
    5763        for (key, rep) in replicate:
    5864            output_DEM = output_DEM.replace(key,rep)
    59         output_DEM = output_DEM[:12]
     65        output_DEM = output_DEM[:10]
    6066        if output_DEM in generate_filename:
    6167            print 'Output_DEM filename (%s) already in use' % output_DEM
    … …  
    6874        gp.ASCIIToRaster_conversion(infile, output_DEM, "FLOAT")
    6975
    70         print 'Process: Define Projection'
    71         gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_53',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    72                                                    ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',135.0],PARAMETER['Scale_Factor',0.9996]"
     76        print 'Process: Define Projection'
     77        # GDA_1994_MGA_Zone_54
     78        gp.DefineProjection_management(output_DEM, "PROJCS['CM_133',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
     79                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',133.0],PARAMETER['Scale_Factor',0.9996]"
    7380                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    74 ##        output_extract = output_DEM + '_E'
     81##        output_extract = output_DEM + 'E'
    7582##        print 'Output Extract ',output_extract
    7683##        print 'Process: Extract by Mask'
    7784##        gp.ExtractByMask_sa(output_DEM, land, output_extract)
     85        # do this bit only if there are only stage asc files
     86        div_file = output_DEM.replace('stage','div_stage')
     87        if not div_file == output_DEM:
     88            print 'divide', output_DEM,' by', max_wave,' and call',div_file
     89            gp.Divide_sa(output_DEM,max_wave,div_file)
    7890
    79 

  • anuga_work/production/australia_ph2/esperance/Arc_asc2raster_GDA94z50.py

    r6474 r6842  
    2020gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
    2121gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    22 gp.overwriteoutput = 0
     22gp.overwriteoutput = 1
    2323
    2424scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\esperance\\"
    2525output_dir="anuga\\outputs\\"
    2626
    27 time_dir1 = '20090302_104336_run_final_0_27255_mhingee'
     27time_dir1 = '20090408_152320_run_final_0_27329_1699_Tb__mhingee'
     28time_dir2 = '20090408_205317_run_final_0_58367_1699_Tb__mhingee'
     29time_dir3 = '20090409_024308_run_final_0_64449_1699_Tb__mhingee'
    2830
    29 time_dirs = [time_dir1] #, time_dir2] #4, time_dir5] #, time_dir3, time_dir4, time_dir5, time_dir6]
    30  
    31 for time_dir in time_dirs:
     31events = [[time_dir1,0.352214],[time_dir2,0.356316],[time_dir3,0.357256]]
    3232
     33for event in events:
     34##for time_dir in time_dirs:
     35    time_dir = event[0]
     36    max_wave = event[1]
     37    print time_dir
     38    print max_wave
    3339    # Local variables...
    34     folder = scenario_dir + output_dir + time_dir +'\\'
     40    folder = scenario_dir + output_dir +  time_dir +'\\'
    3541    raster_gbd = folder + 'raster.gdb'
    36     land = scenario_dir + "map_work\\Busselton.gdb\\Internal_polygons\\initial_conditions_extend_Di"
    37     ocean = scenario_dir + "map_work\\Busselton.gdb\\input_boundaries\\Ocean"
     42##    land = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_condition"
     43##    ocean = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_conditions_ocean"
    3844
    39    
    4045    print 'Process: Create File GDB'
    4146    gp.CreateFileGDB_management(folder, "raster")
    … …  
    4348    gp.Workspace = raster_gbd
    4449
    45     print time_dir
     50    print gp.Workspace
    4651   
    4752    #replication dictionary
    48     replicate = (('esperance', ''),('_', ''),
     53    replicate = (('esperance', ''),('_', ''),('max','_M'),
     54                 ('CBD', 'CDB'),('All',''),
    4955                 ('depth','_depth'),('speed', '_speed'),
    5056                 ('elevation', '_ele_'), ('stage','_stage'))
    5157
    5258    generate_filename = []
    53     input_ascii = glob.glob(folder + '*.asc')
     59    input_ascii = glob.glob(folder + '*stage.asc')
    5460
    5561    for infile in input_ascii:
    … …  
    5763        for (key, rep) in replicate:
    5864            output_DEM = output_DEM.replace(key,rep)
    59         output_DEM = output_DEM[:12]
     65        output_DEM = output_DEM[:10]
    6066        if output_DEM in generate_filename:
    6167            print 'Output_DEM filename (%s) already in use' % output_DEM
    … …  
    6874        gp.ASCIIToRaster_conversion(infile, output_DEM, "FLOAT")
    6975
    70         print 'Process: Define Projection'
     76        print 'Process: Define Projection'
     77        # GDA_1994_MGA_Zone_54
    7178        gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_51',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    7279                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',123.0],PARAMETER['Scale_Factor',0.9996]"
    7380                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    74 ##        output_extract = output_DEM + '_E'
     81##        output_extract = output_DEM + 'E'
    7582##        print 'Output Extract ',output_extract
    7683##        print 'Process: Extract by Mask'
    7784##        gp.ExtractByMask_sa(output_DEM, land, output_extract)
     85        # do this bit only if there are only stage asc files
     86        div_file = output_DEM.replace('stage','div_stage')
     87        if not div_file == output_DEM:
     88            print 'divide', output_DEM,' by', max_wave,' and call',div_file
     89            gp.Divide_sa(output_DEM,max_wave,div_file)
    7890
    79 

  • anuga_work/production/australia_ph2/eucla_motel/Arc_asc2raster_GDA94z50.py

    r6474 r6842  
    2020gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
    2121gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    22 gp.overwriteoutput = 0
     22gp.overwriteoutput = 1
    2323
    2424scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\eucla_motel\\"
    2525output_dir="anuga\\outputs\\"
    2626
    27 time_dir1 = '20090225_153734_run_final_0_65371_mhingee'
     27time_dir1 = '20090408_152411_run_final_0_58332_1837_Tb__mhingee'
     28time_dir2 = '20090409_073114_run_final_0_64449_1837_Tb__mhingee'
     29time_dir3 = '20090409_222301_run_final_0_68755_1837_Tb__mhingee'
    2830
    29 time_dirs = [time_dir1] #, time_dir2] #4, time_dir5] #, time_dir3, time_dir4, time_dir5, time_dir6]
    30  
    31 for time_dir in time_dirs:
     31events = [[time_dir1,0.265183],[time_dir2,0.27574],[time_dir3,0.267985]]
    3232
     33for event in events:
     34##for time_dir in time_dirs:
     35    time_dir = event[0]
     36    max_wave = event[1]
     37    print time_dir
     38    print max_wave
    3339    # Local variables...
    34     folder = scenario_dir + output_dir + time_dir +'\\'
     40    folder = scenario_dir + output_dir +  time_dir +'\\'
    3541    raster_gbd = folder + 'raster.gdb'
    36     land = scenario_dir + "map_work\\Busselton.gdb\\Internal_polygons\\initial_conditions_extend_Di"
    37     ocean = scenario_dir + "map_work\\Busselton.gdb\\input_boundaries\\Ocean"
     42##    land = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_condition"
     43##    ocean = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_conditions_ocean"
    3844
    39    
    4045    print 'Process: Create File GDB'
    4146    gp.CreateFileGDB_management(folder, "raster")
    … …  
    4348    gp.Workspace = raster_gbd
    4449
    45     print time_dir
     50    print gp.Workspace
    4651   
    4752    #replication dictionary
    48     replicate = (('eucla_motel', ''),('_', ''),
     53    replicate = (('eucla_motel', ''),('_', ''),('max','_M'),
     54                 ('CBD', 'CDB'),('All',''),
    4955                 ('depth','_depth'),('speed', '_speed'),
    5056                 ('elevation', '_ele_'), ('stage','_stage'))
    5157
    5258    generate_filename = []
    53     input_ascii = glob.glob(folder + '*.asc')
     59    input_ascii = glob.glob(folder + '*max.asc')
    5460
    5561    for infile in input_ascii:
    … …  
    5763        for (key, rep) in replicate:
    5864            output_DEM = output_DEM.replace(key,rep)
    59         output_DEM = output_DEM[:12]
     65        output_DEM = output_DEM[:10]
    6066        if output_DEM in generate_filename:
    6167            print 'Output_DEM filename (%s) already in use' % output_DEM
    … …  
    6874        gp.ASCIIToRaster_conversion(infile, output_DEM, "FLOAT")
    6975
    70         print 'Process: Define Projection'
     76        print 'Process: Define Projection'
     77        # GDA_1994_MGA_Zone_54
    7178        gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_52',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    7279                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',129.0],PARAMETER['Scale_Factor',0.9996]"
    7380                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    74 ##        output_extract = output_DEM + '_E'
     81##        output_extract = output_DEM + 'E'
    7582##        print 'Output Extract ',output_extract
    7683##        print 'Process: Extract by Mask'
    7784##        gp.ExtractByMask_sa(output_DEM, land, output_extract)
     85        # do this bit only if there are only stage asc files
     86        div_file = output_DEM.replace('stage','div_stage')
     87        if not div_file == output_DEM:
     88            print 'divide', output_DEM,' by', max_wave,' and call',div_file
     89            gp.Divide_sa(output_DEM,max_wave,div_file)
    7890
    79 

  • anuga_work/production/australia_ph2/portland/Arc_asc2raster_GDA94z50.py

    r6716 r6842  
    2020gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
    2121gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    22 gp.overwriteoutput = 0
     22gp.overwriteoutput = 1
    2323
    24 scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\esperance\\"
     24scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\portland\\"
    2525output_dir="anuga\\outputs\\"
    2626
    27 time_dir1 = '20090302_104336_run_final_0_27255_mhingee'
     27time_dir1 = '20090408_152424_run_final_0_27289_1969_Tb__mhingee'
     28time_dir2 = '20090408_194222_run_final_0_58358_1969_Tb__mhingee'
     29time_dir3 = '20090409_001857_run_final_0_64322_1969_Tb__mhingee'
    2830
    29 time_dirs = [time_dir1] #, time_dir2] #4, time_dir5] #, time_dir3, time_dir4, time_dir5, time_dir6]
    30  
    31 for time_dir in time_dirs:
     31events = [[time_dir1,0.430643],[time_dir2,0.462943],[time_dir3,0.418183]]
    3232
     33for event in events:
     34##for time_dir in time_dirs:
     35    time_dir = event[0]
     36    max_wave = event[1]
     37    print time_dir
     38    print max_wave
    3339    # Local variables...
    34     folder = scenario_dir + output_dir + time_dir +'\\'
     40    folder = scenario_dir + output_dir +  time_dir +'\\'
    3541    raster_gbd = folder + 'raster.gdb'
    36     land = scenario_dir + "map_work\\Busselton.gdb\\Internal_polygons\\initial_conditions_extend_Di"
    37     ocean = scenario_dir + "map_work\\Busselton.gdb\\input_boundaries\\Ocean"
     42##    land = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_condition"
     43##    ocean = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_conditions_ocean"
    3844
    39    
    4045    print 'Process: Create File GDB'
    4146    gp.CreateFileGDB_management(folder, "raster")
    … …  
    4348    gp.Workspace = raster_gbd
    4449
    45     print time_dir
     50    print gp.Workspace
    4651   
    4752    #replication dictionary
    48     replicate = (('esperance', ''),('_', ''),
     53    replicate = (('portland', ''),('_', ''),('max','_M'),
     54                 ('CBD', 'CDB'),('All',''),
    4955                 ('depth','_depth'),('speed', '_speed'),
    5056                 ('elevation', '_ele_'), ('stage','_stage'))
    … …  
    5763        for (key, rep) in replicate:
    5864            output_DEM = output_DEM.replace(key,rep)
    59         output_DEM = output_DEM[:12]
     65        output_DEM = output_DEM[:10]
    6066        if output_DEM in generate_filename:
    6167            print 'Output_DEM filename (%s) already in use' % output_DEM
    … …  
    6975
    7076        print 'Process: Define Projection'
    71         gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_51',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    72                                                    ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',123.0],PARAMETER['Scale_Factor',0.9996]"
     77        gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_54',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
     78                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',141.0],PARAMETER['Scale_Factor',0.9996]"
    7379                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    74 ##        output_extract = output_DEM + '_E'
     80##        output_extract = output_DEM + 'E'
    7581##        print 'Output Extract ',output_extract
    7682##        print 'Process: Extract by Mask'
    7783##        gp.ExtractByMask_sa(output_DEM, land, output_extract)
     84        # do this bit only if there are only stage asc files
     85        div_file = output_DEM.replace('stage','div_stage')
     86        print 'divide', output_DEM,' by', max_wave,' and call',div_file
     87        gp.Divide_sa(output_DEM,max_wave,div_file)
    7888
    79 

  • anuga_work/production/australia_ph2/strahan/Arc_asc2raster_GDA94z50.py

    r6619 r6842  
    2020gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
    2121gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
    22 gp.overwriteoutput = 0
     22gp.overwriteoutput = 1
    2323
    24 scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\esperance\\"
     24scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\strahan\\"
    2525output_dir="anuga\\outputs\\"
    2626
    27 time_dir1 = '20090302_104336_run_final_0_27255_mhingee'
     27time_dir1 = '20090408_152513_run_final_0_58337_2044_Tb__mhingee'
     28time_dir2 = '20090408_163848_run_final_0_64214_2044_Tb__mhingee'
     29time_dir3 = '20090408_175114_run_final_0_68779_2044_Tb__mhingee'
    2830
    29 time_dirs = [time_dir1] #, time_dir2] #4, time_dir5] #, time_dir3, time_dir4, time_dir5, time_dir6]
    30  
    31 for time_dir in time_dirs:
     31time_dirs = [time_dir1, time_dir2, time_dir3]
     32events = [[time_dir1,0.593696],[time_dir2,0.566183],[time_dir3,0.561357]]
    3233
     34for event in events:
     35##for time_dir in time_dirs:
     36    time_dir = event[0]
     37    max_wave = event[1]
     38    print time_dir
     39    print max_wave
    3340    # Local variables...
    34     folder = scenario_dir + output_dir + time_dir +'\\'
     41    folder = scenario_dir + output_dir +  time_dir +'\\'
    3542    raster_gbd = folder + 'raster.gdb'
    36     land = scenario_dir + "map_work\\Busselton.gdb\\Internal_polygons\\initial_conditions_extend_Di"
    37     ocean = scenario_dir + "map_work\\Busselton.gdb\\input_boundaries\\Ocean"
     43##    land = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_condition"
     44##    ocean = scenario_dir + "map_work\\\port_hedland.gdb\\outlines\\initial_conditions_ocean"
    3845
    39    
    4046    print 'Process: Create File GDB'
    4147    gp.CreateFileGDB_management(folder, "raster")
    … …  
    4349    gp.Workspace = raster_gbd
    4450
    45     print time_dir
     51    print gp.Workspace
    4652   
    4753    #replication dictionary
    48     replicate = (('esperance', ''),('_', ''),
     54    replicate = (('strahan', ''),('_', ''),('max','_M'),
     55                 ('CBD', 'CDB'),('All',''),
    4956                 ('depth','_depth'),('speed', '_speed'),
    5057                 ('elevation', '_ele_'), ('stage','_stage'))
    … …  
    5764        for (key, rep) in replicate:
    5865            output_DEM = output_DEM.replace(key,rep)
    59         output_DEM = output_DEM[:12]
     66        output_DEM = output_DEM[:10]
    6067        if output_DEM in generate_filename:
    6168            print 'Output_DEM filename (%s) already in use' % output_DEM
    … …  
    6976
    7077        print 'Process: Define Projection'
    71         gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_51',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    72                                                    ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',123.0],PARAMETER['Scale_Factor',0.9996]"
     78        gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_55',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
     79                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',147.0],PARAMETER['Scale_Factor',0.9996]"
    7380                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    74 ##        output_extract = output_DEM + '_E'
     81##        output_extract = output_DEM + 'E'
    7582##        print 'Output Extract ',output_extract
    7683##        print 'Process: Extract by Mask'
    7784##        gp.ExtractByMask_sa(output_DEM, land, output_extract)
     85        # do this bit only if there are only stage asc files
     86        div_file = output_DEM.replace('stage','div_stage')
     87        print 'divide', output_DEM,' by', max_wave,' and call',div_file
     88        gp.Divide_sa(output_DEM,max_wave,div_file)
    7889
    79 

  • anuga_work/production/australia_ph2/sydney/Arc_asc2raster_GDA94z56.py

    r6741 r6842  
    2525output_dir = "anuga\\outputs\\"
    2626
    27 time_dir1 = '20090331_150428_run_final_0_58152_2956_Tb_lfountai'
    28 ##time_dir2 = '20081210_100528_run_final_0_68693_250m_none_kvanputt'
    29 ##time_dirs = [time_dir1, time_dir2]
     27time_dir1 = '20090408_171539_run_final_0_58152_2938_Tb__kvanputt'
     28time_dir2 = '20090409_011331_run_final_0_58222_2938_Tb__kvanputt'
     29time_dir3 = '20090409_085101_run_final_0_58326_2938_Tb__kvanputt'
     30time_dirs = [time_dir1, time_dir2, time_dir3]
    3031
    3132##time_dir1 = '20090311_145740_run_final_0_58152_2956_Tb_jgriffin'
    … …  
    3738
    3839
    39 time_dirs = [time_dir1]#, time_dir4 , time_dir3, time_dir4, time_dir5]#, time_dir6]
     40##time_dirs = [time_dir1]#, time_dir4 , time_dir3, time_dir4, time_dir5]#, time_dir6]
    4041
    4142for time_dir in time_dirs:
    … …  
    5960                 ('_', ''),('Geordie', 'Geo'),('Sorrento', 'Sor'), ('max','M_'),
    6061                 ('Fremantle', 'Fre'),('Rockingham', 'Roc'),('depth','_dep_'),
    61                  ('speed', '_spe_'), ('elevation', '_ele_'), ('stage','_stage'))
     62                 ('speed', '_spe_'), ('elevation', '_ele_'), ('stage','_st'))
    6263
    6364    generate_filename = []
    … …  
    8081
    8182        print 'Process: Define Projection'
    82         gp.DefineProjection_management(output_DEM, "PROJCS['GDA_1994_MGA_Zone_56',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    83                                                    ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',153.0],PARAMETER['Scale_Factor',0.9996]"
     83        gp.DefineProjection_management(output_DEM, "PROJCS['CM_151',GEOGCS['GCS_GDA_1994',DATUM['D_GDA_1994',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
     84                                                   ",PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',10000000.0],PARAMETER['Central_Meridian',151.0],PARAMETER['Scale_Factor',0.9996]"
    8485                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
    8586
Note: See TracChangeset for help on using the changeset viewer.