cartopy icon indicating copy to clipboard operation
cartopy copied to clipboard

Proposal for faster transform option

Open greglucas opened this issue 3 years ago • 0 comments
trafficstars

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.

  1. 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.
  2. 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.
  3. 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

image

proposed (no interpolation along paths): Figure_1

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:

greglucas avatar Nov 17 '22 23:11 greglucas