3

我正在尝试在 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()

覆盖不匹配的图

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

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

谢谢！

4

1 回答 1

0

因此，我设法通过替换 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

于 2015-08-20T16:17:14.480 回答