source: production/pt_hedland_2006/export_results.py @ 3330

import project, os
import sys

from pyvolution.data_manager import sww2dem
from pyvolution.ermapper_grids import read_ermapper_grid
from pyvolution.util import Screen_Catcher
from os import sep

time_dir = '20060706_235036' #test with coarse grid
#time_dir = '20060707_001859' #MSL DLI data
#time_dir = '20060707_003301' #HAT DLI data
#time_dir = '20060707_003424' #LAT DLI data
directory = project.outputdir
name = directory + time_dir +sep + 'source'


#normal screen output is stored in 'outname'
screen_output_name = directory + time_dir + sep + "export_output_elev.txt"
screen_error_name = directory + time_dir + sep + "export_error_elev.txt"

#print 'hello', name

#used to catch screen output to file
sys.stdout = Screen_Catcher(screen_output_name)
sys.stderr = Screen_Catcher(screen_error_name)

print 'output dir:', name
#print 'Which variable do you want to export?'
#which_var = int(raw_input('Stage = 0, Absolute Momentum = 1, Depth = 2, Speed = 3  '  ))
which_var = 4
#sys.stderr.write(sys.stdout.data)
if which_var == 0:  # Stage
    outname = name + '_stage'
    quantityname = 'stage'

if which_var == 1:  # Absolute Momentum
    outname = name + '_momentum'
    quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  #Absolute momentum

if which_var == 2:  # Depth
    outname = name + '_depth'
    quantityname = 'stage-elevation'  #Depth

if which_var == 3:  # Speed
    outname = name + '_speed'
    #quantityname = '((xmomentum/(stage-elevation))**2 + (ymomentum/(stage-elevation))**2)**0.5'  #Speed
    #quantityname = 'xmomentum/(stage-elevation)'  #Speed
    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-30)'  #Speed

if which_var == 4:  # Elevation
    outname = name + '_elevation'
    quantityname = 'elevation'  #Elevation

print 'start sww2dem'
#sys.stderr.write(sys.stdout.data)
sww2dem(name, basename_out = outname,
            quantity = quantityname,
            cellsize = 20,       # Trevor would like this at 25
            # define region for viz purposes
            easting_min = project.e_min_area,
            easting_max = project.e_max_area,
            northing_min = project.n_min_area,
            northing_max = project.n_max_area,        
            reduction = max, #this is because we want max quantityname
            verbose = True,
            format = 'asc')

#sys.stderr.write(sys.stdout.data)

#Check

#data = read_ermapper_grid(name)
#print 'Values from %s are in [%f, %f]' %(name, min(data.flat), max(data.flat))