通过在适当的缩放级别上从Web地图瓦片构建基础地图是获得更好质量的绘图的一种方法。在这里,我演示如何从openstreetmap
Web地图服务器中获取它们。在此示例中,我使用缩放级别10,并获取2个地图瓦片以组合为单个图像数组。其中一个缺点是,组合图像的范围总是大于我们要求的值。以下是可用的代码:
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):
'''Lon./lat. to tile numbers'''
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):
'''Tile numbers to lon./lat.'''
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
def getextents(latmin_deg, lonmin_deg, delta_lat, delta_long, zoom):
'''Return LL and UR, each with (long,lat) of real extent of combined tiles.
latmin_deg: bottom lat of extent
lonmin_deg: left long of extent
delta_lat: extent of lat
delta_long: extent of long, all in degrees
'''
xtile_LL, ytile_LL = deg2num(latmin_deg, lonmin_deg, zoom)
xtile_UR, ytile_UR = deg2num(latmin_deg + delta_lat, lonmin_deg + delta_long, zoom)
lat_NW_LL, lon_NW_LL = num2deg(xtile_LL, ytile_LL, zoom)
lat_NW_LLL, lon_NW_LLL = num2deg(xtile_LL, ytile_LL+1, zoom)
lat_NW_UR, lon_NW_UR = num2deg(xtile_UR, ytile_UR, zoom)
lat_NW_URR, lon_NW_URR = num2deg(xtile_UR+1, ytile_UR, zoom)
minLat = lat_NW_LLL
minLon = lon_NW_LL
maxLat = lat_NW_UR
maxLon = lon_NW_URR
return (minLon, maxLon, minLat, maxLat)
wl = -74.04006
sl = 40.683092
el = -73.834067
nl = 40.88378
lat_deg = sl
lon_deg = wl
d_lat = nl - sl
d_long = el - wl
zoom = 10
timg = getImageCluster(lat_deg, lon_deg, d_lat, d_long, zoom)
latmin_deg, lonmin_deg, delta_lat, delta_long = sl, wl, nl-sl, el-wl
(left, right, bottom, top) = getextents(latmin_deg, lonmin_deg, delta_lat, delta_long, zoom)
m = Basemap(resolution='h',
projection='merc',
area_thresh=50,
lat_0=(bottom + top)/2, lon_0=(left + right)/2,
llcrnrlon=left, llcrnrlat=bottom, urcrnrlon=right, urcrnrlat=top)
fig = plt.figure()
fig.set_size_inches(10, 12)
m.imshow(np.asarray(timg), extent=[left, right, bottom, top], origin='upper' )
m.drawcoastlines(color='gray', linewidth=3.0)
plt.show()
希望对您有所帮助。生成的图表如下:
编辑
为了裁剪图像以得到要绘制的精确区域并不困难,可以使用
PIL
模块进行处理,Numpy的数组切片也可以起到作用。