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)