使 Matplotlib 地图图彼此对齐答案

【问题标题】：Make Matplotlib map plots line up with each other使 Matplotlib 地图图彼此对齐
【发布时间】：2015-11-12 21:09:41
【问题描述】：

我正在尝试在 Python 中生成一个将显示（除其他外）的图形：

A) 从墨卡托投影图像转换的底图

B) 带标签的网格线

我希望人物处于横向墨卡托（或其他球面）投影中。

我已经尝试过 Matplotlib 底图和 Cartopy。 Cartopy 可以做（A），Basemap 可以做（B），但是 Cartopy 只能在 PlateCarree 图上标注网格线，而且 Basemap 不支持使用imshow() 转换图像。

除非有人可以提出其他替代方案，否则我认为解决此问题的最简单方法是在重新投影的图像上覆盖 Basemap 图中的网格线和标签。但是我无法让这两个情节相互对齐。到目前为止我所拥有的：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import cartopy.crs as ccrs

#Setup figure
fig = plt.figure(figsize=(20, 20))

#Set figure limits (lat long)
xlimits = [25, 42]
ylimits = [25, 40]
#Where the projection is centred
centre = [33, 33]

#Image limits in Mercator Eastings and Northings
imxlimits = [25, 42]
imylimits = [25, 40]
#Transform image Limits
imextent = tuple(ccrs.Mercator().transform_points(ccrs.Geodetic(),
    np.array(imxlimits),  np.array(imylimits))[:, 0:2].T.flatten())
#Load image
image = plt.imread(Dir + 'topo.png')

tm = ccrs.TransverseMercator(central_longitude=centre[0], central_latitude=centre[1])
ll = ccrs.Geodetic()

#Setup image axies
ax = fig.add_subplot(111, projection=tm)
ax.set_extent(xlimits + ylimits, ll)          #<--Limits defined here

#Plot image transformed
ax.imshow(image, origin='upper', extent=imextent, transform=ccrs.Mercator())

#Create axes for gridlines
axl = fig.add_subplot(111)
#Make the figure background transparent
axl.patch.set_alpha(0)

#Make basemap instance
m = Basemap(projection='tmerc', resolution='h', ax=axl,
    lat_0=centre[1], lon_0=centre[0],
    llcrnrlon=xlimits[0], llcrnrlat=ylimits[0],    #<--Limits defined here
    urcrnrlon=xlimits[1], urcrnrlat=ylimits[1])

m.drawcoastlines() #to check if the images match up
#Draw gridlines
m.drawparallels(np.arange(20, 50), labels=[False, True, False, False])
m.drawmeridians(np.arange(20, 50), labels=[False, False, False, True])

plt.show()

这会产生大致彼此重叠但不匹配的图。我认为这是因为为第一个图给出的限制可能是针对顶部和底部边缘设置的，而为第二个图给出的（相同）限制是针对右上角和左下角的。

关于如何解决此问题的任何提示？

谢谢！

【问题讨论】：

好问题！也许添加一个cartopy标签？我在将 MPL.Basemap 投影与其他数据对齐时也遇到了问题。我一直认为我没有选择正确的球体或投影参数是我的错误，但我开始觉得 Basemap 在这个部门有一些错误（例如this）。
好的，很高兴听到不仅仅是我有问题。我猜他们可能正在使用不同的方法来转换为轴坐标。

标签： python matplotlib matplotlib-basemap cartopy

【解决方案1】：

所以我设法通过替换 Cartopy 限制的定义方式来使其工作，而不是：

ax.set_extent(xlimits + ylimits, ll)

我先转换成横墨卡托坐标：

tmlims = tm.transform_points(ll, np.array(xlimits),  np.array(ylimits))

然后将范围设置为：

ax.set_extent([tmlims[0][0], tmlims[1][0], tmlims[0][1], tmlims[1][1]], tm)

由于绘图已经在横向墨卡托坐标中，我不知道为什么手动进行转换会有所不同（因为定义 ll 而不是 tm 应该进行相同的转换）。也许我误解了它是如何转换限制的。

在寻找此解决方案的过程中，我还发现了这些代码位，这些代码可能会帮助解决其他“地图不对齐”问题的人：

ax.get_extent(crs) #Return the extents of the Cartopy axies in the given co-ordinates
#
#These work with any matplotlib axies
ax.set_anchor('W') #Anchor the plot to the left (in this case) of the canvas 
ax.set_aspect('equal') #Distort (or not) the aspects of the plot

见http://matplotlib.org/api/axes_api.html和http://scitools.org.uk/cartopy/docs/v0.4/matplotlib/geoaxes.html

【讨论】：