## Tuesday, August 03, 2021

For a project where I was plotting some locations using Matplotlib, I needed a way to get a PNG map from a bounding box lat-long pair. Here's a Python function I came up with.

```from os import listdir, mkdir
from os.path import exists
from PIL import Image
from random import uniform
from time import sleep
from urllib.request import Request, urlopen
import math

# Similar to https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python .
def deg2float(lat_deg, lon_deg, zoom):
n = 2.0 ** zoom
xtile = (lon_deg + 180.0) / 360.0 * n
ytile = (1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n
return (xtile, ytile)

lon_start, lon_end = min(lon1, lon2), max(lon1, lon2)
lat_start, lat_end = max(lat1, lat2), min(lat1, lat2)

# Top left corner of bounding box.
x1, y1 = deg2float(lat_start, lon_start, zoom)
x1i, y1i = math.floor(x1), math.floor(y1)

# Bottom right corner of bounding box.
x2, y2 = deg2float(lat_end, lon_end, zoom)
x2i, y2i = math.ceil(x2), math.ceil(y2)

x_cnt, y_cnt = abs(x1i - x2i), abs(y1i - y2i)
if x_cnt*y_cnt > 250:
err = "Too many tiles. Probably too big an area at too high a zoom level."
err += " See https://operations.osmfoundation.org/policies/tiles/ ."
raise Exception(err)

if not exists("maptiles"):
mkdir("maptiles")

for x in range(x_cnt):
for y in range(y_cnt):
xt, yt = x + x1i, y + y1i
path = "maptiles/{}_{}_{}.png".format(zoom, xt, yt)

if not exists(path):
sleep(uniform(0.5, 1.5))
url = "https://tile.openstreetmap.org/{}/{}/{}.png".format(zoom, xt, yt)
req = Request(url)
ua = "Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:90.0) Gecko/20100101 Firefox/90.0"
resp = urlopen(req)
with open(path, "wb") as f:
f.write(body)

im = Image.open("maptiles/{}_{}_{}.png".format(zoom, x1i, y1i))
tile_w, tile_h = im.size
total_w = x_cnt*tile_w
total_h = y_cnt*tile_h

new_im = Image.new("RGB", (total_w, total_h))

for x in range(x_cnt):
for y in range(y_cnt):
xt, yt = x + x1i, y + y1i
im = Image.open("maptiles/{}_{}_{}.png".format(zoom, xt, yt))
new_im.paste(im, (x*tile_w, y*tile_h))

cropped_w = round((x2 - x1)*tile_w)
cropped_h = round((y2 - y1)*tile_h)
cropped_im = Image.new("RGB", (cropped_w, cropped_h))
translate_x = round(-(x1 - x1i)*tile_w)
translate_y = round(-(y1 - y1i)*tile_h)
cropped_im.paste(new_im, (translate_x, translate_y))
cropped_im.save("map.png")

# Download a map of the SF Bay Area at zoom level 12. Approx 3000*3000px.