source: anuga_work/production/gold_coast_2009/For_DVD/export_results_max.py @ 7306

"""
Generates ascii grids of nominated areas -
Input: sww file from run_perth.py
       boundaries for grids from project.py
Outputs: ascii grids of specified variables
Stored in the 'outputs_dir' folder for respective .sww file

Note:
If producing a grid for the enitre extent cellsize should be greater than 30m
If producing grids for inundation area resolution should be greater than mesh (ie ~22m)
"""

import project, os
import sys
from os.path import join
from anuga.lib.maxasc.maxasc import MaxAsc
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
from anuga.shallow_water.data_manager import sww2dem
from os import sep


directory = project.output_folder

# Directories that contain the sww files of interest found in anuga/outputs:
time_dir1 = 'Event1_HAT'
time_dir2 = 'Event1_MSL'
time_dir3 = 'Event2_HAT'
time_dir4 = 'Event2_MSL'
time_dir5 = 'Event3_HAT'
time_dir6 = 'Event3_MSL'

time_dirs = [time_dir1, time_dir2, time_dir3, time_dir4, time_dir5, time_dir6]   

# sww filename suffixes ie for gold_coast_time_28860_0.sww,
# input 28860 into list, ensuring numbers are in sequential order
time_dirs = [28860, 57720]

# Modify the cellsize variable to set the resolution of the raster you require
# Take into account the mesh size when setting this variable
cellsize = 20 # use 250 if exporting entire model area

#Now set the timestep at which you want the raster generated.
#Either set the actual timestep required or use 'None' to indicate that
#you want the maximum values in the raster over all timesteps
timestep = None

# Set the special areas of interest.  If none, do: area='All'
# Areas identified below are specified by cooridinates in project.py 
area = ['All'] # 'All' means the whole model area

# Specify which properties you would like to export by setting
# one or more key strings from var_equations below
var = ['speed', 'depth']

# Define allowed variable names and associated equations to generate values.
# This would not normally change.

var_equations = {'stage':     'stage',
                 'momentum':  '(xmomentum**2 + ymomentum**2)**0.5',
                 'depth':     'stage-elevation',
                 'speed':     '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)',
                 'elevation': 'elevation' }

######
# Start script, running through variables, area, folder, sww file (determine
# by times). Then creates a maximum asci if there is more than one sww file
######

for which_var in var:
    if which_var not in var_equations:
        print 'Unrecognized variable name: %s' % which_var
        break

    for which_area in area:
        if which_area == 'All':
            easting_min = None
            easting_max = None
            northing_min = None
            northing_max = None
        else:
            try:
                easting_min = eval('project.xmin%s' % which_area)
                easting_max = eval('project.xmax%s' % which_area)
                northing_min = eval('project.ymin%s' % which_area)
                northing_max = eval('project.ymax%s' % which_area)
            except AttributeError:
                print 'Unrecognized area name: %s' % which_area
                break

        for time_dir in time_dirs:
            names = [join(directory, time_dir, project.scenario_name)]
            for time in times:
                names.append(join(directory, time_dir, project.scenario_name)+ '_time_' + str(time) + '_0')
            
            asc_name = []   
                           
            for name in names:
                outname = name + '_' + which_area + '_' + which_var
                quantityname = var_equations[which_var]

                print 'start sww2dem: time_dir=%s' % time_dir
                print 'name=%s' % name
                sww2dem(name, basename_out = outname,
                            quantity = quantityname,
                            timestep = timestep,
                            cellsize = cellsize,      
                            easting_min = easting_min,
                            easting_max = easting_max,
                            northing_min = northing_min,
                            northing_max = northing_max,        
                            reduction = max, 
                            verbose = True,
                            format = 'asc')

                asc_name.append(outname + '.asc')
            
            if len(names) >1:
                maxasc_outname = join(directory, time_dir, project.scenario_name)+'_'+which_area+'_'+which_var+'_max.asc'

                print 'max asc outname ', maxasc_outname
                print 'asc_name ', str(asc_name)
                
                MaxAsc(maxasc_outname, asc_name)