Download OpenStreetMap bounding box PNG

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):
    lat_rad = math.radians(lat_deg)
    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)

def download_map(zoom, lat1, lon1, lat2, lon2):
    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)
                print("Downloading tile {}".format(url))
                req = Request(url)
                ua = "Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:90.0) Gecko/20100101 Firefox/90.0"
                req.add_header("User-Agent", ua)  # OSM seems to not like Python's default UA.
                resp = urlopen(req)
                body = resp.read()
                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.
download_map(12, 38, -122.7, 37.2, -121.7)