Игнорировать пределы проекции при настройке экстента

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

Ответ @ @ 1001 * превосходен.Однако вот альтернативное решение.Рабочий код:

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 2

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 3

Я обнаружил, что пределы проекции в 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?