""" 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 utilities.polygon import poly_xy, is_inside_polygon 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 # 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')