设置范围时忽略投影限制

Question

有没有办法设置图形超出投影限制的范围？

例如，当使用 "Rijksdriehoek" 投影 (EPSG 28992) 时，Cartopy (proj4?) 的限制是错误的，太窄了。

该投影旨在覆盖整个荷兰，但强加的限制甚至导致部分国家被切断。而我宁愿将范围设置得比官方边界稍宽，以提供一些额外的上下文。

很遗憾，set_extent方法报错：

ValueError: Failed to determine the required bounds in projection 
coordinates. Check that the values provided are within the valid range 
(x_limits=[646.3608848793374, 284347.25011780026], 
y_limits=[308289.55751689477, 637111.0245778429]).

set_xlim/set_ylim 方法似乎没有做任何事情，这适用于普通的 matplotlib 轴。

示例代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

projection = ccrs.epsg(28992)

fig, ax = plt.subplots(figsize=(5,10), subplot_kw=dict(projection=projection))

ax.coastlines(resolution='10m')
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land', '10m', facecolor='none', edgecolor='k'))

图形范围自动设置为投影范围：

print(projection.bounds)
print(ax.get_extent())

(646.3608848793374, 284347.25011780026, 308289.55751689477, 637111.0245778429)
(646.3608848793374, 284347.25011780026, 308289.55751689477, 637111.0245778429)

根据有关投影的文档，实际限制应为：(-700 300000 289000 629000)。但就可视化目的而言，即使是那些似乎也有点严格。

参见 "Scope of validity section" 的示例：

https://translate.google.com/translate?hl=en&sl=nl&u=https://nl.wikipedia.org/wiki/Rijksdriehoeksco%25C3%25B6rdinaten

Answer 1

我发现 Cartopy 中的投影限制取自 Proj4 中的投影限制，因此无法立即修复。但是，您可以通过询问参数来定义等效投影... 首先，

>>> import pyepsg
>>> proj4_epsg = pyepsg.get(28992)
>>> print(proj4_epsg.as_proj4())
'+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs'
>>>

然后，例如..

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat
proj_equivalent = ccrs.Stereographic(central_longitude=5.3876388888, central_latitude=52.15616055555,
    false_easting=155000, false_northing=463000, scale_factor=0.9999079)
ax = plt.axes(projection=proj_equivalent)
x0, x1 = -4.7e4, +3.7e5
y0, y1 = 2.6e5, 6.82e5
ax.set_extent((x0, x1, y0, y1), crs=proj_equivalent)
ax.coastlines('50m', color='blue'); ax.gridlines()
ax.add_feature(cfeat.BORDERS, edgecolor='red', linestyle='--')
plt.show()

这会产生这样的情节：

很明显，这里的内置县界非常简陋。另外，我还没有设置正确的椭圆，这需要更多的研究。但它展示了如何突破提供的投影边界的限制。

不知道这里有没有机会反击Proj4？

Answer 2

@pp-mo 的回答非常好。但是，这里有一个替代解决方案。工作代码是：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# subclassing for a modified projection
class nd_prj(ccrs.Stereographic):
    """
    Hacked projection for ccrs.epsg(28992) to get extended plotting limits
    """
    def __init__(self):
        globe = ccrs.Globe(ellipse=u'bessel')
        super(nd_prj, self).__init__(central_latitude=52.15616055555555, \
                     central_longitude=5.38763888888889, \
                     #true_scale_latitude=52.0, \
                     scale_factor=0.9999079, \
                     false_easting=155000, false_northing=463000, globe=globe)

    @property
    def x_limits(self):
        return (500, 300000)   # define the values you need

    @property
    def y_limits(self):
        return (300000, 650000) # define the values you need

projection = nd_prj()  # make use of the projection
fig, ax = plt.subplots(figsize=(5,10), subplot_kw=dict(projection=projection))

ax.coastlines(resolution='10m')
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land', \
                                            '10m', facecolor='none', edgecolor='k'))
plt.show()

结果图：

希望这有用。

Answer 3

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

这里是 "custom extent" 投影 class 的稍微灵活一点的版本。这也应该使其适用于其他预测。例如，对于横跨赤道的国家/地区的 UTM 投影。范围仍然需要手动输入，它可以扩展为按百分比扩展默认 proj4 范围。

class ProjectCustomExtent(ccrs.Projection):

    def __init__(self, epsg=28992, extent=[-200000, 500000, 200000, 700000]):

        xmin, xmax, ymin, ymax = extent

        self.xmin = xmin
        self.xmax = xmax
        self.ymin = ymin
        self.ymax = ymax

        super().__init__(ccrs.epsg(epsg).proj4_params)

    @property
    def boundary(self):

        coords = ((self.x_limits[0], self.y_limits[0]),
                  (self.x_limits[0], self.y_limits[1]),
                  (self.x_limits[1], self.y_limits[1]),
                  (self.x_limits[1], self.y_limits[0]))

        return ccrs.sgeom.LineString(coords)

    @property
    def bounds(self):
        xlim = self.x_limits
        ylim = self.y_limits
        return xlim[0], xlim[1], ylim[0], ylim[1]

    @property
    def threshold(self):
        return 1e5

    @property
    def x_limits(self):
        return self.xmin, self.xmax

    @property
    def y_limits(self):
        return self.ymin, self.ymax

获取新投影：

projection = ProjectCustomExtent(epsg=28992, extent=[-300000, 500000, -100000, 800000])

结果：

fig, ax = plt.subplots(figsize=(10,15), subplot_kw=dict(projection=projection), facecolor='w')

ax.coastlines(resolution='10m')
ax.add_feature(cfeature.NaturalEarthFeature('cultural', 'admin_0_boundary_lines_land', '10m', 
                                            facecolor='none', edgecolor='k'), label='Stereo', zorder=999, lw=1, linestyle='-')


ax.set_extent([-100000, 400000, 200000, 700000], crs=projection)

Answer 4

@Rutger_Kassies 的回答非常好，但是对于最新版本的 Cartopy 和 PyProj，该解决方案存在一些问题，我已经修复并简化了它：

import cartopy.crs as ccrs
from cartopy.crs import Projection
from pyproj import CRS

class ProjectCustomExtent(Projection):
    def __init__(self, epsg, extent):
        self.xmin, self.xmax, self.ymin, self.ymax = extent
        super().__init__(CRS.from_epsg(epsg).to_string())

    @Projection.boundary.getter
    def boundary(self):
        coords = ((self.x_limits[0], self.y_limits[0]),
                  (self.x_limits[0], self.y_limits[1]),
                  (self.x_limits[1], self.y_limits[1]),
                  (self.x_limits[1], self.y_limits[0]))
        return ccrs.sgeom.LineString(coords)

    @Projection.x_limits.getter
    def x_limits(self):
        return self.xmin, self.xmax

    @Projection.y_limits.getter
    def y_limits(self):
        return self.ymin, self.ymax

# example
projection = ProjectCustomExtent(epsg=3005, extent=[xmin, xmax, ymin, ymax])

设置范围时忽略投影限制

Ignore projection limits when setting the extent

python

matplotlib

cartopy