Python 的简易 OpenStreetMap 磁贴显示
Easy OpenStreetMap tile displaying for Python
我想在我的 python 代码中包含开放街道地图 (OSM)。
我已经阅读了很多关于 OSM 的网页。但不幸的是，关于我最好使用哪个包，我有点迷茫。
我正在寻找一种在我的应用程序中获取 OSM 图像的简单方法。作为我的起点，我在想这样的事情：
import matplotlib.pyplot as plt
# Pseudo - Code for required function 'GetOSMImage'
Map = GetOSMImage(lat,long,delta_lat,delta_long)
imgplot = plt.imshow(Map)
稍后我想在这个plt中添加我的额外数据。 （我知道我需要处理投影等）
我不need/want:
- 在我自己的网站上展示
- 将我的数据上传到某个 Internet 服务器
- 缩放、滚动等交互功能（首先）
- 手动处理和渲染来自 OSM
的 .xml 数据
- 首先我不想定义渲染风格的每一个细节。我 hope/expect 存在一些默认样式。
你有好的起点吗？
还是我低估了这个话题的复杂性？
没那么复杂。可以从 this link 中获得一点指导，其中详细解释了图块的复杂性。
这里很难复现，但总的来说还是得
- 通过formula
确定您需要的瓷砖
- 从他们的服务器加载它们（有一定的地图样式选择）
- 可能在两个方向上连接它们
- 然后显示它们。
请注意，您可能还存在必须解决的纵横比问题...
根据您的意见，我能够实现我的目标。这是我为其他人编写的代码，他们正在寻找 OSM 的起点。 （当然还有很大的改进空间）。
更新
请尊重开放街道地图的使用政策！
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.)
更多详情见：https://operations.osmfoundation.org/policies/tiles/
代码如下：
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()
在 BerndGit 的精彩回答的基础上，我添加了一个稍微修改过的版本，它允许与图块一起显示其他内容（使用底图）。顺便说一句，我遇到了一个专用的库，geotiler (http://wrobell.it-zone.org/geotiler/intro.html)，但它需要 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()
使用python 3.6.5 你需要稍微修改header:
import matplotlib.pyplot as plt
import numpy as np
import math
import urllib3
from io import StringIO
from PIL import Image
只需使用pip install
并注意 PIL 必须像 pip install Pillow
一样实现
另一种获取组合 openstreetmap 图像的方法（使用 python3，惊人的 mercantile 库和并行获取）：
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')
以下也是根据BerndGit的精彩回答整理而来。我必须进行一些修改才能使其与 Python 3.6.7 一起使用。将它们张贴在这里以防对其他人有所帮助。
设置需要的Pillow，将urllib替换为requests，将io/StringIO替换为io/ByesIO
import requests
from io import BytesIO
然后只需要修改getImageCluster()函数中图片的下载方式：
imgstr = requests.get(imgurl)
tile = Image.open(BytesIO(imgstr.content))
非常感谢 BerndGit 不厌其烦地发布原件。
尚未设法让 Etna 的修改底图版本正常工作。必须为底图的 PROJ_LIB 错误添加导出路径：
export PROJ_LIB=/path/to/your/instalation/of/anaconda/share/proj/
（Basemap import error in PyCharm —— KeyError: 'PROJ_LIB' 处的解决方案）
并在尝试绘制时出现设置属性错误。它也发生在使用底图教程 (https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#plot) 时，但不同之处在于，尽管存在错误，但数据的散点仍然绘制为地图顶部的图层。使用 OSM tiles，无法让数据层显示在地图顶部。必须单独导出每个图层，然后使用图像编辑软件进行组合。
编辑：OpenStreetMap 声明他们的图块服务器不能免费使用，并且遵循使用政策：
https://operations.osmfoundation.org/policies/tiles/
请在使用示例之前阅读此内容。
由于我在 Python 3.8 中实现代码时遇到问题，我将一些答案合并在一起并修改了代码。现在它对我有用，我没有收到任何错误。
当我在 Python 3 中尝试 运行 来自 BerndGit 的原始代码时，我不得不进行与他的回答中描述的连接点相同的更改。我替换了
import urllib2
import StringIO
和
import requests
from io import BytesIO
因为 urllib2 库不再适用于 Python 3。您必须使用 urllib.request 或请求。
然后我不得不从 getImageCluster 函数中更改这两行
imgstr = urllib2.urlopen(imgurl).read()
tile = Image.open(StringIO.StringIO(imgstr))
至
imgstr = requests.get(imgurl)
tile = Image.open(BytesIO(imgstr.content))
之后我可以 运行 代码没有错误，但它仍然无法下载图像。结果我总是得到一块黑色的瓷砖。通过一些研究我了解到，在使用请求时伪造用户代理很重要，因为网站可以判断请求来自 Python 并可能阻止它。以下网站对此进行了描述：
https://www.scrapehero.com/how-to-fake-and-rotate-user-agents-using-python-3/
所以我按照网站的建议在 getImageCluster 函数的开头添加了这一行：
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"}
现在我们需要将这些 headers 包含到请求调用中：
imgstr = requests.get(imgurl, headers=headers)
现在整个代码如下所示：
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()
我想在我的 python 代码中包含开放街道地图 (OSM)。
我已经阅读了很多关于 OSM 的网页。但不幸的是，关于我最好使用哪个包，我有点迷茫。
我正在寻找一种在我的应用程序中获取 OSM 图像的简单方法。作为我的起点，我在想这样的事情：
import matplotlib.pyplot as plt
# Pseudo - Code for required function 'GetOSMImage'
Map = GetOSMImage(lat,long,delta_lat,delta_long)
imgplot = plt.imshow(Map)
稍后我想在这个plt中添加我的额外数据。 （我知道我需要处理投影等）
我不need/want:
- 在我自己的网站上展示
- 将我的数据上传到某个 Internet 服务器
- 缩放、滚动等交互功能（首先）
- 手动处理和渲染来自 OSM 的 .xml 数据
- 首先我不想定义渲染风格的每一个细节。我 hope/expect 存在一些默认样式。
你有好的起点吗？ 还是我低估了这个话题的复杂性？
没那么复杂。可以从 this link 中获得一点指导，其中详细解释了图块的复杂性。
这里很难复现，但总的来说还是得
- 通过formula 确定您需要的瓷砖
- 从他们的服务器加载它们（有一定的地图样式选择）
- 可能在两个方向上连接它们
- 然后显示它们。
请注意，您可能还存在必须解决的纵横比问题...
根据您的意见，我能够实现我的目标。这是我为其他人编写的代码，他们正在寻找 OSM 的起点。 （当然还有很大的改进空间）。
更新
请尊重开放街道地图的使用政策！
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.)
更多详情见：https://operations.osmfoundation.org/policies/tiles/
代码如下：
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()
在 BerndGit 的精彩回答的基础上，我添加了一个稍微修改过的版本，它允许与图块一起显示其他内容（使用底图）。顺便说一句，我遇到了一个专用的库，geotiler (http://wrobell.it-zone.org/geotiler/intro.html)，但它需要 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()
使用python 3.6.5 你需要稍微修改header:
import matplotlib.pyplot as plt
import numpy as np
import math
import urllib3
from io import StringIO
from PIL import Image
只需使用pip install
并注意 PIL 必须像 pip install Pillow
另一种获取组合 openstreetmap 图像的方法（使用 python3，惊人的 mercantile 库和并行获取）：
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')
以下也是根据BerndGit的精彩回答整理而来。我必须进行一些修改才能使其与 Python 3.6.7 一起使用。将它们张贴在这里以防对其他人有所帮助。
设置需要的Pillow，将urllib替换为requests，将io/StringIO替换为io/ByesIO
import requests
from io import BytesIO
然后只需要修改getImageCluster()函数中图片的下载方式：
imgstr = requests.get(imgurl)
tile = Image.open(BytesIO(imgstr.content))
非常感谢 BerndGit 不厌其烦地发布原件。
尚未设法让 Etna 的修改底图版本正常工作。必须为底图的 PROJ_LIB 错误添加导出路径：
export PROJ_LIB=/path/to/your/instalation/of/anaconda/share/proj/
（Basemap import error in PyCharm —— KeyError: 'PROJ_LIB' 处的解决方案）
并在尝试绘制时出现设置属性错误。它也发生在使用底图教程 (https://basemaptutorial.readthedocs.io/en/latest/plotting_data.html#plot) 时，但不同之处在于，尽管存在错误，但数据的散点仍然绘制为地图顶部的图层。使用 OSM tiles，无法让数据层显示在地图顶部。必须单独导出每个图层，然后使用图像编辑软件进行组合。
编辑：OpenStreetMap 声明他们的图块服务器不能免费使用，并且遵循使用政策：
https://operations.osmfoundation.org/policies/tiles/
请在使用示例之前阅读此内容。
由于我在 Python 3.8 中实现代码时遇到问题，我将一些答案合并在一起并修改了代码。现在它对我有用，我没有收到任何错误。
当我在 Python 3 中尝试 运行 来自 BerndGit 的原始代码时，我不得不进行与他的回答中描述的连接点相同的更改。我替换了
import urllib2
import StringIO
和
import requests
from io import BytesIO
因为 urllib2 库不再适用于 Python 3。您必须使用 urllib.request 或请求。
然后我不得不从 getImageCluster 函数中更改这两行
imgstr = urllib2.urlopen(imgurl).read()
tile = Image.open(StringIO.StringIO(imgstr))
至
imgstr = requests.get(imgurl)
tile = Image.open(BytesIO(imgstr.content))
之后我可以 运行 代码没有错误，但它仍然无法下载图像。结果我总是得到一块黑色的瓷砖。通过一些研究我了解到，在使用请求时伪造用户代理很重要，因为网站可以判断请求来自 Python 并可能阻止它。以下网站对此进行了描述：
https://www.scrapehero.com/how-to-fake-and-rotate-user-agents-using-python-3/
所以我按照网站的建议在 getImageCluster 函数的开头添加了这一行：
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"}
现在我们需要将这些 headers 包含到请求调用中：
imgstr = requests.get(imgurl, headers=headers)
现在整个代码如下所示：
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()
