aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorBryan Brattlof <hello@bryanbrattlof.com>2020-11-12 18:26:05 -0500
committerBryan Brattlof <hello@bryanbrattlof.com>2020-11-12 18:46:37 -0500
commit57c5bc143911052e3a3e6428059f2e3ba7aa24f6 (patch)
treec7944a7c290e32497ccf2055608dc6e6e9c49b2d
parent41cb4761013bb76c0ccc56d48d4512e35f5251c6 (diff)
downloadnorta-57c5bc143911052e3a3e6428059f2e3ba7aa24f6.tar.gz
norta-57c5bc143911052e3a3e6428059f2e3ba7aa24f6.tar.bz2
add code to create article visuals
This included the example code and code to create all visuals for my 'Adding OpenStreetMaps to Matplotlib' article.
-rw-r--r--add-osm-to-mpl.py212
-rw-r--r--readme.rst24
2 files changed, 232 insertions, 4 deletions
diff --git a/add-osm-to-mpl.py b/add-osm-to-mpl.py
new file mode 100644
index 0000000..13b275b
--- /dev/null
+++ b/add-osm-to-mpl.py
@@ -0,0 +1,212 @@
+import pandas as pd
+import os.path
+
+DATA_FILE = 'data/bus.csv'
+assert os.path.isfile(DATA_FILE) is True,\
+ "dataset not found! did you run 'prepare-data.py?"
+
+# load the dataset
+df = pd.read_csv(
+ 'data/bus.csv',
+ dtype={
+ 'epoch': 'str',
+ 'vid': 'category',
+ 'lat': 'float32',
+ 'lon': 'float32',
+ 'hdg': 'Int16',
+ 'des': 'category',
+ 'dly': 'boolean',
+ 'pdist': 'float32'
+ },
+ parse_dates=[
+ 'epoch'
+ ],
+)
+df.set_index('epoch')
+
+# clean the data
+_len = len(df.index)
+df = df[(df['lat'] < 30.386759) & ( df['lat'] > 29.587366) &
+ (df['lon'] > -90.874932) & (df['lon'] < -89.513523)]
+print("Removed {} Rows".format(_len - len(df.index)))
+
+# narrow down to just route 16
+df = df[(df['des'].str.contains('S. Claiborne Canal St via Poydras St.', na=False))|
+ (df['des'].str.contains('S. Claiborne Poydras St. to S. Carrollton Ave.'))]
+
+#################################
+# Main Visual
+#################################
+from matplotlib import pyplot as plt
+import basemap
+
+avg_lat, avg_lon = df.lat.mean(), df.lon.mean()
+std_lat, std_lon = df.lat.std(), df.lon.std()
+
+top = avg_lat + (9 * std_lat)
+rgt = avg_lon + (4 * std_lon)
+bot = avg_lat - (9 * std_lat)
+lef = avg_lon - (4 * std_lon)
+
+img = basemap.image(top, rgt, bot, lef, zoom=13,
+ url="http://c.tile.stamen.com/toner/{z}/{x}/{y}.png")
+
+fig, ax = plt.subplots(figsize=(12,8))
+ax.scatter(df.lon, df.lat, alpha=0.1, c='red', s=1)
+ax.imshow(img, extent=(lef, rgt, bot, top), aspect= 'equal')
+plt.savefig('visuals/route-16-plot.png')
+
+#################################
+# Self-Portrait Visual
+#################################
+from matplotlib import pyplot as plt
+
+img = plt.imread('visuals/self-portrait.png')
+plt.imshow(img)
+
+plt.savefig('visuals/self-portrait-plot.png')
+
+#################################
+# Single Tile Visual
+#################################
+URL = "https://tile.openstreetmap.org/{z}/{x}/{y}.png".format
+
+import math
+TILE_SIZE = 256
+
+def point_to_pixels(lon, lat, zoom):
+ """convert gps coordinates to web mercator"""
+ r = math.pow(2, zoom) * TILE_SIZE
+ lat = math.radians(lat)
+
+ x = int((lon + 180.0) / 360.0 * r)
+ y = int((1.0 - math.log(
+ math.tan(lat) + (1.0 / math.cos(lat))) / math.pi) / 2.0 * r)
+
+ return x, y
+
+zoom = 16
+x, y = point_to_pixels(-90.064279, 29.95863, zoom)
+
+x_tile, y_tile = int(x / TILE_SIZE), int(y / TILE_SIZE)
+
+from io import BytesIO
+from PIL import Image
+import requests
+
+# format the url
+url = URL(x=x_tile, y=y_tile, z=zoom)
+
+# make the request
+with requests.get(url) as resp:
+ img = Image.open(BytesIO(resp.content))
+
+# plot the tile
+plt.imshow(img)
+plt.savefig('visuals/french-quarter-plot.png')
+
+#################################
+# Multi Tile Visual
+#################################
+top, bot = df.lat.max(), df.lat.min()
+lef, rgt = df.lon.min(), df.lon.max()
+
+zoom = 13
+x0, y0 = point_to_pixels(lef, top, zoom)
+x1, y1 = point_to_pixels(rgt, bot, zoom)
+
+x0_tile, y0_tile = int(x0 / TILE_SIZE), int(y0 / TILE_SIZE)
+x1_tile, y1_tile = math.ceil(x1 / TILE_SIZE), math.ceil(y1 / TILE_SIZE)
+
+assert (x1_tile - x0_tile) * (y1_tile - y0_tile) < 50, "That's too many tiles!"
+
+from itertools import product
+
+# full size image we'll add tiles to
+img = Image.new('RGB', (
+ (x1_tile - x0_tile) * TILE_SIZE,
+ (y1_tile - y0_tile) * TILE_SIZE))
+
+# loop through every tile inside our bounded box
+for x_tile, y_tile in product(range(x0_tile, x1_tile),
+ range(y0_tile, y1_tile)):
+
+ with requests.get(URL(x=x_tile, y=y_tile, z=zoom)) as resp:
+ tile_img = Image.open(BytesIO(resp.content))
+
+ # add each tile to the full size image
+ img.paste(
+ im=tile_img,
+ box=((x_tile - x0_tile) * TILE_SIZE,\
+ (y_tile - y0_tile) * TILE_SIZE))
+
+plt.figure(figsize=(12,8))
+plt.imshow(img)
+plt.savefig('visuals/the-basemap.png', bbox_inches ="tight")
+
+
+#################################
+# Crop Lines Visual
+# Don't include this in article
+#################################
+from PIL import ImageDraw
+img2 = img.copy()
+
+# find the mercator coordinates of the top-left corner of all
+# the tiles we downloaded from OpenStreetMap
+x, y = x0_tile * TILE_SIZE, y0_tile * TILE_SIZE
+
+# draw a red rectangle of the part of the image we want to keep
+draw = ImageDraw.Draw(img2)
+draw.rectangle(
+ (x0 - x, # left
+ y0 - y, # top
+ x1 - x, # right
+ y1 - y), # bottom
+ outline=(255, 0, 0),
+ width=5)
+
+# draw black rectangles around each individual tile
+for x_tile, y_tile in product(range(0, (x1_tile - x0_tile)),
+ range(0, (y1_tile - y0_tile))):
+
+ x, y = x_tile * TILE_SIZE, y_tile * TILE_SIZE
+ draw.rectangle(
+ (x - 2, # left
+ y - 2, # top
+ (x + TILE_SIZE), # right
+ (y + TILE_SIZE)), # bottom
+ outline=(0, 0, 0),
+ width=2)
+
+plt.figure(figsize=(12,8))
+plt.imshow(img2)
+plt.savefig('visuals/basemap-cropping-lines.png', bbox_inches ="tight")
+
+#################################
+# Cropping The Multi Tile Visual
+#################################
+x, y = x0_tile * TILE_SIZE, y0_tile * TILE_SIZE
+
+img = img.crop((
+ x0 - x, # left
+ y0 - y, # top
+ x1 - x, # right
+ y1 - y)) # bottom
+
+plt.figure(figsize=(12,8))
+plt.imshow(img)
+plt.savefig('visuals/basemap-cropped.png', bbox_inches ="tight")
+
+#################################
+# Final Visual
+#################################
+fig, ax = plt.subplots(figsize=(12,12))
+
+ax.set_ylim(bot, top)
+ax.set_xlim(lef, rgt)
+
+ax.scatter(df.lon, df.lat, alpha=0.1, c='red', s=1)
+ax.imshow(img, extent=(lef, rgt, bot, top), aspect="equal")
+
+plt.savefig('visuals/final-visual.png')
diff --git a/readme.rst b/readme.rst
index d4fb185..38706c8 100644
--- a/readme.rst
+++ b/readme.rst
@@ -10,6 +10,13 @@ you are running at least Python 3.6. All other dependencies are listed in
$ pip install -r requirements.txt
+The current versions that work on *my computer*\ :sup:`tm` are:
+
+- matplotlib=3.3.2
+- pandas=1.1.3
+- pillow=7.0.0
+- requests=2.24.0
+
The Data-set
############
@@ -79,6 +86,10 @@ dependencies
It assumes you have `requests <https://requests.readthedocs.io/en/master/>`__, and
`pillow <https://python-pillow.org/>`__ libraries installed.
+.. code-block::
+
+ $ pip install requests pillow
+
usage
-----
@@ -95,7 +106,7 @@ the ``top``, ``rgt``, ``bot``, ``lef`` arguments along with the appropriate
img = basemap.image(top, rgt, bot, lef, zoom=13)
-The module will return a ``Pillow.Image`` object that you can add to your
+The module will return a ``Pillow.Image()`` object that you can add to your
matplotlib visuals like this:
.. code-block:: python
@@ -112,15 +123,20 @@ You can also use ``url`` to specify which tile servers you want to use:
img = basemap.image(top, rgt, bot, lef, zoom=13,
url="http://c.tile.stamen.com/toner/{z}/{x}/{y}.png")
-Any extra arguments to needed for the ``url`` argument can be passed along as key word arguments in the ``basemap.image()`` function. For example:
+Any extra arguments to format the ``url`` argument can be passed along as key word arguments in the ``basemap.image()`` function. For example:
.. code-block:: python
img = basemap.image(top, rgt, bot, lef, zoom=13, api=API_KEY
url="http://tileserver.example.com/{api}/{z}/{x}/{y}.png")
-The full write-up of how this module works is available at
-https://bryanbrattlof.com/adding-openstreetmaps-to-matplotlib/
+add-osm-to-mpl.py
+#################
+
+`add-osm-to-mpl.py <https://git.bryanbrattlof.com/norta/tree/add-osm-to-mpl.py>`__:
+holds all the example code and code to generate the visuals I used in my `Adding
+OpenStreetMaps To MatplotLib <https://bryanbrattlof.com/
+adding-openstreetmaps-to-matplotlib/>`__ article.
Contributing
############