在 cartopy 中匹配 shapefile 的投影 [英] Match projection of shapefile in cartopy

问题描述

我正在尝试使用matplotlib和cartopy制作Choropleth贴图，显然我首先需要为其绘制shapefile.但是，我没有设法这样做，即使有人问过类似的问题 新的 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 中的第一个设置:

在[34]中:市政[0].bounds出[34]: (32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775)

但是其他元素有不同的界限.您可以基于所有 municipalities 范围的最小/最大值将限制刷新到实际图形.

  bnd = np.array([i.i.bounds for i inboundities 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)

I am trying to make a Choropleth map using matplotlib and cartopy for which I obviously need to plot a shapefile first. However, I did not manage to do so, even though a similar question has been asked here and here. I suspect either the projection or the bounds to be misspecified.

My shapefile has the projection

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]]

and can be downloaded here where I am talking about vg250_2010-01-01.utm32w.shape.ebenen/vg250_ebenen-historisch/de1001/vg250_gem.shp

My code is

#!/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()

where I obtained the bounds using the corresponding method from fiona. Python throws an error

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]

This doesn't make sense to me. Also, experimenting with ccrs.UTM() gives me a plot showing a white area. I'd appreciate it if anyone can tell me how to fix this. Thank you!

解决方案

I have found two issues. One is an incorrect specification of limits in your call to set_extent, the documentation specifies [x0,x1,y0,y1] should be the input, you seem to have given [x0,y0,x1,y1]. The other issue seems to be a limitation in cartopy, as best I can tell. It looks like projections outside the limits listed in the error message will always fail, and those limits are hardcoded. You can just edit the source (this line in their latest release), changing -2e7 to -4e7, likewise for the upper bound. After these fixes, your plot is generated without issue: The new set_extent line:

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

You may also want to set central_longitude=9.0 in your TransverseMercator, that seems to be what's specified in your shapefile.

I would recommend contacting the developers about this, they might have a good reason for setting those bounds, or they might have a better workaround, or maybe they'll widen the bounds in a later release!

Update

Your bounds also seem to have been set based on only the first of the municipalities:

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

But the other elements have different bounds. You can get limits flushed to the actual drawing based on min/max values of the bounds of all municipalities.

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 的投影的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！