In [1]:
# Allow to import without installing
import sys
sys.path.insert(0, "..")
Let's replicate an OSMnx example by finding all the buildings in Piedmont.
Step one is to find the region which forms the town of Piedmont. As usual, the best way to solve this problem is to play around on www.openstreetmap.org until you find the correct element. In this case, it's a "way" which outlines the town. A bit of browsing around the wiki finds that the important tag is "boundary":"administrative". See http://wiki.openstreetmap.org/w/index.php?title=Tag:boundary%3Dadministrative
In [2]:
import os
filename = os.path.join("..", "..", "..", "Data", "california-latest.db")
In [3]:
import osmdigest.sqlite as sq
db = sq.OSM_SQLite(filename)
In [4]:
all_admin_boundaries = db.search_relation_tags({"boundary": "administrative", "name":"Piedmont"})
all_admin_boundaries
Out[4]:
In [5]:
all_admin_boundaries = [db.complete_relation(rel) for rel in all_admin_boundaries]
In [6]:
import osmdigest.geometry as geometry
import geopandas as gpd
%matplotlib inline
import matplotlib.pyplot as plt
In [11]:
piedmont = geometry.dataframe_from_elements(all_admin_boundaries)
piedmont
Out[11]:
In [12]:
piedmont.plot()
Out[12]:
In [13]:
piedmont.bounds
Out[13]:
In [14]:
import os, sys
try:
os.remove("piedmont.db")
except Exception as ex:
print("{}/{}".format(type(ex), ex), file=sys.stderr)
pass
bounds = piedmont.bounds.ix[0]
sq.extract(db, bounds.minx, bounds.maxx, bounds.miny, bounds.maxy, "piedmont.db")
In [15]:
pied_db = sq.OSM_SQLite("piedmont.db")
frame = geometry.dataframe_from_elements((db.complete_way(way)
for way in pied_db.search_way_tag_keys({"building"})), polygonise=True)
In [16]:
frame[:5]
Out[16]:
In [17]:
frame.plot(figsize=(15,10), color="black")
Out[17]:
In [18]:
# Idea from OSMNX. Much faster to use descartes and patches.
from descartes import PolygonPatch
from matplotlib.collections import PatchCollection
def manual_plot(frame):
bgcolor = "#202020"
fig, ax = plt.subplots(figsize=(15,10))#, facecolor=bgcolor)
ax.set_facecolor(bgcolor)
patches = []
for geo in frame['geometry']:
if geo.geometryType() == "Polygon":
patches.append(PolygonPatch(geo))
elif geo.geometryType() == "MultiPolygon":
for subpolygon in geo:
patches.append(PolygonPatch(subpolygon))
else:
raise Exception(geo.geometryType())
pc = PatchCollection(patches, facecolor="white", edgecolor="white", linewidth=0, alpha=1)
ax.add_collection(pc)
left, bottom, right, top = frame.total_bounds
ax.set_xlim((left, right))
ax.set_ylim((bottom, top))
ax.set_aspect('equal')
manual_plot(frame)
In [19]:
frame = gpd.GeoDataFrame({"geometry":frame.geometry, "osm_id":frame.osm_id})
frame[:5]
Out[19]:
In [20]:
piedmont
Out[20]:
In [21]:
piedmont_reduced = gpd.GeoDataFrame({"geometry":piedmont.geometry})
In [22]:
piedmont_frame = gpd.overlay(piedmont_reduced, frame, how="intersection")
piedmont_frame.crs = {'init':'epsg:4326'}
piedmont_frame = piedmont_frame.to_crs({"init":"epsg:3310"})
In [23]:
manual_plot(piedmont_frame)
In [ ]: