使 Matplotlib 地图绘图相互对齐

Question

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

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

B) 标记的网格线

我希望图形在横轴墨卡托（或其他球面）投影中。

Matplotlib Basemap 和 Cartopy 我都试过了。 Cartopy可以做（A），Basemap可以做（B），但是Cartopy只能在PlateCarree图上标注网格线，Basemap不支持使用imshow().

进行图片转换

除非有人可以提出另一种选择，否则我认为解决此问题的最简单方法是将底图绘图中的网格线和标签叠加在重新投影的图像上。但是我无法让这两个地块相互对齐。我目前拥有的：

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()

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

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

谢谢！

Answer 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 应该进行相同的转换）。也许我误解了它是如何转换限制的。

在寻找这个解决方案的过程中，我还发现了这些代码，它们可能会帮助人们解决其他 'Maps not lining up' 问题：

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 and http://scitools.org.uk/cartopy/docs/v0.4/matplotlib/geoaxes.html

使 Matplotlib 地图绘图相互对齐

Make Matplotlib map plots line up with each other

python

matplotlib

matplotlib-basemap

cartopy