In [1]:
import sys, os
sys.path.insert(0, os.path.join("..", ".."))
In [2]:
%matplotlib inline
import matplotlib.pyplot as plt
import geopandas as gpd
import os
import descartes
import open_cp.plot
import matplotlib.collections
In [3]:
datadir = os.path.join("/media", "disk", "Data")
areas = gpd.GeoDataFrame.from_file(os.path.join(datadir, "Chicago_Areas.geojson"))
These are the canonical source for USA data.
The downloads are all shapefiles, so having a look at them in QGis is a good idea...
There are other files you can obtain as well, but the "roads" one seems good enough.
In [5]:
tiger_path = os.path.join("/media", "disk", "TIGER Data")
In [5]:
filename = os.path.join(tiger_path, "tl_2016_17031_roads")
tiger_frame = gpd.GeoDataFrame.from_file(filename)
tiger_frame.head()
Out[5]:
In [6]:
fig, ax = plt.subplots(figsize=(12,12))
lc = matplotlib.collections.LineCollection(open_cp.plot.lines_from_geometry(tiger_frame.geometry),
color="black", linewidth=0.5)
ax.add_collection(lc)
patches = open_cp.plot.patches_from_geometry(areas.geometry,
fc="blue", ec="none", alpha=0.5, zorder=5)
ax.add_collection(matplotlib.collections.PatchCollection(patches))
ax.set_aspect(1.0)
xmin, ymin, xmax, ymax = areas.total_bounds
xd, yd = xmax - xmin, ymax - ymin
ax.set(xlim=[xmin - xd/10, xmax + xd/10], ylim=[ymin - yd/10, ymax + yd/10])
None
In [6]:
filename = os.path.join(tiger_path, "tl_2016_17031_edges")
tiger_frame = gpd.GeoDataFrame.from_file(filename)
tiger_frame.head()
Out[6]:
In [7]:
fig, ax = plt.subplots(figsize=(12,12))
lc = matplotlib.collections.LineCollection(open_cp.plot.lines_from_geometry(tiger_frame.geometry),
color="black", linewidth=0.5)
ax.add_collection(lc)
patches = open_cp.plot.patches_from_geometry(areas.geometry,
fc="blue", ec="none", alpha=0.5, zorder=5)
ax.add_collection(matplotlib.collections.PatchCollection(patches))
ax.set_aspect(1.0)
xmin, ymin, xmax, ymax = areas.total_bounds
xd, yd = xmax - xmin, ymax - ymin
ax.set(xlim=[xmin - xd/10, xmax + xd/10], ylim=[ymin - yd/10, ymax + yd/10])
None
The easiest option appears to be:
In [7]:
ordnancedir = os.path.join("/media", "disk", "OS_OpenMap")
openroadsdir = os.path.join(ordnancedir, "oproad_essh_gz", "data")
file = "SE_RoadLink.shp"
In [8]:
frame = gpd.GeoDataFrame.from_file(os.path.join(openroadsdir, file))
frame.head()
Out[8]:
In [9]:
fig, ax = plt.subplots(figsize=(12,12))
lines = open_cp.plot.lines_from_geometry(frame.geometry)
# We get `z` coords, but always 0
for line in lines:
for coord in line:
assert coord[2] == 0
lines = [ [coord[:2] for coord in line] for line in lines ]
lc = matplotlib.collections.LineCollection(lines, color="black", linewidth=0.5)
ax.add_collection(lc)
ax.set_aspect(1.0)
xmin, ymin, xmax, ymax = frame.total_bounds
xd, yd = xmax - xmin, ymax - ymin
ax.set(xlim=[xmin - xd/20, xmax + xd/20], ylim=[ymin - yd/20, ymax + yd/20])
None
For example, my office is at (53.80417, -1.55494).
Now perform the calculation ourself:
transform method of pyproj. CreditFor OS VectorMap raster images, each image 4000x4000 corresponding to 10000m square.
For OS OpenMap Local we have 4 further squares per image; our residual is in the SE corner.
In [10]:
import pyproj
bng = pyproj.Proj(init="epsg:27700")
wgs84 = pyproj.Proj(init="epsg:4326")
pyproj.transform(wgs84, bng, -1.55494, 53.80417)
Out[10]:
In [11]:
import PIL.Image
In [12]:
vectormap_raster = os.path.join(ordnancedir, "OS VectorMap District (Full Colour Raster) SE", "data")
image1 = PIL.Image.open(os.path.join(vectormap_raster, "SE23.tif"))
In [13]:
x, y = (9408.57842795923, 4299.8353044173)
x = x / 10000 * 4000
y = 4000 - y / 10000 * 4000
image1.crop((int(x-100), int(y-100), int(x+100), int(y+100)))
Out[13]:
In [14]:
localmap_raster = os.path.join(ordnancedir, "OS OpenMap Local (Full Colour Raster) SE", "data")
image2 = PIL.Image.open(os.path.join(localmap_raster, "SE23SE.tif"))
In [15]:
x, y = (4408.57842795923, 4299.8353044173)
y = 5000 - y
image2.crop((int(x-100), int(y-100), int(x+100), int(y+100)))
Out[15]:
OSM has the advantage of covering the whole planet, and being free. A downside is that it's distinctly harder to access and process.
However, OSMNX is a bit opaque, and downloads data. We prefer to work with off-line data:
We then have 2 options:
We again use Illinois as an example: http://download.geofabrik.de/north-america/us/illinois.html
In [16]:
osm_path = os.path.join("/media", "disk", "OSM_Data")
file = os.path.join(osm_path, "illinois", "gis.osm_roads_free_1.shp")
frame = gpd.GeoDataFrame.from_file(file)
frame.head()
Out[16]:
In [17]:
# Perform some cleaning
frame = frame[frame.geometry.map(lambda x : x is not None)]
bds = frame.bounds
bad = bds.minx < -100
frame = frame[~bad]
In [18]:
# This is, erm, messy...
list(set(frame.fclass))[:10]
Out[18]:
In [19]:
fig, ax = plt.subplots(figsize=(12,12))
lc = matplotlib.collections.LineCollection(open_cp.plot.lines_from_geometry(frame.geometry),
color="black", linewidth=0.5)
ax.add_collection(lc)
patches = open_cp.plot.patches_from_geometry(areas.geometry,
fc="blue", ec="none", alpha=0.5, zorder=5)
ax.add_collection(matplotlib.collections.PatchCollection(patches))
ax.set_aspect(1.0)
xmin, ymin, xmax, ymax = areas.total_bounds
xd, yd = xmax - xmin, ymax - ymin
ax.set(xlim=[xmin - xd/10, xmax + xd/10], ylim=[ymin - yd/10, ymax + yd/10])
None
In [21]:
fig, ax = plt.subplots(ncols=2, figsize=(16,10))
ax[0].set_aspect(1)
lc = matplotlib.collections.LineCollection(open_cp.plot.lines_from_geometry(frame.geometry),
color="black", linewidth=0.5)
ax[0].add_collection(lc)
ax[0].set(title="OpenStreetMap")
ax[1].set_aspect(1)
lc = matplotlib.collections.LineCollection(open_cp.plot.lines_from_geometry(tiger_frame.geometry),
color="black", linewidth=0.5)
ax[1].add_collection(lc)
ax[1].set(title="TIGER")
for a in ax:
a.set(xlim=[-87.7, -87.6], ylim=[41.7, 41.8])
In [ ]: