add module to generate OSM basemaps

this is a simple script to ease making visuals for other projects I'm
working on with this dataset.

2 files changed, 137 insertions, 6 deletions

diff --git a/basemap.py b/basemap.py
new file mode 100644
index 0000000..86e2ea7
--- /dev/null
+++ b/basemap.py
@@ -0,0 +1,70 @@
+
+from itertools import product
+from io import BytesIO
+from PIL import Image
+import requests
+import math
+
+URL = "https://tile.openstreetmap.org/{z}/{x}/{y}.png"
+RADIUS = 6371   # earth's average radius in km
+TILE_SIZE = 256 # each tile's size in pixels
+
+
+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
+
+
+def tile(session=None, **url_format_args):
+    """download tile as PIL.Image from Tile Server API"""
+    url = url_format_args.get('url', URL)
+    if not session:
+        session = requests
+
+    with session.get(url.format(**url_format_args)) as resp:
+        return Image.open(BytesIO(resp.content))
+
+
+def image(top, right, bottom, left, zoom=14, **tile_args):
+    """return an osm map from the given bounding box points"""
+
+    # convert gps coordinates to web mercator
+    x0, y0 = point_to_pixels(left, top, zoom)
+    x1, y1 = point_to_pixels(right, bottom, zoom)
+
+    # conver pixel coordinates to tiles
+    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!"
+
+    # crate a single large image to add the tiles to
+    img = Image.new('RGB', (
+        int(x1_tile - x0_tile) * TILE_SIZE,
+        int(y1_tile - y0_tile) * TILE_SIZE))
+
+    # download all tiles
+    with requests.Session() as session:
+        for x_tile, y_tile in product(
+                range(x0_tile, x1_tile), range(y0_tile, y1_tile)):
+
+            img.paste(
+                im=tile(session, x=x_tile, y=y_tile, z=zoom, **tile_args),
+                box=(int(x_tile - x0_tile) * TILE_SIZE,
+                     int(y_tile - y0_tile) * TILE_SIZE))
+
+    # crop image to the given bounding box
+    x, y = x0_tile * TILE_SIZE, y0_tile * TILE_SIZE
+    return img.crop((
+        int(x0 - x),  # left
+        int(y0 - y),  # top
+        int(x1 - x),  # right
+        int(y1 - y))) # bottom
diff --git a/readme.rst b/readme.rst
index 2d47453..d4fb185 100644
--- a/readme.rst
+++ b/readme.rst
@@ -25,8 +25,8 @@ the API. This is on the larger end of what I feel comfortable publishing online.
 So if you wish for a copy please `send a DM or email me
 <https://bryanbrattlof.com/connect/>`__ and I'll gladly send it to you.
 
-Preparing The Data
-##################
+prepare-data.py
+###############
 
 `prepare-data.py <https://git.bryanbrattlof.com/norta/tree/prepare-data.py>`__:
 is a small script to convert the ``bus.log.tar.gz`` file into a CSV file named
@@ -53,12 +53,73 @@ is a small script to convert the ``bus.log.tar.gz`` file into a CSV file named
     )
     df.set_index('epoch')
 
-Base-map
-########
+usage
+-----
 
-**TODO**
+With the ``bus.log.tar.gz`` inside the ``data`` directory, simply run the python
+script to generate the ``bus.csv`` file inside the ``data`` directory like so:
 
-The full write-up is available at
+.. code-block::
+
+   $ python prepare-data.py
+
+There is no data cleaning involved. This script will only decompress and convert
+the log file into a csv file.
+
+basemap.py
+##########
+
+`basemap.py <https://git.bryanbrattlof.com/norta/tree/basemap.py>`__: is a simple
+module that downloads and combines OpenStreetMap tiles into one large image you
+can add to your matplotlib visuals.
+
+dependencies
+------------
+
+It assumes you have `requests <https://requests.readthedocs.io/en/master/>`__, and
+`pillow <https://python-pillow.org/>`__ libraries installed.
+
+usage
+-----
+
+To use this module, pass the bounding box in GPS coordinates of the area into
+the ``top``, ``rgt``, ``bot``, ``lef`` arguments along with the appropriate
+``zoom`` level:
+
+.. code-block:: python
+
+   import basemap
+
+   top, bot = df.lat.max(), df.lat.min()
+   lef, rgt = df.lon.min(), df.lon.max()
+
+   img = basemap.image(top, rgt, bot, lef, zoom=13)
+
+The module will return a ``Pillow.Image`` object that you can add to your
+matplotlib visuals like this:
+
+.. code-block:: python
+
+   from matplotlib import pyplot as plt
+   fig, ax = plt.subplots()
+   ax.imshow(img, extent=(lef, rgt, bot, top), aspect= 'equal')
+   plt.show()
+
+You can also use ``url`` to specify which tile servers you want to use:
+
+.. code-block:: python
+
+   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:
+
+.. 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/
 
 Contributing