source: anuga_work/production/australia_ph2/ceduna/Arc_asc2raster_GDA94z50.py @ 6842

Last change on this file since 6842 was 6842, checked in by myall, 15 years ago

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.

File size: 4.0 KB
Line 
1# ---------------------------------------------------------------------------
2# This python script is an ArcGIS script that can only be run on a computer
3# with and ArcGIS licence and version 2.4.1 python.
4# This script is designed to read in .asc files and deliever rasters with
5# projection (GDA94z50) held in a file geodatabase (called raster)
6# written by Kristy Van Putten and Ross Wilson
7# ---------------------------------------------------------------------------
8
9# Import system modules
10import sys, string, os, arcgisscripting, glob, os.path
11
12# Create the Geoprocessor object
13gp = arcgisscripting.create()
14
15# Check out any necessary licenses
16gp.CheckOutExtension("spatial")
17
18# Load required toolboxes...
19gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Spatial Analyst Tools.tbx")
20gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Conversion Tools.tbx")
21gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")
22gp.overwriteoutput = 1
23
24scenario_dir = "\\\\nas2\\gemd\\georisk_models\\inundation\\data\\australia_ph2\\ceduna\\"
25output_dir="anuga\\outputs\\"
26
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'
30
31events = [[time_dir1,0.314774],[time_dir2,0.27307],[time_dir3,0.315495]]
32
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
39    # Local variables...
40    folder = scenario_dir + output_dir +  time_dir +'\\'
41    raster_gbd = folder + 'raster.gdb'
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"
44
45    print 'Process: Create File GDB'
46    gp.CreateFileGDB_management(folder, "raster")
47
48    gp.Workspace = raster_gbd
49
50    print gp.Workspace
51   
52    #replication dictionary
53    replicate = (('ceduna', ''),('_', ''),('max','_M'),
54                 ('CBD', 'CDB'),('All',''),
55                 ('depth','_depth'),('speed', '_speed'),
56                 ('elevation', '_ele_'), ('stage','_stage'))
57
58    generate_filename = []
59    input_ascii = glob.glob(folder + '*max.asc')
60
61    for infile in input_ascii:
62        output_DEM = os.path.basename(infile)[:-4]
63        for (key, rep) in replicate:
64            output_DEM = output_DEM.replace(key,rep)
65        output_DEM = output_DEM[:10]
66        if output_DEM in generate_filename:
67            print 'Output_DEM filename (%s) already in use' % output_DEM
68            sys.exit(10)
69        generate_filename.append(output_DEM)
70
71        print 'Output DEM ',output_DEM
72       
73        print 'Process: ASCII to Raster'
74        gp.ASCIIToRaster_conversion(infile, output_DEM, "FLOAT")
75
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]"
80                                                   ",PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]]")
81##        output_extract = output_DEM + 'E'
82##        print 'Output Extract ',output_extract
83##        print 'Process: Extract by Mask'
84##        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)
90
Note: See TracBrowser for help on using the repository browser.