-
Notifications
You must be signed in to change notification settings - Fork 359
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
GeoAxes.set_extent on non-cylindrical projections #697
Comments
I agree the behaviour may not be immediately obvious. One advantage of keeping the axes rectangular is that we can make use of the pan/zoom functionality given to us by matplotlib. If we took the extent provided, and used that as the axes extent limits explicitly, we are precluding that functionality. As for the example you gave, I agree - we can make this a huge amount easier with very little effort. For the record, your example can be simplified if we just make use of matplotlib's Path machinery (much of which was added to for cartopy). Your example drops out to be:
|
Proposal:
Because of argument renaming, this would be a breaking change. |
Since cartopy has changed significantly since this issue was opened, I want to clarify a few things.
cartopy/lib/cartopy/mpl/geoaxes.py Lines 857 to 858 in 2148e13
which leads to: Lines 206 to 207 in 2148e13
which leads to: Line 627 in d12c86c
The problem behind #1367 is that the interpolator doesn't use enough points: The issue behind #1362 and an issue I came across with transverse Mercator is that interpolation of the bounding geometry when a Plate Carree CRS is passed does not lead to an extent in projected space that most people expect. Note that [-180, 180, -90, 90] does not work, perhaps because of a numerical issue. When I pass a region in Plate Carre corresponding to the globe, I expect a plot in projected space with as much as possible of the globe. The issue behind #1617 is indeed that no interpolation is done. |
Continuing here #1804. Code to reproduceimport matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy.crs as ccrs
import cartopy.mpl.ticker as ctk
lon1, lon2, lat1, lat2 = [-20, 20, 50, 80]
rect = mpath.Path([[lon1, lat1], [lon2, lat1],
[lon2, lat2], [lon1, lat2], [lon1, lat1]]).interpolated(50)
name='NearsidePerspective'
proj=ccrs.NearsidePerspective(central_longitude=(lon1+lon2)*0.5,
central_latitude=(lat1+lat2)*0.5)
fig, (ax1, ax2) = plt.subplots(1,2, subplot_kw={'projection':proj,})
proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax1) - ax1.transData
rect_in_target = proj_to_data.transform_path(rect)
ax1.set_boundary(rect_in_target)
ax1.set_extent([lon1, lon2, lat1, lat2], crs=ccrs.PlateCarree())
ax1.coastlines()
gl=ax1.gridlines(draw_labels=True, x_inline=False, y_inline=False, linestyle='dashed')
gl.top_labels=False
gl.right_labels=False
gl.rotate_labels=False
gl.xlocator=ctk.LongitudeLocator(4)
gl.ylocator=ctk.LatitudeLocator(6)
gl.xformatter=ctk.LongitudeFormatter(zero_direction_label=True)
gl.yformatter=ctk.LatitudeFormatter()
ax1.set_title(name+'\nset_extent()')
proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax2) - ax2.transData
rect_in_target = proj_to_data.transform_path(rect)
ax2.set_boundary(rect_in_target)
ax2.set_xlim(rect_in_target.vertices[:,0].min(), rect_in_target.vertices[:,0].max())
ax2.set_ylim(rect_in_target.vertices[:,1].min(), rect_in_target.vertices[:,1].max())
ax2.coastlines()
gl=ax2.gridlines(draw_labels=True, x_inline=False, y_inline=False, linestyle='dashed')
gl.top_labels=False
gl.right_labels=False
gl.rotate_labels=False
gl.xlocator=ctk.LongitudeLocator(4)
gl.ylocator=ctk.LatitudeLocator(6)
gl.xformatter=ctk.LongitudeFormatter(zero_direction_label=True)
gl.yformatter=ctk.LatitudeFormatter()
ax2.set_title(name+'\nset_xlim()/set_ylim()')
fig.tight_layout()
plt.show() |
Concerning your example #1367 , again this kind of behaviour is better addressed using Code to reproduceimport numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
def fix_extent(ax, extent):
mlon = np.mean(extent[:2])
mlat = np.mean(extent[2:])
xtrm_data = np.array([[extent[0], mlat], [mlon, extent[2]], [extent[1], mlat], [mlon, extent[3]]])
proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
xtrm = proj_to_data.transform(xtrm_data)
ax.set_xlim(xtrm[:,0].min(), xtrm[:,0].max())
ax.set_ylim(xtrm[:,1].min(), xtrm[:,1].max())
extent=(-180, 180, -60, 60)
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
ax.gridlines()
fix_extent(ax, extent)
plt.show() |
I disagree---your work and #1367 reveal a bug in |
After spending more time delving into this cartopy/lib/cartopy/mpl/geoaxes.py Lines 836 to 839 in 82a6b66
with: x1, x2, y1, y2 = extents
rect = mpath.Path([[x1, y1], [x2, y1],
[x2, y2], [x1, y2], [x1, y1]], closed=True).interpolated(50)
domain_in_crs = cpatch.path_to_geos(rect)[0] This change gives the same result as my previous workaround but using the usual For the sake of completeness I would like to add one more information. The request of a closed |
@pgf Your solution using I agree with @kdpenner that the real solution would probably be to make Failing that, I'm completely fine with the approach to use Matplotlib's |
@dopplershift with the caveat that I know no cython/c++ and haven't used GEOS directly: the important parameter may be Lines 446 to 447 in d12c86c
Line 614 in d12c86c
so a first pass at this would be to alter the projection's |
@kdpenner @dopplershift I've already experimented with the Line 449 in d12c86c
Reducing from 1.0e-6 to 1.0e-8 results in some more points, but not enough to fix this issue:Original value, 1.0e-6 :1.0e-8 :Smaller values do not add new points unfortunatelly. |
I have been looking into
So reducing |
@blazing216 yup, that's my understanding as well. Can you tell if In this thread we are talking about projecting a cartopy/lib/cartopy/mpl/geoaxes.py Lines 837 to 839 in f019487
with a |
@kdpenner Just to clarify, for the Robinson projection I've changed: Lines 1874 to 1877 in 7c75f5a
to: self._threshold = 1e4
@property
def threshold(self):
return self._threshold and then modified #1367 with: print(ax.projection._threshold)
ax.projection._threshold=1.e-2
print(ax.projection._threshold)
ax.set_extent(extent, ccrs.PlateCarree()) Trying several values didn't fix this particular case. |
@pgf thanks for changing My development branch fixes the issue: and I haven't changed |
@kdpenner This makes sense. Indeed sounds like my previous experience:
Maybe the issues is related to differences between |
@kdpenner I think so. Under the hood, both
See the below the beginning of
|
OK, to clarify, Lines 641 to 645 in 7c75f5a
The reason that an extent defined as a Lines 351 to 352 in 4ca93c1
Oops, I looked into #1797 and found a different way to fix the problem. I will follow up in that thread. |
I found cartopy/lib/cartopy/mpl/geoaxes.py Lines 857 to 863 in 7c75f5a
I added the following lines before calling The advantage of the new It also does not break the pan/zoom functionality. But we cannot see the content beyond the I think we have to think about this: what do we want if crs is not None:
path, = cpatch.geos_to_path(projected)
self.patch.set_boundary(path, self.transData)
self.spines['geo'].set_boundary(path, self.transData) Code to produce the example: import cartopy.crs as ccrs
import cartopy.feature as feat
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import shapely.geometry as sgeom
import cartopy.mpl.patch as cpatch
import matplotlib.patches as patches
high_res_proj = ccrs.LambertConformal()
high_res_proj.threshold = 1e3
xlim = [-120, -60]
ylim = [60, 80]
# with smaller threshold, it is no need to interpolate
rect = mpath.Path([[xlim[0], ylim[0]],
[xlim[1], ylim[0]],
[xlim[1], ylim[1]],
[xlim[0], ylim[1]],
[xlim[0], ylim[0]],
])#.interpolated(20)
plt.figure(figsize=(6*2,2*2))
ax = plt.subplot(131, projection=high_res_proj)
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
facecolor=feat.COLORS['land'],
edgecolor='black',
linewidth=1.2))
proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)
ax.add_patch(patches.PathPatch(rect_in_target, edgecolor='r',
facecolor='none', lw=4, zorder=2))
ax.set_boundary(rect_in_target)
# Notice the ugly hack to stop any further clipping - this is
# the same problem as #363.
ax.set_extent([xlim[0], xlim[1], ylim[0] - 2, ylim[1]])
plt.title('Path')
ax = plt.subplot(132, projection=high_res_proj)
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
facecolor=feat.COLORS['land'],
edgecolor='black',
linewidth=1.2))
ax.add_patch(patches.PathPatch(rect_in_target, edgecolor='r',
facecolor='none', lw=4, zorder=2))
ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]])
plt.title('New set_extent(crs=None)')
ax = plt.subplot(133, projection=high_res_proj)
ax.gridlines()
ax.add_feature(feat.NaturalEarthFeature('physical', 'land', '50m',
facecolor=feat.COLORS['land'],
edgecolor='black',
linewidth=1.2))
ax.add_patch(patches.PathPatch(rect_in_target, edgecolor='r',
facecolor='none', lw=4, zorder=2))
ax.set_extent([xlim[0], xlim[1], ylim[0], ylim[1]],
crs=ccrs.PlateCarree())
plt.title('New set_extent(crs=PlateCarree())')
plt.tight_layout()
plt.savefig('boundary_using_set_extent.png', dpi=150) |
@blazing216 where do you add the code above? It breaks a script of mine if I add it here: cartopy/lib/cartopy/mpl/geoaxes.py Line 871 in 7c75f5a |
Yes, I added to where you pointed. I am just sharing an idea that you can add |
I think I am not the only one who finds it a bit confusing that while
GeoAxes.set_extent
accepts acrs
, it only really has an effect on the "corners". That is, it essentially re-projects the corners and then takes the largest x/y boundary that encases them. It is probably fine in something like Mercator or PlateCarree where latitude and longitude are at right-angles.But if you pick something conical like Albers or Lambert Conformal, then lines of constant latitude or longitude are not parallel to the figure edges:
There is a relevant example for how to set a more complex boundary, but there are various subtleties. For example, you cannot just create a rectangle because of #363. I can add this snippet to go along with the star example, but I'm wondering if this should be available as a builtin method?
The text was updated successfully, but these errors were encountered: