设置范围时忽略投影限制

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)。但即使是那些看起来也有点严格的可视化目的。

参见例如“有效性范围”部分：

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

这里是“自定义范围”投影类的稍微灵活的版本。这也应该使它适用于其他预测。例如，在一个跨越赤道的国家的 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])