[2860] | 1 | from pylab import * |
---|
| 2 | from geospatial_data import * |
---|
| 3 | from pmesh.create_mesh import convert_points_from_latlon_to_utm |
---|
| 4 | from coordinate_transforms.redfearn import degminsec2decimal_degrees |
---|
| 5 | from os import sep |
---|
| 6 | from utilities.polygon import inside_polygon |
---|
| 7 | import project |
---|
[2920] | 8 | from utilities.polygon import plot_polygons, poly_xy |
---|
[2860] | 9 | |
---|
| 10 | datadir = project.datadir |
---|
| 11 | ion() |
---|
| 12 | |
---|
[2922] | 13 | plot_data = False |
---|
[2920] | 14 | plot_with_poly = True |
---|
[2860] | 15 | |
---|
[2922] | 16 | def read_community_file(filename1, filename2, polygon): |
---|
| 17 | from coordinate_transforms.redfearn import redfearn |
---|
| 18 | from utilities.polygon import inside_polygon |
---|
| 19 | fid1 = open(filename1) |
---|
| 20 | lines = fid1.readlines() |
---|
| 21 | fid1.close() |
---|
| 22 | xlocs = [] |
---|
| 23 | ylocs = [] |
---|
| 24 | community = [] |
---|
| 25 | population = [] |
---|
| 26 | line1 = lines[0] |
---|
| 27 | line11 = line1.split(',') |
---|
| 28 | fid2 = open(filename2,'w') |
---|
| 29 | s = 'Easting, Northing, Community, Population\n' |
---|
| 30 | fid2.write(s) |
---|
| 31 | for i in range(len(line11)): |
---|
| 32 | if line11[i].strip('\n').strip(' ') == 'LATITUDE': lat_index = i |
---|
| 33 | if line11[i].strip('\n').strip(' ') == 'LONGITUDE': lon_index = i |
---|
| 34 | if line11[i].strip('\n').strip(' ') == 'COMMUNITY': name_index = i |
---|
| 35 | if line11[i].strip('\n').strip(' ') == 'POPULATION': pop_index = i |
---|
| 36 | for line in lines[1:]: |
---|
| 37 | fields = line.split(',') |
---|
| 38 | lat = float(fields[lat_index]) |
---|
| 39 | lon = float(fields[lon_index]) |
---|
| 40 | z, easting, northing = redfearn(lat,lon) |
---|
| 41 | v = inside_polygon([easting,northing],polygon) |
---|
[2947] | 42 | if z == 50 and v == True: |
---|
[2922] | 43 | xlocs.append(easting) |
---|
| 44 | ylocs.append(northing) |
---|
| 45 | community = fields[name_index] |
---|
| 46 | population = fields[pop_index] |
---|
| 47 | s = '%f, %f, %s, %s\n' %(easting, northing, community, population) |
---|
| 48 | fid2.write(s) |
---|
| 49 | print community, easting, northing |
---|
| 50 | |
---|
| 51 | return xlocs, ylocs, community, population |
---|
| 52 | |
---|
| 53 | def read_file(filename1): |
---|
| 54 | fid1 = open(filename1) |
---|
| 55 | lines = fid1.readlines() |
---|
| 56 | fid1.close() |
---|
| 57 | x = [] |
---|
| 58 | y = [] |
---|
| 59 | line1 = lines[0] |
---|
| 60 | for line in lines[1:]: |
---|
| 61 | fields = line.split(',') |
---|
| 62 | x.append(float(fields[0])) |
---|
| 63 | y.append(float(fields[1])) |
---|
| 64 | |
---|
| 65 | return x, y |
---|
| 66 | |
---|
[2908] | 67 | if plot_data == True: |
---|
| 68 | file1 = datadir+'onslow_onshore_30m_dted.pts' |
---|
| 69 | file2 = datadir+'onslow_offshore_points.xya' |
---|
[2860] | 70 | |
---|
[2908] | 71 | G1 = Geospatial_data(file_name = file1) |
---|
| 72 | G2 = Geospatial_data(file_name = file2) |
---|
[2860] | 73 | |
---|
[2908] | 74 | pts1 = G1.get_data_points(absolute = True) |
---|
| 75 | pts2 = G2.get_data_points(absolute = True) |
---|
[2860] | 76 | |
---|
[2908] | 77 | if plot_with_poly == True: |
---|
| 78 | figname = 'onslow_data_poly' |
---|
[2920] | 79 | x1, y1 = poly_xy(project.polyAll) |
---|
| 80 | x2, y2 = poly_xy(project.poly_onslow) |
---|
| 81 | x3, y3 = poly_xy(project.poly_coast) |
---|
| 82 | x4, y4 = poly_xy(project.poly_region) |
---|
[2908] | 83 | figure(1) |
---|
| 84 | plot(pts1[:,0],pts1[:,1],'g.', |
---|
[2920] | 85 | pts2[:,0],pts2[:,1],'b.', |
---|
| 86 | x1,y1,'r-',x2,y2,'r-',x3,y3,'r-',x4,y4,'r-') |
---|
| 87 | xlabel('Eastings') |
---|
| 88 | ylabel('Northings') |
---|
[2908] | 89 | savefig(figname) |
---|
[2860] | 90 | else: |
---|
[2908] | 91 | figure(2) |
---|
| 92 | plot(pts1[:,0],pts1[:,1],'g.', |
---|
| 93 | pts2[:,0],pts2[:,1],'b.') |
---|
[2920] | 94 | xlabel('Eastings') |
---|
| 95 | ylabel('Northings') |
---|
[2908] | 96 | savefig('onslow_data_extent') |
---|
| 97 | else: |
---|
| 98 | |
---|
| 99 | figure(3) |
---|
| 100 | figname = 'onslow_polys' |
---|
[2920] | 101 | vec = plot_polygons([project.polyAll, |
---|
| 102 | project.poly_onslow, project.poly_coast, project.poly_region], |
---|
[2914] | 103 | figname, |
---|
| 104 | verbose = True) |
---|
[2922] | 105 | figure(4) |
---|
| 106 | x1, y1 = poly_xy(project.polyAll) |
---|
| 107 | figname = 'onslow_communities' |
---|
| 108 | x, y, community, population= \ |
---|
| 109 | read_community_file(project.community_filename, \ |
---|
| 110 | project.community_scenario, project.polyAll) |
---|
| 111 | #coastx, coasty = read_file(project.coast_filename) |
---|
| 112 | #plot(x1,y1,'r-',x,y,'r+', coastx,coasty,'g.') |
---|
| 113 | plot(x1,y1,'r-',x,y,'r+') |
---|
| 114 | xlabel('Eastings') |
---|
| 115 | ylabel('Northings') |
---|
| 116 | savefig(figname) |
---|
[2860] | 117 | |
---|
[2908] | 118 | close('all') |
---|