| 1 | from pylab import * |
|---|
| 2 | from anuga.geospatial_data.geospatial_data import * |
|---|
| 3 | from anuga.pmesh.create_mesh import convert_points_from_latlon_to_utm |
|---|
| 4 | from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees |
|---|
| 5 | from os import sep |
|---|
| 6 | from anuga.utilities.polygon import inside_polygon |
|---|
| 7 | import project |
|---|
| 8 | from anuga.utilities.polygon import plot_polygons, poly_xy, polygon_area |
|---|
| 9 | |
|---|
| 10 | datadir = project.datadir |
|---|
| 11 | ion() |
|---|
| 12 | |
|---|
| 13 | print 'Area of overall polygon: ', polygon_area(project.polyAll) |
|---|
| 14 | |
|---|
| 15 | plot_data = False |
|---|
| 16 | plot_with_poly = True |
|---|
| 17 | |
|---|
| 18 | def 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 | |
|---|
| 55 | def 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 | |
|---|
| 71 | if 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') |
|---|
| 101 | else: |
|---|
| 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) |
|---|
| 139 | close('all') |
|---|
