Easy OpenStreetMap tile displaying for Python

Question

I want to include the open street map (OSM) in my python code.

I have read through lots of webpages regarding to OSM. But unfortunately I'm a bit lost, regarding which package I use best.

I'm looking for an easy way to get an OSM image in my app. As I starting point I'm thinking of something like:

import matplotlib.pyplot as plt

# Pseudo - Code for required function 'GetOSMImage'
Map = GetOSMImage(lat,long,delta_lat,delta_long)

imgplot = plt.imshow(Map)

Later I want to add plot my additional data in this plt. (I'm aware that I'll need to deal with projections etc.)

What I don't need/want:

• Display on my own website
• To upload my data to some Internet Server
• interactive features like zooming, scrolling (in the first place)
• manually process and render the .xml data from OSM
• In the first place I don't want to define every detail of the rendering style. I hope/expect that there exists some default styles.

Do you have a good starting point for me? Or do I underestimate the complexity of this topic?

Answer 1

Based on your input, I was able to achive my target. Here is my code for others, which are searching a starting point to OSM. (Of course there is still much room for improvements).

Update

Please respect the usage policy of Open Street Map!

OpenStreetMap data is free for everyone to use. Our tile servers are not.

Requirements

• Heavy use (e.g. distributing an app that uses tiles from openstreetmap.org) is forbidden without prior permission from the Operations Working Group. See below for alternatives.
• Clearly display license attribution.
• Do not actively or passively encourage copyright infringement.
• Calls to /cgi-bin/export may only be triggered by direct end-user action. (For example: “click here to export”.) The export call is an expensive (CPU+RAM) function to run and will frequently reject when server is under high load.
• Recommended: Do not hardcode any URL at tile.openstreetmap.org as doing so will limit your ability to react quickly if the service is disrupted or blocked.
• Recommended: add a link to https://www.openstreetmap.org/fixthemap to allow your users to report and fix problems in our data.

Technical Usage Requirements

• Valid HTTP User-Agent identifying application. Faking another app’s User-Agent WILL get you blocked.
• If known, a valid HTTP Referer.
• DO NOT send no-cache headers. (“Cache-Control: no-cache”, “Pragma: no-cache” etc.)
• Cache Tile downloads locally according to HTTP Expiry Header, alternatively a minimum of 7 days.
• Maximum of 2 download threads. (Unmodified web browsers’ download thread limits are acceptable.)

More details see: https://operations.osmfoundation.org/policies/tiles/

Here is the code:

import matplotlib.pyplot as plt
import numpy as np

import math
import urllib2
import StringIO
from PIL import Image



def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)
  
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)
  
  
    
def getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom):
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)
    
    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl=smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = urllib2.urlopen(imgurl).read()
                tile = Image.open(StringIO.StringIO(imgstr))
                Cluster.paste(tile, box=((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None

    return Cluster
    
   
  
if __name__ == '__main__':
    
    a = getImageCluster(38.5, -77.04, 0.02,  0.05, 13)
    fig = plt.figure()
    fig.patch.set_facecolor('white')
    plt.imshow(np.asarray(a))
    plt.show()
    

Answer 2

Building up on BerndGit's nice answer, I add a slightly modified version which allows to display other contents together with the tiles (using Basemap). Btw I've come across a dedicated library, geotiler (http://wrobell.it-zone.org/geotiler/intro.html), but it requires Python 3.

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np

import math
import urllib2
import StringIO
from PIL import Image

def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)

def num2deg(xtile, ytile, zoom):
  """
  http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
  This returns the NW-corner of the square. 
  Use the function with xtile+1 and/or ytile+1 to get the other corners. 
  With xtile+0.5 & ytile+0.5 it will return the center of the tile.
  """
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)

def getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom):
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax = deg2num(lat_deg, lon_deg, zoom)
    xmax, ymin = deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)

    bbox_ul = num2deg(xmin, ymin, zoom)
    bbox_ll = num2deg(xmin, ymax + 1, zoom)
    #print bbox_ul, bbox_ll

    bbox_ur = num2deg(xmax + 1, ymin, zoom)
    bbox_lr = num2deg(xmax + 1, ymax +1, zoom)
    #print bbox_ur, bbox_lr

    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) )
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl=smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = urllib2.urlopen(imgurl).read()
                tile = Image.open(StringIO.StringIO(imgstr))
                Cluster.paste(tile, box=((xtile-xmin)*255 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None

    return Cluster, [bbox_ll[1], bbox_ll[0], bbox_ur[1], bbox_ur[0]]

if __name__ == '__main__':
    lat_deg, lon_deg, delta_lat,  delta_long, zoom = 45.720-0.04/2, 4.210-0.08/2, 0.04,  0.08, 14
    a, bbox = getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom)

    fig = plt.figure(figsize=(10, 10))
    ax = plt.subplot(111)
    m = Basemap(
        llcrnrlon=bbox[0], llcrnrlat=bbox[1],
        urcrnrlon=bbox[2], urcrnrlat=bbox[3],
        projection='merc', ax=ax
    )
    # list of points to display (long, lat)
    ls_points = [m(x,y) for x,y in [(4.228, 45.722), (4.219, 45.742), (4.221, 45.737)]]
    m.imshow(a, interpolation='lanczos', origin='upper')
    ax.scatter([point[0] for point in ls_points],
               [point[1] for point in ls_points],
               alpha = 0.9)
    plt.show()

Answer 3

Edit: OpenStreetMap states that their tile servers are not free to use and are under a usage policy:
https://operations.osmfoundation.org/policies/tiles/
Please read this before using the example.

As I had problems implementing the code in Python 3.8, I merged a few of the answers together and modified the code. Now it works for me and I don't get any errors.
When I tried to run the original code from BerndGit in Python 3, I had to make the same changes as Joining Dots described in his answer. I replaced

 import urllib2
 import StringIO

with

import requests
from io import BytesIO

because the urllib2 library doesn't work with Python 3 anymore. You have to use urllib.request or requests.
Then I had to change these two lines from the getImageCluster function

imgstr = urllib2.urlopen(imgurl).read()
tile = Image.open(StringIO.StringIO(imgstr))

to

imgstr = requests.get(imgurl)
tile = Image.open(BytesIO(imgstr.content))

After that I could run the code without errors but it still couldn't download the images. I always got a black tile as a result. Through some research I learned, that it is important to fake a user agent while using requests, as the website could tell that the request is coming from Python and may block it. The following website describes this:
https://www.scrapehero.com/how-to-fake-and-rotate-user-agents-using-python-3/
So I followed the suggestions from the website which resulted in adding this line right at the beginning of the getImageCluster function:

headers = {"User-Agent":"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36"}

Now we need to include these headers into the requests call:

imgstr = requests.get(imgurl, headers=headers)

The whole code looks like this now:

import matplotlib.pyplot as plt
import numpy as np
import math
import requests
from io import BytesIO
from PIL import Image
   
    
    
def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)
    
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)
   
def getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom):
    headers = {"User-Agent":"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36"}
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)
   
    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl = smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = requests.get(imgurl, headers=headers)
                tile = Image.open(BytesIO(imgstr.content))
                Cluster.paste(tile, box = ((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None
   
    return Cluster
    
    
    
if __name__ == '__main__':
    
    a = getImageCluster(38.5, -77.04, 0.02,  0.05, 13)
    fig = plt.figure()
    fig.patch.set_facecolor('white')
    plt.imshow(np.asarray(a))
    plt.show()

The result is the following:

Answer 4

It is not so very complex. A little bit of guidance can be obtained from this link, where the complexity of tiles are explained in detail.

It can hardly be reproduced here, but in general you have to

• determine the tiles you need by formula
• load them from their server (there is a certain choice of map styles)
• possibly concatenate them in both directions
• and then display them.

Be aware that you possibly have aspect ratio issues which you must solve as well...

Answer 5

Using python 3.6.5 you need to modify the header a bit:

import matplotlib.pyplot as plt
import numpy as np
import math
import urllib3
from io import StringIO
from PIL import Image

simply use pip install and be aware that PIL has to be implemented like pip install Pillow

Answer 6

Yet another way to get combined openstreetmap image (with python3, amazing mercantile library and parallel fetching):

import multiprocessing
import random
import io
import mercantile
import urllib.request
import PIL.Image

def _download_tile(tile: mercantile.Tile):
    """
    Helper function for downloading associated image
    """
    server = random.choice(['a', 'b', 'c'])
    url = 'http://{server}.tile.openstreetmap.org/{zoom}/{x}/{y}.png'.format(
        server=server,
        zoom=tile.z,
        x=tile.x,
        y=tile.y
    )
    response = urllib.request.urlopen(url)
    img = PIL.Image.open(io.BytesIO(response.read()))

    return img, tile    

def get_image(west, south, east, north, zoom):
    """
    return glued tiles as PIL image
    :param west: west longitude in degrees
    :param south: south latitude in degrees
    :param east: east longitude in degrees
    :param north: north latitude in degrees
    :param zoom: wanted size
    :return: Image
    """
    tiles = list(mercantile.tiles(west, south, east, north, zoom))

    tile_size = 256
    min_x = min_y = max_x = max_y = None

    for tile in tiles:
        min_x = min(min_x, tile.x) if min_x is not None else tile.x
        min_y = min(min_y, tile.y) if min_y is not None else tile.y
        max_x = max(max_x, tile.x) if max_x is not None else tile.x
        max_y = max(max_y, tile.y) if max_y is not None else tile.y

    out_img = PIL.Image.new(
        'RGB',
        ((max_x - min_x + 1) * tile_size, (max_y - min_y + 1) * tile_size)
    )

    pool = multiprocessing.Pool(8)
    results = pool.map(_download_tile, tiles)
    pool.close()
    pool.join()

    for img, tile in results:
        left = tile.x - min_x
        top = tile.y - min_y
        bounds = (left * tile_size, top * tile_size, (left + 1) * tile_size, (top + 1) * tile_size)
        out_img.paste(img, bounds)

    return out_img   

if __name__ == '__main__':
    # get combined image and save to file
    get_image(west=103, south=51, east=110, north=56, zoom=8).save('osm_image.png')

Answer 7

The following is also based on BerndGit's wonderful answer. I had to do some modifications to get it working with Python 3.6.7. Posting them here in case it helps others.

Set-up required Pillow, and replacing urllib with requests, and replacing io/StringIO with io/ByesIO

import requests
from io import BytesIO

And then just needed to modify how the image is downloaded in the getImageCluster() function:

imgstr = requests.get(imgurl)
tile = Image.open(BytesIO(imgstr.content))

Big thanks to BerndGit for going to the trouble of posting the original.

Haven't managed to get Etna's modified Basemap version working yet. Had to add in an export path for the PROJ_LIB error for Basemap:

export PROJ_LIB=/path/to/your/instalation/of/anaconda/share/proj/

(solution at Basemap import error in PyCharm —— KeyError: 'PROJ_LIB')

And getting a set attribute error when trying to plot. It occurs using the Basemap tutorial too (https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#plot) but with the difference that the scatter of data does still plot as a layer on top of the map despite the error. With the OSM tiles, cannot get the data layer to show on top of the map. Having to export each layer individually and then combine using image editing software.

Answer 8

Sorry to wake an old question but came across it when i tried to do the exact same thing.

You can use the contextily package which does just that.
From the docs:

contextily is a small Python ... package to retrieve tile maps from the internet. It can add those tiles as basemap to matplotlib figures or write tile maps to disk into geospatial raster files.

Showing just a slice of a map

Sample code:

import contextily as cx
import matplotlib.pyplot as plt

west, south, east, north = (long, lat - delta_lat, long + delta_long, lat)
im, bbox = cx.bounds2img(west, south, east, north, ll=True, source=cx.providers.OpenStreetMap.Mapnik)
plt.imshow(im)
plt.show()

And then a window will open with a shown figure of the map slice your requested:

Showing a map with embedded data

You mentioned that later you'd like to add data embedded on the map, so you'll need a slightly different approach. For showing a polygon you'd need to use geopandas with a geojson, then let it plot itself, and then add a basemap with contextily:

import contextily as cx
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import shape

geojson = {'type': 'Polygon', 'coordinates': [[[34, 31], [32, 31], [34, 33]]]}
gdf = gpd.GeoDataFrame({'geometry':[shape(geojson)]}, crs='EPSG:4326')
ax = gdf.plot(facecolor='none', edgecolor='red', linewidth=2)
cx.add_basemap(ax=ax, crs='EPSG:4326', source=cx.providers.OpenStreetMap.Mapnik)
plt.show()

Resulting with a window showing a polygon on map: