Skip to content
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

Proposal for faster transform option #2104

Open
greglucas opened this issue Nov 17, 2022 · 0 comments
Open

Proposal for faster transform option #2104

greglucas opened this issue Nov 17, 2022 · 0 comments

Comments

@greglucas
Copy link
Contributor

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:
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

1 participant