Easy display of the OpenStreetMap panel for Python

Question

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

I read a lot of web pages related to OSM. But unfortunately, I am a little lost as to which package I use best.

I am looking for an easy way to get an OSM image in my application. As a starting point, I think 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 a plot to my additional data in this plt. (I know that I need to deal with projections, etc.)

What I don't need / need:

• Show on my own website
• To upload my data to any internet server
• interactive features like scaling, scrolling (first)
• manually process and display .xml data from OSM
• First, I don’t want to define every detail of the rendering style. I hope / expect some default styles exist.

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

Answer 1

Based on your input, I was able to achieve my goal. Here is my code for others looking for a starting point for OSM. (Of course, there is still much room for improvement).

 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

By creating a good BerndGit answer, I am adding a slightly modified version that allows other content to be displayed along with tiles (using Basemap). Btw I came across a dedicated library, geotiler ( http://wrobell.it-zone.org/geotiler/intro.html ), but this 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

It is not that difficult. A bit of guidance can be obtained from this link, which explains in detail the complexity of the tiles.

It's hard to reproduce here, but overall you have to

• determine the tiles you need
• download them from your server (there is a certain choice of map styles)
• it is possible to combine them in both directions
• and then display them.

Remember that you may have aspect ratios that you must also solve ...

Answer 4

Using python 3.6.5 you need to change 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 

just use pip install and keep in mind that PIL should be implemented as pip install Pillow

Answer 5

Another way to get a combined openstreetmap image (using python3, awesome product 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 6

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

Requires Pillow configuration, replacing urllib with queries, and replacing io / StringIO with io / ByesIO

 import requests from io import BytesIO 

And then you just need to change the way the image is loaded into the getImageCluster () function:

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

Many thanks to BerndGit for publishing the original.

Failed to get Etna a modified version of Basemap to work. I had to add the PROJ_LIB error for the base card to the export path:

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

(the solution is in the import error of the base map in PyCharm —— KeyError: & # 39; PROJ_LIB & # 39; )

And getting the error of the set attribute when trying to build. This also happens with the help of the basemap tutorial ( https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#plot ), but with the difference that the data scatter still appears as a layer on top of the map, despite the error . With OSM tiles, I cannot display a data layer on top of a map. The need to export each layer separately, and then merge using image editing software.