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.

view plainprint?

1. from os import listdir, mkdir  
2. from os.path import exists  
3. from PIL import Image  
4. from random import uniform  
5. from time import sleep  
6. from urllib.request import Request, urlopen  
7. import math  
8.   
9. # Similar to https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python .  
10. def deg2float(lat_deg, lon_deg, zoom):  
11.     lat_rad = math.radians(lat_deg)  
12.     n = 2.0 ** zoom  
13.     xtile = (lon_deg + 180.0) / 360.0 * n  
14.     ytile = (1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n  
15.     return (xtile, ytile)  
16.   
17. def download_map(zoom, lat1, lon1, lat2, lon2):  
18.     lon_start, lon_end = min(lon1, lon2), max(lon1, lon2)  
19.     lat_start, lat_end = max(lat1, lat2), min(lat1, lat2)  
20.   
21.     # Top left corner of bounding box.  
22.     x1, y1 = deg2float(lat_start, lon_start, zoom)  
23.     x1i, y1i = math.floor(x1), math.floor(y1)  
24.   
25.     # Bottom right corner of bounding box.  
26.     x2, y2 = deg2float(lat_end, lon_end, zoom)  
27.     x2i, y2i = math.ceil(x2), math.ceil(y2)  
28.   
29.     x_cnt, y_cnt = abs(x1i - x2i), abs(y1i - y2i)  
30.     if x_cnt*y_cnt > 250:  
31.         err = "Too many tiles. Probably too big an area at too high a zoom level."  
32.         err += " See https://operations.osmfoundation.org/policies/tiles/ ."  
33.         raise Exception(err)  
34.   
35.     if not exists("maptiles"):  
36.         mkdir("maptiles")  
37.   
38.     for x in range(x_cnt):  
39.         for y in range(y_cnt):  
40.             xt, yt = x + x1i, y + y1i  
41.             path = "maptiles/{}_{}_{}.png".format(zoom, xt, yt)  
42.   
43.             if not exists(path):  
44.                 sleep(uniform(0.5, 1.5))  
45.                 url = "https://tile.openstreetmap.org/{}/{}/{}.png".format(zoom, xt, yt)  
46.                 print("Downloading tile {}".format(url))  
47.                 req = Request(url)  
48.                 ua = "Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:90.0) Gecko/20100101 Firefox/90.0"  
49.                 req.add_header("User-Agent", ua)  # OSM seems to not like Python's default UA.  
50.                 resp = urlopen(req)  
51.                 body = resp.read()  
52.                 with open(path, "wb") as f:  
53.                     f.write(body)  
54.   
55.     im = Image.open("maptiles/{}_{}_{}.png".format(zoom, x1i, y1i))  
56.     tile_w, tile_h = im.size  
57.     total_w = x_cnt*tile_w  
58.     total_h = y_cnt*tile_h  
59.   
60.     new_im = Image.new("RGB", (total_w, total_h))  
61.   
62.     for x in range(x_cnt):  
63.         for y in range(y_cnt):  
64.             xt, yt = x + x1i, y + y1i  
65.             im = Image.open("maptiles/{}_{}_{}.png".format(zoom, xt, yt))  
66.             new_im.paste(im, (x*tile_w, y*tile_h))  
67.   
68.     cropped_w = round((x2 - x1)*tile_w)  
69.     cropped_h = round((y2 - y1)*tile_h)  
70.     cropped_im = Image.new("RGB", (cropped_w, cropped_h))  
71.     translate_x = round(-(x1 - x1i)*tile_w)  
72.     translate_y = round(-(y1 - y1i)*tile_h)  
73.     cropped_im.paste(new_im, (translate_x, translate_y))  
74.     cropped_im.save("map.png")  
75.   
76. # Download a map of the SF Bay Area at zoom level 12. Approx 3000*3000px.  
77. download_map(12, 38, -122.7, 37.2, -121.7)  

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)