Łatwe wyświetlanie kafelków OpenStreetMap dla Python

Question

Chcę dołączyć otwartą mapę ulic (OSM) w moim kodzie Pythona.

Przeczytałem wiele stron internetowych dotyczących OSM. Ale niestety jestem trochę zagubiony, w odniesieniu do którego pakietu używam najlepiej.

Szukam łatwego sposobu na uzyskanie obrazu OSM w mojej aplikacji. Jak począwszy punkt mam na myśli coś takiego:

import matplotlib.pyplot as plt 

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

imgplot = plt.imshow(Map) 

Później chcę dodać działce moich danych dodatkowych w tym PLT. (Jestem świadomy, że będę musiał radzić sobie z występami itp)

Co nie muszę/chcą:

• Pokaż na własnej stronie
• Aby przesłać swoje dane na niektóre Internet Server
• interaktywne funkcje, takie jak powiększanie, przewijanie (w pierwszej kolejności)
• ręcznie przetwarzać i renderować dane .xml z OSM
• w pierwszej kolejności nie chcę definiować każdy szczegół stylu renderowania . Mam nadzieję/oczekuję, że istnieją pewne domyślne style.

Czy masz dobry punkt wyjścia dla mnie? Czy mogę nie docenić złożoności tego tematu?

Answer 1

Na podstawie danych wejściowych, udało mi się achive mój cel. Oto mój kod dla innych, którzy szukają punktu wyjścia do OSM. (Corse nadal ma wiele do zrobienia).

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

To nie jest takie skomplikowane. Trochę wskazówek można uzyskać z linku this, gdzie szczegółowo objaśniono złożoność płytek.

Może trudno być powielana tutaj, ale generalnie trzeba

• określić płytek trzeba przez formula
• załadować ich z serwera (istnieje pewien wybór stylów map)
• Ewentualnie połączyć je w obu kierunkach
• , a następnie wyświetlić je.

Należy pamiętać, że może mieć problemy z ratio aspektem, który trzeba rozwiązać, jak również ...

Answer 3

Opierając się na miłą odpowiedź BerndGit za dodam nieco zmodyfikowaną wersję, która pozwala na wyświetlanie innych treści wraz z płytkami (używając bazowa). Btw ja natrafiłem dedykowany biblioteki, geotiler (http://wrobell.it-zone.org/geotiler/intro.html), ale wymaga Pythona 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()