source: anuga_work/production/onslow_2006/make_gauges.py @ 3650

""" Make a line of gauges from near the boundary to near centre polygon
for Onslow study to look at MOST wave and compare to ANUGA output
Output from here is used as the gauge file input in MOST_timeseries.py
"""
import project
from pylab import plot, xlabel, ylabel, title, ion, axis, savefig, close, text
from anuga.utilities.polygon import poly_xy, is_inside_polygon
from Numeric import arange
from anuga.pmesh.create_mesh import convert_points_from_latlon_to_utm

x_bound, y_bound = poly_xy(project.polyAll)

MOST_south = project.south
MOST_north = project.north
MOST_east = project.east
MOST_west = project.west
p0 = [MOST_south, MOST_west]
p1 = [MOST_south, MOST_east]
p2 = [MOST_north, MOST_east]
p3 = [MOST_north, MOST_west]
MOSTpoly, zone = convert_points_from_latlon_to_utm([p0, p1, p2, p3])

MOSTymax = 0
MOSTymin = 7800000.0
MOSTxmax = 0
MOSTxmin = 7800000.0
for i, point in enumerate(MOSTpoly):
    x = point[0]
    y = point[1]
    if y > MOSTymax: MOSTymax = y
    if y < MOSTymin: MOSTymin = y
    if x > MOSTxmax: MOSTxmax = x
    if x < MOSTxmin: MOSTxmin = x

MOSTymax = MOSTymax-3000.0
# nominated point in centre region for Onslow
x2 = 305000.0
y2 = 7605000.0
testy = 7625000.0

d0 = project.d0
d1 = project.d1
d2 = project.d2
d3 = project.d3

d = [d0, d1, d2, d3]
ion()
fid = open(project.gauge_comparison,'w')
s = 'Easting,Northing,Name,Elevation \n'
fid.write(s)
count = 0
miny = 7580000.0
maxy = 7800000.0
xplot = []
yplot = []
names = []
for i, point in enumerate(d):
    x = point[0]
    y = point[1]
    m = (y-y2)/(x-x2)
    c = y2-m*x2
    if m > 0:
        plotx = arange(MOSTxmin,MOSTxmax,1000.0)
        for j, xpoint in enumerate(plotx):
            ypoint = m*xpoint+c
            test = is_inside_polygon([xpoint,ypoint],project.polyAll,closed=True)
            if test == True:
                if ypoint > testy and ypoint < MOSTymax:
                    xplot.append(xpoint)
                    yplot.append(ypoint)
                    thisname = 'Point%d' %count
                    names.append(thisname)
                    count += 1
                    s = '%.2f,%.2f,%s,%.2f \n' %(xpoint,ypoint,thisname,0.0)
                    fid.write(s)
    else:
        plotx = arange(MOSTxmin,MOSTxmax,5000.0)
        for j, xpoint in enumerate(plotx):
            ypoint = m*xpoint+c
            test = is_inside_polygon([xpoint,ypoint],project.polyAll,closed=True)
            if test == True:
                if ypoint > testy and ypoint < MOSTymax:
                    xplot.append(xpoint)
                    yplot.append(ypoint)
                    thisname = 'Point%d' %count
                    names.append(thisname)
                    count += 1
                    s = '%.2f,%.2f,%s,%.2f \n' %(xpoint,ypoint,thisname,0.0)
                    fid.write(s)

    plot(x_bound, y_bound, xplot, yplot,'b+')
    for i, string in enumerate(names):
        if xplot[i] > 240000 and xplot[i] < 340000:
            text(xplot[i],yplot[i],string,fontsize=7)

axis([240000,340000, miny, maxy])
savefig('boundarycomparison')
close('all')