People often complain about how slow Cartopy is lately since we moved over to pyproj and not the direct C bindings we were using. If a user knows what they are doing and wants something fast or for a small domain with no wrapped coordinates, we can help them out and give a much faster option, so perhaps we should do that. Below is a quick proposal for how this could work.
Feedback on this would be most welcome! This is just a thought at this point, but perhaps one that others have better ideas for.
Proposal
Add a fast-path to CRS.project_geometry()
that immediately transforms all of the points in the geometry to the destination projection. This is what basemap does and it is much faster than handling interpolations/bisections between points to make smoother lines. It is equivalent to the fast_transform
argument that we added to the contour()
functions.
There are a few ways we could implement this and I'm not sure which would be best / most useful.
- A global configuration attribute somewhere:
cartopy.interpolate_paths = False
.
Benefit: easy to change behavior
Downside: can't individually control which paths/projections are affected by this.
- Add more
fast_transform
keyword decorators to plotting functions.
Benefit: we already have this written
Downside: would require some more work on which x/y arguments to transform for each plotting function, and would also probably require some work on the add_patch()
etc... areas of the code that aren't directly plotting functions.
- Add a new keyword / property to CRS and subclasses that could control this on a per-CRS basis.
Benefit: user could decide which projection and artists to interpolate or not (i.e. you could have both
pc_slow = PlateCarree(interpolate_paths=True) and
pc_fast = PlateCarree(interpolate_paths=False)`)
Downside: Would require more work on keeping CRS's consistent and subclassing of them would add another somewhat arbitrary keyword that isn't CRS-specific, bur rather Cartopy-specific.
Issues with this approach
There are some downsides to this associated with long paths and wrapped coordinates. For example, the "logo" example:
current (with interpolated paths): https://scitools.org.uk/cartopy/docs/latest/gallery/miscellanea/logo.html
proposed (no interpolation along paths):
Code to reproduce
Example of some minimal code that could be adapted/expanded upon.
diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py
index 30c4170e..a77dbb4f 100644
--- a/lib/cartopy/crs.py
+++ b/lib/cartopy/crs.py
@@ -19,6 +19,7 @@ import math
import warnings
import numpy as np
+import shapely.ops
import shapely.geometry as sgeom
from pyproj import Transformer
from pyproj.exceptions import ProjError
@@ -123,7 +124,7 @@ class CRS(_CRS):
#: Whether this projection can handle ellipses.
_handles_ellipses = True
- def __init__(self, proj4_params, globe=None):
+ def __init__(self, proj4_params, globe=None, interpolate_paths=True):
"""
Parameters
----------
@@ -137,7 +138,11 @@ class CRS(_CRS):
globe: :class:`~cartopy.crs.Globe` instance, optional
If omitted, the default Globe instance will be created.
See :class:`~cartopy.crs.Globe` for details.
-
+ interpolate_paths: bool, default True
+ Whether to interpolate paths in this CRS. Setting this
+ to False can significantly improve performance, but
+ paths will not wrap and may look jagged if there
+ are not enough points in the initial geometry.
"""
# for compatibility with pyproj.CRS and rasterio.crs.CRS
try:
@@ -187,6 +192,7 @@ class CRS(_CRS):
init_items.append(f'+{k}')
self.proj4_init = ' '.join(init_items) + ' +no_defs'
super().__init__(self.proj4_init)
+ self.interpolate_paths = interpolate_paths
def __eq__(self, other):
if isinstance(other, CRS) and self.proj4_init == other.proj4_init:
@@ -801,6 +807,11 @@ class Projection(CRS, metaclass=ABCMeta):
elif not isinstance(src_crs, CRS):
raise TypeError('Source CRS must be an instance of CRS'
' or one of its subclasses, or None.')
+
+ if not getattr(src_crs, 'interpolate_paths', True):
+ pyproj_trans = _get_transformer_from_crs(src_crs, self).transform
+ return shapely.ops.transform(pyproj_trans, geometry)
+
geom_type = geometry.geom_type
method_name = self._method_map.get(geom_type)
if not method_name:
Type: Question