1 | """ |
---|
2 | Generates ascii grids of nominated areas - |
---|
3 | Input: sww file from run_perth.py |
---|
4 | boundaries for grids from project.py |
---|
5 | Outputs: ascii grids of specified variables |
---|
6 | Stored in the 'outputs_dir' folder for respective .sww file |
---|
7 | |
---|
8 | Note: |
---|
9 | If producing a grid for the enitre extent cellsize should be greater than 30m |
---|
10 | If producing grids for inundation area resolution should be greater than mesh (ie ~22m) |
---|
11 | """ |
---|
12 | |
---|
13 | import glob |
---|
14 | import os |
---|
15 | import sys |
---|
16 | from os.path import join |
---|
17 | from os import sep |
---|
18 | import anuga |
---|
19 | from anuga_parallel import myid, numprocs, distribute, finalize, barrier |
---|
20 | |
---|
21 | import time |
---|
22 | t0 = time.time() |
---|
23 | |
---|
24 | rootDir = r"/short/w85/Shane/Output" |
---|
25 | scenario_name = "Busselton" |
---|
26 | event = "B2Full" |
---|
27 | |
---|
28 | time_dirs = [#"final_"+event+"_sfm547_np1",\ |
---|
29 | # "final_"+event+"_sfm547_np1_From28800",\ |
---|
30 | # "final_"+event+"_sfm547_np1_From57600",\ |
---|
31 | # "final_"+event+"_sfm547_np1_From86400",\ |
---|
32 | # "final_"+event+"_sfm547_np1_From108000",\ |
---|
33 | "final_"+event+"_sfm547_np1_From115200"] |
---|
34 | |
---|
35 | #Modify the cellsize value to set the size of the raster you require |
---|
36 | #Take into account mesh size when aplying this paramater |
---|
37 | cellsize = 25 #250 |
---|
38 | |
---|
39 | #Now set the timestep at which you want the raster generated. |
---|
40 | #Either set the actual timestep required or use 'None' to indicate that |
---|
41 | #you want the maximum values in the raster over all timesteps |
---|
42 | timestep = None |
---|
43 | |
---|
44 | # Set the special areas of interest. If none, do: area=['All'] |
---|
45 | # Areas identified below are specified by cooridinates in project.py |
---|
46 | |
---|
47 | #area = ['South', 'NW', 'Hobart'] ### JANE - CHECK THIS ##### |
---|
48 | area = ['All'] |
---|
49 | |
---|
50 | #one or more key strings from var_equations below |
---|
51 | #var = ['stage','depth','momentum','elevation'] |
---|
52 | vars = ['depth','momentum'] |
---|
53 | |
---|
54 | ###### |
---|
55 | # Define allowed variable names and associated equations to generate values. |
---|
56 | # This would not normally change. |
---|
57 | ###### |
---|
58 | var_equations = {'stage': 'stage', |
---|
59 | 'momentum': '(xmomentum**2 + ymomentum**2)**0.5', |
---|
60 | 'depth': 'stage-elevation', |
---|
61 | 'speed': '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)', |
---|
62 | 'elevation': 'elevation' } |
---|
63 | |
---|
64 | ###### |
---|
65 | # Start script, running through variables, area, folder, sww file (determine by times) |
---|
66 | # Then creates a maximum asci if there is more than one sww file |
---|
67 | ###### |
---|
68 | |
---|
69 | proc = 0 |
---|
70 | Needed = True |
---|
71 | |
---|
72 | for dirID in range(len(time_dirs)): |
---|
73 | |
---|
74 | if myid == proc: |
---|
75 | time_dir = time_dirs[dirID] |
---|
76 | var = vars[0] |
---|
77 | |
---|
78 | elif myid == proc + 1: |
---|
79 | time_dir = time_dirs[dirID] |
---|
80 | var = vars[1] |
---|
81 | |
---|
82 | proc+=2 |
---|
83 | |
---|
84 | if myid >= len(time_dirs)*2: |
---|
85 | print "I am not needed: Proc ", myid |
---|
86 | Needed = False |
---|
87 | |
---|
88 | if Needed: |
---|
89 | |
---|
90 | if var not in var_equations: |
---|
91 | print 'Unrecognized variable name: %s' % which_var |
---|
92 | sys.exit() |
---|
93 | |
---|
94 | for which_area in area: |
---|
95 | if which_area == 'All': |
---|
96 | easting_min = None |
---|
97 | easting_max = None |
---|
98 | northing_min = None |
---|
99 | northing_max = None |
---|
100 | else: |
---|
101 | try: |
---|
102 | easting_min = eval('project.xmin%s' % which_area) |
---|
103 | easting_max = eval('project.xmax%s' % which_area) |
---|
104 | northing_min = eval('project.ymin%s' % which_area) |
---|
105 | northing_max = eval('project.ymax%s' % which_area) |
---|
106 | except AttributeError: |
---|
107 | print 'Unrecognized area name: %s' % which_area |
---|
108 | break |
---|
109 | |
---|
110 | names = [] |
---|
111 | Name = glob.glob(join(rootDir,event,time_dir,scenario_name) +'*.sww') |
---|
112 | names.append(Name) |
---|
113 | |
---|
114 | for names in names: |
---|
115 | |
---|
116 | asc_name = [] |
---|
117 | print 'start sww2dem: time_dir = %s' % time_dir |
---|
118 | for name in names: |
---|
119 | |
---|
120 | test = name |
---|
121 | name = test[:-4] |
---|
122 | outname = name + '_' + which_area + '_' + var + str(cellsize) + "m" |
---|
123 | quantityname = var_equations[var] |
---|
124 | |
---|
125 | print 'start sww2dem: name = %s' % os.path.basename(name) |
---|
126 | print outname |
---|
127 | anuga.sww2dem(name+'.sww', outname+'.asc', |
---|
128 | quantity = quantityname, |
---|
129 | cellsize = cellsize, |
---|
130 | easting_min = easting_min, |
---|
131 | easting_max = easting_max, |
---|
132 | northing_min = northing_min, |
---|
133 | northing_max = northing_max, |
---|
134 | reduction = timestep, |
---|
135 | verbose = True) |
---|
136 | |
---|
137 | asc_name.append(outname + '.asc') |
---|
138 | |
---|
139 | barrier() |
---|
140 | |
---|
141 | if myid == 0: |
---|
142 | |
---|
143 | Duration = (int((time.time()-t0)/3600), int(((time.time()-t0)%3600)/60), int((time.time()-t0)%60)) |
---|
144 | |
---|
145 | print "Finished" |
---|
146 | print "Time taken:" |
---|
147 | print "Hours: %i, Minutes: %i, Seconds: %i" %Duration |
---|
148 | print "Output written to ~", join(rootDir,scenario_name) |
---|