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

Last change on this file since 3769 was 3769, checked in by ole, 18 years ago

Changed references to convert_points_from_latlon_to_utm to the new convert_from_latlan_to_utm as introduced in changeset:3739

File size: 4.8 KB
Line 
1from pylab import *
2from anuga.geospatial_data.geospatial_data import *
3from anuga.pmesh.create_mesh import convert_from_latlon_to_utm
4from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
5from os import sep
6from anuga.utilities.polygon import inside_polygon
7import project
8from anuga.utilities.polygon import plot_polygons, poly_xy, polygon_area
9
10datadir = project.datadir
11ion()
12
13print 'Area of overall polygon: ', polygon_area(project.polyAll)
14
15plot_data = False
16plot_with_poly = True
17
18def read_community_file(filename1, filename2, polygon):
19    from anuga.coordinate_transforms.redfearn import redfearn
20    from anuga.utilities.polygon import inside_polygon
21    fid1 = open(filename1)
22    lines = fid1.readlines()
23    fid1.close()
24    xlocs = []
25    ylocs = []
26    community = []
27    population = []
28    line1 = lines[0]
29    line11 = line1.split(',')
30    fid2 = open(filename2,'w')
31    s = 'Easting, Northing, Community, Population\n'
32    fid2.write(s)
33    for i in range(len(line11)):
34        if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i
35        if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i
36        if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i
37        if line11[i].strip('\n').strip(' ') == 'POPULATION': pop_index = i
38    for line in lines[1:]:
39        fields = line.split(',')
40        lat = float(fields[lat_index])
41        lon = float(fields[lon_index])
42        z, easting, northing = redfearn(lat,lon)
43        v = inside_polygon([easting,northing],polygon)
44        if z == 50 and v == True:
45            xlocs.append(easting)
46            ylocs.append(northing)
47            community = fields[name_index]
48            population = fields[pop_index]
49            s = '%f, %f, %s, %s\n' %(easting, northing, community, population)
50            fid2.write(s)
51            print community, easting, northing
52   
53    return xlocs, ylocs, community, population
54
55def read_file(filename1):
56    fid1 = open(filename1)
57    lines = fid1.readlines()
58    fid1.close()
59    x = []
60    y = []
61    name = []
62    line1 = lines[0]
63    for line in lines[1:]:
64        fields = line.split(',')
65        x.append(float(fields[0]))
66        y.append(float(fields[1]))
67        if len(line) > 2: name.append(fields[2])
68   
69    return x, y, name
70
71if plot_data == True:
72    file1 = datadir+'onslow_onshore_30m_dted.pts'
73    file2 = datadir+'onslow_offshore_points.xya'
74
75    G1 = Geospatial_data(file_name = file1)
76    G2 = Geospatial_data(file_name = file2)
77
78    pts1 = G1.get_data_points(absolute = True)
79    pts2 = G2.get_data_points(absolute = True)
80
81    if plot_with_poly == True:
82        figname = 'onslow_data_poly'
83        x1, y1 = poly_xy(project.polyAll)
84        x2, y2 = poly_xy(project.poly_onslow)
85        x3, y3 = poly_xy(project.poly_coast)
86        x4, y4 = poly_xy(project.poly_region)
87        figure(1)
88        plot(pts1[:,0],pts1[:,1],'g.',
89             pts2[:,0],pts2[:,1],'b.',
90             x1,y1,'r-',x2,y2,'r-',x3,y3,'r-',x4,y4,'r-')
91        xlabel('Eastings')
92        ylabel('Northings')
93        savefig(figname)
94    else:
95        figure(2)
96        plot(pts1[:,0],pts1[:,1],'g.',
97             pts2[:,0],pts2[:,1],'b.')
98        xlabel('Eastings')
99        ylabel('Northings')
100        savefig('onslow_data_extent')
101else:
102   
103    figure(3)
104    figname = 'onslow_polys_test'
105    vec = plot_polygons([project.polyAll,
106                        project.poly_onslow, project.poly_coast],#, project.poly_region],
107                        figname,
108                        verbose = True)
109    figure(4)
110    x1, y1 = poly_xy(project.polyAll)
111    figname = 'onslow_communities'
112    x, y, community, population= \
113           read_community_file(project.community_filename, \
114                               project.community_scenario, project.polyAll)
115    #coastx, coasty = read_file(project.coast_filename)
116    #plot(x1,y1,'r-',x,y,'r+', coastx,coasty,'g.')
117    plot(x1,y1,'r-',x,y,'r+')
118    xlabel('Eastings')
119    ylabel('Northings')
120    savefig(figname)
121    figure(4)
122    figname = 'onslow_gauges'
123    x, y, name = read_file(project.gauge_filename)
124    plot(x1,y1,'r-',x,y,'r+')
125    savefig(figname)
126
127    figure(5)
128    figname = 'onslow_new_boundary_test'
129    gaugex, gaugey, name = read_file(project.gauges50)
130    x50, y50 = poly_xy(project.bounding_poly50)
131    x_onslow, y_onslow = poly_xy(project.poly_onslow)
132    x_coast, y_coast = poly_xy(project.poly_coast)
133    #x25, y25 = poly_xy(project.bounding_poly25)
134    #plot(x50,y50,'r-',x25,y25,'r-',x1,y1,'g-')
135    plot(x50,y50,'r-',x1,y1,'g-', gaugex, gaugey, 'b+', x_onslow, y_onslow,'r-', x_coast, y_coast,'r-')
136    for i, string in enumerate(name):
137        text(gaugex[i],gaugey[i],string,fontsize=7)
138    savefig(figname)
139close('all')
Note: See TracBrowser for help on using the repository browser.