在cartopy中匹配shapefile的投影

Question

我正在尝试使用 matplotlib 和 cartopy 制作 Choropleth 地图，我显然需要先为其绘制 shapefile。然而，我没有设法这样做，即使有人问过类似的问题here and here。我怀疑投影或边界指定有误。

我的 shapefile 有投影

PROJCS["WGS_1984_UTM_Zone_32Nz",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",32500000],
    PARAMETER["False_Northing",0],
    PARAMETER["Central_Meridian",9],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0],
    UNIT["Meter",1]]

可以下载here我说的是vg250_2010-01-01.utm32w.shape.ebenen/vg250_ebenen-historisch/de1001/vg250_gem.shp

我的密码是

#!/usr/local/bin/python
# -*- coding: utf8 -*-

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

fname = 'path/vg250_gem.shp'
proj = ccrs.TransverseMercator(central_longitude=0.0,central_latitude=0.0,
                           false_easting=32500000.0,false_northing=0.0,
                           scale_factor=0.9996)
municipalities = list(shpreader.Reader(fname).geometries())
ax = plt.axes(projection=proj)
plt.title('Deutschland')
ax.add_geometries(municipalities,proj,edgecolor='black',facecolor='gray',alpha=0.5)
ax.set_extent([32458044.649189778*0.9, 5556418.748046352*1.1, 32465287.307457082*0.9, 5564153.5456742775*1.1],proj)
plt.show()

我使用 fiona 的相应方法获得了边界。 Python 抛出错误

Traceback (most recent call last):
  File "***/src/analysis/test.py", line 16, in <module>
ax.set_extent([32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775],proj)
  File "/usr/local/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py", line 652, in set_extent
ylim=self.projection.y_limits))
ValueError: Failed to determine the required bounds in projection coordinates. Check that the values provided are within the valid range (x_limits=[-20000000.0, 20000000.0], y_limits=[-10000000.0, 10000000.0]).
[Finished in 53.9s with exit code 1]

这对我来说没有意义。此外，使用 ccrs.UTM() 进行试验会给我一个显示白色区域的图。如果有人能告诉我如何解决这个问题，我将不胜感激。谢谢！

Answer 1

我发现了两个问题。一个是您在调用 set_extent 时指定的限制不正确，documentation 指定 [x0,x1,y0,y1] 应该是输入，您似乎已经给出了 [x0,y0,x1,y1]。据我所知，另一个问题似乎是 cartopy 的限制。看起来超出错误消息中列出的限制的投影总是会失败，并且这些限制是硬编码的。您可以只编辑源 (this line in their latest release)，将 -2e7 更改为 -4e7，对于上限也是如此。在这些修复之后，您的情节生成没有问题：新的 set_extent 行：

ax.set_extent([32458044.649189778*0.975, 32465287.307457082*1.025,5556418.748046352*0.9, 556415,3.5456742775*1.1],proj)

您可能还想在 TransverseMercator 中设置 central_longitude=9.0，这似乎是您的 shapefile 中指定的内容。

我建议就此联系开发人员，他们可能有设置这些界限的充分理由，或者他们可能有更好的解决方法，或者他们可能会在以后的版本中扩大界限！

更新

您的界限似乎也仅根据 municipalities 中的第一个设置：

In [34]: municipalities[0].bounds
Out[34]: (32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775)

但其他元素有不同的界限。您可以根据所有 municipalities.

边界的 min/max 值将限制刷新到实际图形中

bnd = np.array([i.bounds for i in municipalities])
x0,x1 = np.min(bnd[:,0]),np.max(bnd[:,2])
y0,y1 = np.min(bnd[:,1]),np.max(bnd[:,3])
ax.set_extent([x0,x1,y0,y1],proj)

在cartopy中匹配shapefile的投影

Match projection of shapefile in cartopy

projection

matplotlib

cartopy