source: trunk/anuga_work/development/Busselton/BP3Hydro/setup_model.py @ 8449

"""
Module to import the project.py 'configuration' file and perform
sanity checks plus a quick check on mesh generation.

Also callable as a stand-alone program, mainly to view the results
of the mesh generation.
"""

import os

from os.path import join, exists
from anuga.geometry.polygon import read_polygon, number_mesh_triangles
from anuga.file.csv_file import load_csv_as_polygons as csv2polygons
import project



#-------------------------------------------------------------------------------
# Sanity checks - check that required files/directories, etc, exist.
#-------------------------------------------------------------------------------

# flag - we check many things and then don't proceed if anything wrong
sanity_error = False               # checked at bottom of this file

# Test that environment variables are defined.
#if os.getenv(project.ENV_INUNDATIONHOME) is None:
#    print "Environment variable '%s' is not set" % project.ENV_INUNDATIONHOME
#    sanity_error = True

#-------------------------------------------------------------------------------
# Directory Structure
#-------------------------------------------------------------------------------

# check folders generated from environment variables.
if not exists(project.home):
    print "Sorry, data directory '%s' doesn't exist" % project.home
    sanity_error = True
       
if not exists(project.anuga_folder):
    print "Sorry, ANUGA directory '%s' doesn't exist" % project.anuga_folder
    sanity_error = True

if not exists(project.topographies_folder):
    print ("Sorry, topo directory '%s' doesn't exist"
           % project.topographies_folder)
    sanity_error = True

if not exists(project.polygons_folder):
    print ("Sorry, polygon directory '%s' doesn't exist"
           % project.polygons_folder)
    sanity_error = True

# if not exists(project.rain_folder):
#    print ("Sorry, rain directory '%s' doesn't exist"
#           % project.rain_folder)
#    sanity_error = True

#if not exists(project.output_folder):
#    print "Sorry, outputs directory '%s' doesn't exist" % project.output_folder
#    sanity_error = True

if not exists(project.hydrographs_folder):
    print "Sorry, hydrographs directory '%s' doesn't exist" % project.hydrographs_folder
    sanity_error = True
        
#if not exists(project.ocean_folder):
#    print "Sorry, ocean directory '%s' doesn't exist" % project.ocean_folder
#    sanity_error = True

#-------------------------------------------------------------------------------
# Determine type of run
#-------------------------------------------------------------------------------
# parameters in this can be changed if desired.
if project.setup == 'trial':
    project.scale_factor = 100      # multiplies the triangles in mesh by 100    
    project.time_thinning = 96
#    project.yieldstep = 240         # data is only stored at 240 sec intervals
elif project.setup == 'basic': 
    project.scale_factor = 4
    project.time_thinning = 12
#    project.yieldstep = 120
elif project.setup == 'final': 
    project.scale_factor = 1        
    project.time_thinning = 4
#    project.yieldstep = 60
else:
    print ("Sorry, you must set the 'setup' variable to one of:"
           '   trial - coarsest mesh, fast\n'
           '   basic - coarse mesh\n'
           '   final - fine mesh, slowest\n'
           '\n'
           "'setup' was set to '%s'" % project.setup)
    sanity_error = True

#-------------------------------------------------------------------------------
# Check for errors detected above.
#-------------------------------------------------------------------------------

if sanity_error:
    msg = 'You must fix the above errors before continuing.'
    raise Exception, msg

#-------------------------------------------------------------------------------
# Reading polygons and creating interior regions
#-------------------------------------------------------------------------------

# Create list of land polygons with initial conditions
##project.land_initial_conditions = []
### if it's a list, then it's a list of land condition filenames
### else it's a string - a single file, multiple land conditions
##if isinstance(project.land_initial_conditions_filename, list):
##    for filename, MSL in project.land_initial_conditions_filename:
##        polygon = read_polygon(join(project.polygons_folder, filename))
##        project.land_initial_conditions.append([polygon, MSL])
##elif isinstance(project.land_initial_conditions_filename, basestring):
##    polygons, MSL = csv2polygons(join(project.polygons_folder,
##                                      project.land_initial_conditions_filename))
##    for i, key in enumerate(polygons):
##        if i%100==0: print i
##        poly = polygons[key]
##        land = float(MSL[key])
##        project.land_initial_conditions.append([poly, land])
##else:
##    msg = ('project.land_initial_conditions_filename must be a list or '
##           'string, got %s'
##           % type(project.land_initial_conditions_filename))
##    raise Exception, msg
### Create list of interior polygons with scaling factor
project.interior_regions = []

if project.PriorityArea_filename is not None:
    polygons, maxareas = csv2polygons(project.PriorityAreas)
    print 'Creating %d internal polygons' % len(polygons)
    #def create_polygon_function(polygons, geo_reference=None):
    project.interior_regions = []
    for i, key in enumerate(polygons):
        if i%100==0: print i
        poly = polygons[key]
        maxarea = float(maxareas[key])
        project.interior_regions.append([poly,
                                         maxarea*project.scale_factor])

for filename, maxarea in project.interior_regions_data:
    polygon = read_polygon(join(project.polygons_folder, filename))
    project.interior_regions.append([polygon,
                                     maxarea*project.scale_factor])

# Initial bounding polygon for data clipping 
project.bounding_polygon = read_polygon(join(project.polygons_folder, 
                                             project.bounding_polygon_filename))
project.bounding_maxarea = project.bounding_polygon_maxarea*project.scale_factor

project.num_boundary_segments = len(project.bounding_polygon) - 1

# Estimate the number of triangles                     
trigs_min = number_mesh_triangles(project.interior_regions,
                                  project.bounding_polygon,
                                  project.bounding_maxarea)

#print 'min estimated number of triangles', trigs_min