Easy OpenStreetMap tile displaying for Python

问题:

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?

回答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 corse there is still much room for improvements).

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()

回答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()

回答3:

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...

回答4:

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

回答5:

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')

回答6:

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.