source: anuga_work/production/dampier_2006/get_timeseries_bf.py @ 4631

"""
Generate images of "gauges" for production run

Inputs:

production dirs: dictionary of production directories with a
                 association to that simulation run, eg high tide,
                 low tide and MSL.
                    
Outputs:

* figures used for report stored in the report_figure directory

"""
from os import getcwd, sep, altsep, mkdir, access, F_OK
import project
from anuga.abstract_2d_finite_volumes.util import sww2timeseries
from anuga.shallow_water import Domain
from anuga.pmesh.mesh_interface import create_mesh_from_regions 
from anuga.shallow_water import File_boundary
#from anuga.fit_interpolate import Interpolation_function
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher


start_screen_catcher(project.temp_dir_name, print_to_screen=True)

# Create report directory
reportdir = getcwd()+sep+'report'+sep
if access(reportdir,F_OK) == 0:
    mkdir (reportdir)
    
# sww file created from URS boundary
production_dirs = {#'20070206_005106_run': '10000 yr wave height'}
                   'boundaries': 'URS boundary condition'}
                   
meshes_dir_name = project.meshes_dir_name+'.msh'
print 'create mesh' 
create_mesh_from_regions(project.poly_all,
                             boundary_tags=None,
                             maximum_triangle_area=500000,
                             interior_regions=None,
                             filename=meshes_dir_name,
                             use_cache=False,
                             verbose=True)
                             
domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
#print"domain id", id(domain)

print domain.statistics()
boundaries_dir_name=project.boundaries_dir_name

domain.set_quantity('stage', 2.4)
print 'Reading Boundary file', boundaries_dir_name
Bf = File_boundary(project.boundaries_dir_name3 + '.sww',
                  domain, time_thinning=1, use_cache=True, verbose=True)
#print "domain id", id(domain)
print "Bf id", id(Bf)

print "dir,bf", dir(Bf)
#print 'finished reading boundary file', Bf.get_midpoint_coords()
print 'finished reading boundary file', Bf.midpoint_coordinates.shape[0],Bf.midpoint_coordinates[0]

# interpolation_points[0]

#number_of_points = Bf.midpoint_coordinates.shape[0]
number_of_points = (50,200,324)

#number_of_points = Bf.interpolation_points.shape[0]
print "number_of_points", number_of_points
dir(Bf.F)
#for id in range(number_of_points)[10:]:
for id in number_of_points:
    print "id", id
    x,y=Bf.midpoint_coordinates[id]
#    x,y=Bf.interpolation_points[id]
    print"x,y",x,y
    t_boundary = Bf.F.precomputed_values["stage"][:,id]
    print "t_boundary", t_boundary

#t_file = 


import sys   
sys.exit()



is_parallel = False
if is_parallel == True:
    nodes = 2
    
# Generate figures
swwfiles = {}

if is_parallel == True:

    for i in range(nodes):
        print 'Sending node %d of %d' %(i,nodes)
        swwfiles = {}
        #reportname = report_name + '_%s' %(i)
        for label_id in production_dirs.keys():
            if label_id == 'boundaries':
                file_loc = project.boundaries_in_dir
                swwfile = file_loc + sep + '20061102_221245_build' + sep + 'dampier.sww'
                print 'hi'
            #else:
            #file_loc = project.output_dir + label_id + sep
            #sww_extra = '_P%s_%s' %(i,nodes)
            #swwfile = file_loc + project.scenario_name + sww_extra + '.sww'
            swwfiles[swwfile] = label_id

        texname, elev_output = sww2timeseries(swwfiles,
                                              project.gauges_dir_name,
                                              production_dirs,
                                              report = False,
                                              #reportname = reportname,
                                              plot_quantity = ['stage', 'momentum'],
                                              generate_fig = False,
                                              surface = False,
                                              time_min = None,
                                              time_max = None,
                                              title_on = False,
                                              verbose = True)
            
            #latex_output.append(texname)
    
else:
    
    for label_id in production_dirs.keys():
        
        if label_id == 'boundaries':
                file_loc = project.boundaries_in_dir
                swwfile = project.boundaries_dir_name3 + '.sww'
        else:
            file_loc = project.output_dir + label_id + sep
            swwfile = file_loc + project.scenario_name + '.sww'
        swwfiles[swwfile] = label_id
        
    texname, elev_output = sww2timeseries(swwfiles,
                                          project.gauges_dir_name_test,
                                          production_dirs,
                                          report = False,
                                          #reportname = report_name,
                                          plot_quantity = ['stage', 'momentum'],
                                          surface = False,
                                          time_min = None,
                                          time_max = None,
                                          title_on = False,
                                          verbose = True)

    fid = open(filename)
    lines = fid.readlines()
    fid.close()
    polygon = []
    for line in lines:
        fields = line.split(split)
        polygon.append( [float(fields[0]), float(fields[1])] )