Existe-t-il un moyen de définir l'étendue de la figure au-delà des limites de projection?
Par exemple, lors de l'utilisation de la projection "Rijksdriehoek" (EPSG 28992), les limites de Cartopy (proj4?) sont également erronées étroit.
Cette projection est conçue pour couvrir tous les Pays-Bas, mais les limites imposées entraînent même la coupure d'une partie du pays. Alors que je préfère définir l'étendue légèrement plus large que les limites officielles pour fournir un contexte supplémentaire.
Malheureusement, le set_extent renvoie une erreur:
print(projection.bounds) print(ax.get_extent()) (646.3608848793374, 284347.25011780026, 308289.55751689477, 637111.0245778429) (646.3608848793374, 284347.25011780026, 308289.55751689477, 637111.0245778429)
Les méthodes set_xlim / set_ylim ne semblent rien faire , qui fonctionnerait pour un axe matplotlib normal.
Exemple de code:
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'))
L'étendue de la figure est automatiquement définie aux limites de la projection:
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]).
Selon la documentation sur la projection, les limites réelles devraient être: (-700 300000 289000 629000) . Mais même ceux-ci semblent un peu stricts à des fins de visualisation.
Voir par exemple la section "Périmètre de validité" sur:
3 Réponses :
J'ai trouvé que les limites de projection dans Cartopy sont tirées de celles de Proj4, donc il n'y a pas de solution immédiate. Cependant, vous pouvez définir une projection équivalente en interrogeant les paramètres ... Premièrement,
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()
Ensuite, par exemple ..
>>> 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' >>>
Ce qui donne un tracé comme celui-ci:
De toute évidence, les limites de comté intégrées sont très grossières ici. De plus, je n'ai pas configuré l'ellipse correcte, ce qui nécessite un peu plus de recherche. Mais cela montre comment briser la limitation des limites de projection fournies.
Je ne sais pas s'il y a une opportunité ici de repousser Proj4?
La réponse de @ pp-mo est excellente. Cependant, voici une solution alternative. Le code de travail est:
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()
Le tracé résultant:
J'espère que cela vous sera utile .
Merci! Après quelques recherches sur Github, je suis également arrivé à la conclusion qu'une projection personnalisée est probablement le meilleur moyen de gérer cela. Voir ma réponse ci-dessous pour une légère variation à ce sujet. Initialiser la projection avec ccrs.epsg (28992) est très bien dans ce cas, puisque la définition elle-même est correcte, mais les limites ne le sont pas.
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)
Question stupide: à l'échelle de cette limitation (c'est-à-dire une carte de tout le pays), observeriez-vous une différence avec toute autre projection qui n'a pas ces limites?
Oui, cela serait bien sûr possible, en utilisant par exemple
UTM 31N. Mais il faudrait également reprojeter les données raster que j'ai, ce qui est moins optimal. C'est une solution de contournement.Encore une fois, je me trompe peut-être ici, mais j'avais à l'esprit que la différence entre les systèmes de coordonnées serait de 100 mètres au maximum, donc une reprojection pourrait même ne pas être nécessaire? (Cela est basé sur Gauß-Krüger contre UMT en Allemagne, je ne sais pas pour "Rijksdriehoek")