source: anuga_work/production/busselton/build_boundary.py @ 5575

from anuga.utilities.numerical_tools import ensure_numeric
from Scientific.IO.NetCDF import NetCDFFile
from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
    compress,zeros,fabs,allclose,ones
from anuga.utilities.polygon import inside_polygon,read_polygon
from os import sep
from time import localtime, strftime, gmtime
from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
import os
import project


######################################
# Create urs boundary from mux2files #
######################################

#print project.muxhome
dir=os.path.join(project.muxhome,'mux')
#print dir 
prefix=os.path.join(dir,'Java-00')
suffix='-z.grd'

urs_filenames={}
for filename in os.listdir(dir):
    basename = filename[:-7]

    urs_filenames[basename] = True


# Make list of basenames  
urs_filenames = [os.path.join(dir, basename) for basename in urs_filenames.keys()]

# AS per David Burbidge email on friday 4th July the mag 9.3 event
# has 1m worth of slip on each sub fault therefore mutliple each unit
# soucre by the slip (10.4544) and sum the 44 time series together to
# get the time series for this event at the points on your boundary. 

weights=10.4544*ones(len(urs_filenames),Float)

scenario_name=project.scenario_name
order_filename=os.path.join(project.boundaries_dir, 'thinned_bounding_order.txt')

print 'reading', order_filename
# Create ordered sts file
print 'creating sts file'

urs2sts(urs_filenames,
        basename_out=os.path.join(project.boundaries_dir,scenario_name),
        ordering_filename=order_filename,
        weights=weights,
        mean_stage=project.tide,
        verbose=True)

# Read in boundary from ordered sts file
urs_boundary_polygon=create_sts_boundary(scenario_name)