source: production/karratha_2005/run_karratha.py @ 1801

"""Convert from Arcview ASCII DEMs via native netcdf dem format
to native pts netcdf format for use with least_squares fits
"""

import os
import time


from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary
from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
     dem2pts, ferret2sww
from pyvolution.pmesh2domain import pmesh_to_domain_instance 
from caching import cache
import project

#Data preparation
#Convert ASC 2 DEM 2 PTS using source data and store result in source data
demname = project.demname

cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
      dependencies = [demname + '.asc'],
      verbose = True)
      #evaluate = True)

cache(dem2pts, demname, {'verbose': True},
      dependencies = [demname + '.dem'],      
      verbose = True)


#Convert MOST boundary
source_dir = project.boundarydir


#Data
from pyvolution.data_manager import ferret2sww


south = project.south 
north = project.north
west = project.west
east = project.east

#Origin of existing dem (FIXME: Temporary measure)
mesh_origin = (50, 421544.35127423, 7677669.5257159)

ferret2sww(source_dir+project.boundary_basename,
           project.boundary_basename, 
           verbose=True,
           minlat=south-1, maxlat=north+1,
           minlon=west-1, maxlon=east+1,
           origin = mesh_origin,
           zscale = 1,
           fail_on_NaN = False,
           inverted_bathymetry = True),



#Setup domain
mesh = project.meshname + '.msh'
domain = cache(pmesh_to_domain_instance, (mesh, Domain),
               dependencies = [mesh],                     
               verbose = True)               


domain.set_name(project.basename)
domain.set_datadir(project.outputdir)
domain.store = True

print "Number of triangles = ", len(domain)
print 'The extent is ', domain.get_extent()


#Setup IC
domain.set_quantity('stage', 0)
domain.set_quantity('elevation',
                    filename = demname + '.pts',
                    use_cache = True,
                    verbose = True)

#Setup BC
Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)

Br = Reflective_boundary(domain)
domain.set_boundary( {'wall': Br} )


#Run
for t in domain.evolve(yieldstep = 1, finaltime = 100):
    domain.write_time()