diff --git a/map-state.py b/map-state.py new file mode 100644 index 0000000..4b1eeba --- /dev/null +++ b/map-state.py @@ -0,0 +1,211 @@ +import overpy +import numpy +import pickle +from matplotlib import pyplot +from shapely.geometry import LineString +from shapely.geometry import shape +from shapely.geometry import box +from shapely.ops import polygonize_full +from descartes import PolygonPatch +from fiona import collection + +features = collection("~/Desktop/gshhg-shp-2/GSHHS_shp/f/GSHHS_f_L1.shp") +northAmericaPolygon = shape(features[3]['geometry']) + +api = overpy.Overpass() + +stateRelation = api.query("rel(165475);(._;>;);out;") +fnwrRelation = api.query("rel(3947664);(._;>;);out;") +motorwaysInBoundingBox = api.query("way(32.120,-125.222,42.212,-113.928)[highway=motorway];(._;>);out;") +trunkInBoundingBox = api.query("way(32.120,-125.222,42.212,-113.928)[highway=trunk];(._;>);out;") + +stateLineStrings = [] +for way in stateRelation.ways: + lineString = [] + for node in way.nodes: + lineString.append((node.lon,node.lat)) + stateLineStrings.append(LineString(lineString)) + +polygons, dangles, cuts, invalids = polygonize_full(stateLineStrings) +statePolygonWithOcean = polygons.geoms[3] + +statePolygon = statePolygonWithOcean.intersection(northAmericaPolygon) + +fnwrLineStrings = [] +for way in fnwrRelation.ways: + lineString = [] + for node in way.nodes: + lineString.append((node.lon,node.lat)) + fnwrLineStrings.append(LineString(lineString)) + +Fpolygons, dangles, cuts, invalids = polygonize_full(fnwrLineStrings) + +motorwayLineStrings = [] +for way in motorwaysInBoundingBox.ways: + line = [] + for node in way.nodes: + line.append((node.lon,node.lat)) + wayLineString = LineString(line) + if statePolygonWithOcean.contains(wayLineString): motorwayLineStrings.append(wayLineString) + +for way in trunkInBoundingBox.ways: + line = [] + for node in way.nodes: + line.append((node.lon,node.lat)) + wayLineString = LineString(line) + if statePolygonWithOcean.contains(wayLineString): motorwayLineStrings.append(wayLineString) + +islandPolygons = [] +for i in range(0,len(polygons.geoms)): + if (i != 3) and (i != 4): + islandPolygons.append(polygons.geoms[i]) + +for geom in Fpolygons.geoms: + islandPolygons.append(geom) + +fig = pyplot.figure(figsize=(100,100)) +ax = fig.add_subplot(111) + +for line in motorwayLineStrings: + x, y = line.xy + ax.plot(x, y, color='#000000', linewidth=1, zorder=2) + +patch = PolygonPatch(statePolygon, fc='#FFFFFF', ec='#000000', zorder=1) +ax.add_patch(patch) + +for polygon in islandPolygons: + patch = PolygonPatch(polygon, fc='#FFFFFF', ec='#000000', zorder=1) + ax.add_patch(patch) + +fig.savefig('test2.png') + +##if we want more roads: +hgvInBoundingBox = api.query("way(32.120,-125.222,42.212,-113.928)[hgv=designated];(._;>);out;") +motorwayLineStrings = [] +for way in hgvInBoundingBox.ways: + line = [] + for node in way.nodes: + line.append((node.lon,node.lat)) + wayLineString = LineString(line) + if statePolygonWithOcean.contains(wayLineString): motorwayLineStrings.append(wayLineString) + + +##pickle +output = open('data.pkl', 'wb') +pickle.dump(islandPolygons,output) +pickle.dump(statePolygon,output) +pickle.dump(motorwayLineStrings,output) +output.close() + +pkl_file = open('data.pkl', 'rb') +islandPolygons = pickle.load(pkl_file) +statePolygon = pickle.load(pkl_file) +motorwayLineStrings = pickle.load(pkl_file) +pkl_file.close() + + +##multiprocessing +from multiprocessing import Pool + +def buildLineString(way): + lineString = [] + for node in way.nodes: + lineString.append((node.lon,node.lat)) + return lineString + +p = Pool(8) + +for way in stateRelation.ways: + stateLineStrings = p.map(buildLineString,stateRelation.ways) + +features = collection("~/Desktop/gshhg-shp-2/GSHHS_shp/f/GSHHS_f_L1.shp") +northAmericaPolygon = shape(features[3]['geometry']) + +api = overpy.Overpass() + +stateRelation = api.query("rel(165475);(._;>;);out;") +fnwrRelation = api.query("rel(3947664);(._;>;);out;") +motorwaysInBoundingBox = api.query("way(32.120,-125.222,42.212,-113.928)[highway=motorway];(._;>);out;") + +stateLineStrings = [] +for way in stateRelation.ways: + lineString = [] + for node in way.nodes: + lineString.append((node.lon,node.lat)) + stateLineStrings.append(LineString(lineString)) + +## debug info +from geopandas import GeoSeries +from geopandas import GeoDataFrame + +for i in range(0,len(islandPolygons)): + ax.text(islandPolygons[i].centroid.x,islandPolygons[i].centroid.y,i,color = 'k', weight = 'bold') + +ax.plot([-124.482003,-114.1307816],[32.5295236,42.009499]) + +from shapely.geometry import MultiLineString +test = MultiLineString(highwayLineStrings) +patch = PolygonPatch(countyPolygon, fc='#6699cc', ec='#6699cc', alpha=0.5, zorder=2) +ax.add_patch(patch) + +object.__dict__ +result.nodes +result.nodes.[0].lat +result.ways +result.ways[0]._node_ids[0] +result.relations +result.relations[0].members[0].ref + +len(streetsInBB.ways[0].get_nodes(resolve_missing=True)) + +for way in aboveIslands: + gray = gray + .025 + for node in way.nodes: + ax.scatter(node.lon, node.lat, color=[gray,gray,gray], s=100, zorder=1) + +states = [shapely.geometry.shape(f['geometry']) for f in features] + +http://overpass.osm.rambler.ru/cgi/interpreter?data=%5Bout:json%5D;relation(396479);out; + +Farallon National Wildlife Refuge (3947664) +Santa Catalina Island (237602) +San Clemente Island (237603) +Anacapa Island (3635899) +Santa Cruz Island (237600) +Santa Rosa Island (237601) +Way: Santa Barbara Island (40501085) +Way: Sutil Island (40501070) +Way: San Nicolas Island (40500976) +Way: San Miguel Island (40500912) + +islandRelationNumbers = [3947664,237602,237603,3635899,237600,237601] +islandWayNumbers = [40501085,40501070,40500976,40500912] + +islandRelations = [] +for relation in islandRelationNumbers: + response = api.query("rel(" + relation + ");(._;>;);out;") + islandRelations.append(response) + +islandWays = [] +for ways in islandWayNumbers: + response = api.query("way(" + way + ");(._;>;);out;") + islandWays.append(response) + +islandPolygons = [] +for relation in islandRelations: + islandLineStrings = [] + for way in relation.ways: + lineString = [] + for node in way.nodes: + lineString.append((node.lon,node.lat)) + islandLineStrings.append(LineString(lineString)) + polygons, dangles, cuts, invalids = polygonize_full(stateLineStrings) + islandPolygons.append(polygons.geoms[0]) + +for way in islandWays: + lineString = [] + for node in way.nodes: + lineString.append((node.lon,node.lat)) + islandLineStrings.append(LineString(lineString)) + polygons, dangles, cuts, invalids = polygonize_full(stateLineStrings) + islandPolygons.append(polygons.geoms[0])