source: anuga_work/production/onslow_2006/plot_data_extent.py @ 4354

from pylab import *
from anuga.geospatial_data.geospatial_data import *
from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
from os import sep
from anuga.utilities.polygon import inside_polygon
import project
from anuga.utilities.polygon import plot_polygons, poly_xy, polygon_area

datadir = project.datadir
ion()

print 'Area of overall polygon: ', polygon_area(project.polyAll)

plot_data = False
plot_with_poly = True

def read_community_file(filename1, filename2, polygon):
    from anuga.coordinate_transforms.redfearn import redfearn
    from anuga.utilities.polygon import inside_polygon
    fid1 = open(filename1)
    lines = fid1.readlines()
    fid1.close()
    xlocs = []
    ylocs = []
    community = []
    population = []
    line1 = lines[0]
    line11 = line1.split(',')
    fid2 = open(filename2,'w')
    s = 'Easting, Northing, Community, Population\n'
    fid2.write(s)
    for i in range(len(line11)):
        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
        if line11[i].strip('\n').strip(' ') == 'POPULATION': pop_index = i
    for line in lines[1:]:
        fields = line.split(',')
        lat = float(fields[lat_index])
        lon = float(fields[lon_index])
        z, easting, northing = redfearn(lat,lon)
        v = inside_polygon([easting,northing],polygon)
        if z == 50 and v == True:
            xlocs.append(easting)
            ylocs.append(northing)
            community = fields[name_index]
            population = fields[pop_index]
            s = '%f, %f, %s, %s\n' %(easting, northing, community, population)
            fid2.write(s)
            print community, easting, northing
    
    return xlocs, ylocs, community, population

def read_file(filename1):
    fid1 = open(filename1)
    lines = fid1.readlines()
    fid1.close()
    x = []
    y = []
    name = []
    line1 = lines[0]
    for line in lines[1:]:
        fields = line.split(',')
        x.append(float(fields[0]))
        y.append(float(fields[1]))
        if len(line) > 2: name.append(fields[2])
    
    return x, y, name

if plot_data == True:
    file1 = datadir+'onslow_onshore_30m_dted.pts'
    file2 = datadir+'onslow_offshore_points.xya'

    G1 = Geospatial_data(file_name = file1)
    G2 = Geospatial_data(file_name = file2)

    pts1 = G1.get_data_points(absolute = True)
    pts2 = G2.get_data_points(absolute = True)

    if plot_with_poly == True:
        figname = 'onslow_data_poly'
        x1, y1 = poly_xy(project.polyAll)
        x2, y2 = poly_xy(project.poly_onslow)
        x3, y3 = poly_xy(project.poly_coast)
        x4, y4 = poly_xy(project.poly_region)
        figure(1)
        plot(pts1[:,0],pts1[:,1],'g.',
             pts2[:,0],pts2[:,1],'b.',
             x1,y1,'r-',x2,y2,'r-',x3,y3,'r-',x4,y4,'r-')
        xlabel('Eastings')
        ylabel('Northings')
        savefig(figname)
    else:
        figure(2)
        plot(pts1[:,0],pts1[:,1],'g.',
             pts2[:,0],pts2[:,1],'b.')
        xlabel('Eastings')
        ylabel('Northings')
        savefig('onslow_data_extent')
else:
    
    figure(3)
    figname = 'onslow_polys_test'
    vec = plot_polygons([project.polyAll,
                        project.poly_onslow, project.poly_coast],#, project.poly_region],
                        figname,
                        verbose = True)
    figure(4)
    x1, y1 = poly_xy(project.polyAll)
    figname = 'onslow_communities'
    x, y, community, population= \
           read_community_file(project.community_filename, \
                               project.community_scenario, project.polyAll)
    #coastx, coasty = read_file(project.coast_filename)
    #plot(x1,y1,'r-',x,y,'r+', coastx,coasty,'g.')
    plot(x1,y1,'r-',x,y,'r+')
    xlabel('Eastings')
    ylabel('Northings')
    savefig(figname)
    figure(4)
    figname = 'onslow_gauges'
    x, y, name = read_file(project.gauge_filename)
    plot(x1,y1,'r-',x,y,'r+')
    savefig(figname)

    figure(5)
    figname = 'onslow_new_boundary_test'
    gaugex, gaugey, name = read_file(project.gauges50)
    x50, y50 = poly_xy(project.bounding_poly50)
    x_onslow, y_onslow = poly_xy(project.poly_onslow)
    x_coast, y_coast = poly_xy(project.poly_coast)
    #x25, y25 = poly_xy(project.bounding_poly25)
    #plot(x50,y50,'r-',x25,y25,'r-',x1,y1,'g-')
    plot(x50,y50,'r-',x1,y1,'g-', gaugex, gaugey, 'b+', x_onslow, y_onslow,'r-', x_coast, y_coast,'r-')
    for i, string in enumerate(name):
        text(gaugex[i],gaugey[i],string,fontsize=7)
    savefig(figname)
close('all')