source: production/onslow_2006/make_gauges.py @ 3171

""" Make a line of gauges from outside boundary to near centre polygon
for Onslow study to look at MOST wave and compare to ANUGA output
"""
import project
from pylab import plot, xlabel, ylabel, title, ion, axis, savefig, close, text
from utilities.polygon import poly_xy
from Numeric import arange
from 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
print 'MOST clipping: ', MOSTxmin, MOSTxmax, MOSTymin, MOSTymax
# 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
            if ypoint > testy and ypoint <= d0[1]:
            #if ypoint >= testy and ypoint <= MOSTymax:
                if xpoint > MOSTxmin and xpoint < MOSTxmax:
                    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
            if ypoint > testy and ypoint < MOSTymax: #testy
                if xpoint > MOSTxmin and xpoint < MOSTxmax:
                    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')