From cc309ee2f7945ad76d47439d15a8c6e746d1fbb3 Mon Sep 17 00:00:00 2001
From: Elliott Sales de Andrade <quantum.analyst@gmail.com>
Date: Sat, 29 Dec 2018 18:58:14 -0500
Subject: [PATCH 01/10] Use shapely for prepared geometry handling.

---
 lib/cartopy/trace.pyx | 79 ++++++++++++++++---------------------------
 1 file changed, 30 insertions(+), 49 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index f33dbd55b..376164015 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -32,29 +32,24 @@ cdef extern from "geos_c.h":
     ctypedef struct GEOSGeometry:
         pass
     ctypedef struct GEOSCoordSequence
-    ctypedef struct GEOSPreparedGeometry
-    GEOSContextHandle_t GEOS_init_r() nogil
     GEOSCoordSequence *GEOSCoordSeq_create_r(GEOSContextHandle_t, unsigned int, unsigned int) nogil
-    GEOSGeometry *GEOSGeom_createPoint_r(GEOSContextHandle_t, GEOSCoordSequence *) nogil
     GEOSGeometry *GEOSGeom_createLineString_r(GEOSContextHandle_t, GEOSCoordSequence *) nogil
     GEOSGeometry *GEOSGeom_createCollection_r(GEOSContextHandle_t, int, GEOSGeometry **, unsigned int) nogil
     GEOSGeometry *GEOSGeom_createEmptyCollection_r(GEOSContextHandle_t, int) nogil
-    void GEOSGeom_destroy_r(GEOSContextHandle_t, GEOSGeometry *) nogil
     GEOSCoordSequence *GEOSGeom_getCoordSeq_r(GEOSContextHandle_t, GEOSGeometry *) nogil
     int GEOSCoordSeq_getSize_r(GEOSContextHandle_t handle, const GEOSCoordSequence* s, unsigned int *size) nogil
     int GEOSCoordSeq_getX_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double *) nogil
     int GEOSCoordSeq_getY_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double *) nogil
     int GEOSCoordSeq_setX_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
     int GEOSCoordSeq_setY_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
-    const GEOSPreparedGeometry *GEOSPrepare_r(GEOSContextHandle_t handle, const GEOSGeometry* g) nogil
-    char GEOSPreparedCovers_r(GEOSContextHandle_t, const GEOSPreparedGeometry*, const GEOSGeometry*) nogil
-    char GEOSPreparedDisjoint_r(GEOSContextHandle_t, const GEOSPreparedGeometry*, const GEOSGeometry*) nogil
-    void GEOSPreparedGeom_destroy_r(GEOSContextHandle_t handle, const GEOSPreparedGeometry* g) nogil
     cdef int GEOS_MULTILINESTRING
 
 import re
 import warnings
 
+import shapely.geometry as sgeom
+import shapely.prepared as sprep
+from shapely.geos import lgeos
 from pyproj import Geod, Transformer
 from pyproj.exceptions import ProjError
 import shapely.geometry as sgeom
@@ -252,22 +247,14 @@ cdef enum State:
     POINT_NAN
 
 
-cdef State get_state(const Point &point, const GEOSPreparedGeometry *gp_domain,
-                     GEOSContextHandle_t handle):
+cdef State get_state(const Point &point, object gp_domain):
     cdef State state
-    cdef GEOSCoordSequence *coords
-    cdef GEOSGeometry *g_point
 
     if isfinite(point.x) and isfinite(point.y):
         # TODO: Avoid create-destroy
-        coords = GEOSCoordSeq_create_r(handle, 1, 2)
-        GEOSCoordSeq_setX_r(handle, coords, 0, point.x)
-        GEOSCoordSeq_setY_r(handle, coords, 0, point.y)
-        g_point = GEOSGeom_createPoint_r(handle, coords)
-        state = (POINT_IN
-                 if GEOSPreparedCovers_r(handle, gp_domain, g_point)
-                 else POINT_OUT)
-        GEOSGeom_destroy_r(handle, g_point)
+        g_point = sgeom.Point((point.x, point.y))
+        state = POINT_IN if gp_domain.covers(g_point) else POINT_OUT
+        del g_point
     else:
         state = POINT_NAN
     return state
@@ -278,7 +265,7 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
                             double t_end, const Point &p_end,
                             Interpolator interpolator, double threshold,
                             GEOSContextHandle_t handle,
-                            const GEOSPreparedGeometry *gp_domain,
+                            object gp_domain,
                             bool inside) except *:
     """
     Return whether the given line segment is suitable as an
@@ -305,8 +292,6 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
     cdef double along
     cdef double separation
     cdef double hypot
-    cdef GEOSCoordSequence *coords
-    cdef GEOSGeometry *g_segment
 
     # This could be optimised out of the loop.
     if not (isfinite(p_start.x) and isfinite(p_start.y)):
@@ -385,26 +370,23 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
             # TODO: Re-use geometries, instead of create-destroy!
 
             # Create a LineString for the current end-point.
-            coords = GEOSCoordSeq_create_r(handle, 2, 2)
-            GEOSCoordSeq_setX_r(handle, coords, 0, p_start.x)
-            GEOSCoordSeq_setY_r(handle, coords, 0, p_start.y)
-            GEOSCoordSeq_setX_r(handle, coords, 1, p_end.x)
-            GEOSCoordSeq_setY_r(handle, coords, 1, p_end.y)
-            g_segment = GEOSGeom_createLineString_r(handle, coords)
+            g_segment = sgeom.LineString([
+                (p_start.x, p_start.y),
+                (p_end.x, p_end.y)])
 
             if inside:
-                valid = GEOSPreparedCovers_r(handle, gp_domain, g_segment)
+                valid = gp_domain.covers(g_segment)
             else:
-                valid = GEOSPreparedDisjoint_r(handle, gp_domain, g_segment)
+                valid = gp_domain.disjoint(g_segment)
 
-            GEOSGeom_destroy_r(handle, g_segment)
+            del g_segment
 
     return valid
 
 
 cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
                  GEOSContextHandle_t handle,
-                 const GEOSPreparedGeometry *gp_domain, const State &state,
+                 object gp_domain, const State &state,
                  Interpolator interpolator, double threshold,
                  double &t_min, Point &p_min, double &t_max, Point &p_max) except *:
     cdef double t_current
@@ -460,7 +442,7 @@ cdef void _project_segment(GEOSContextHandle_t handle,
                            const GEOSCoordSequence *src_coords,
                            unsigned int src_idx_from, unsigned int src_idx_to,
                            Interpolator interpolator,
-                           const GEOSPreparedGeometry *gp_domain,
+                           object gp_domain,
                            double threshold, LineAccumulator lines) except *:
     cdef Point p_current, p_min, p_max, p_end
     cdef double t_current, t_min, t_max
@@ -484,7 +466,7 @@ cdef void _project_segment(GEOSContextHandle_t handle,
         print("   ", p_end.x, ", ", p_end.y)
 
     t_current = 0.0
-    state = get_state(p_current, gp_domain, handle)
+    state = get_state(p_current, gp_domain)
 
     cdef size_t old_lines_size = lines.size()
     while t_current < 1.0 and (lines.size() - old_lines_size) < 100:
@@ -517,7 +499,7 @@ cdef void _project_segment(GEOSContextHandle_t handle,
             else:
                 t_current = t_max
                 p_current = p_max
-                state = get_state(p_current, gp_domain, handle)
+                state = get_state(p_current, gp_domain)
                 if state == POINT_IN:
                     lines.new_line()
 
@@ -528,14 +510,14 @@ cdef void _project_segment(GEOSContextHandle_t handle,
             else:
                 t_current = t_max
                 p_current = p_max
-                state = get_state(p_current, gp_domain, handle)
+                state = get_state(p_current, gp_domain)
                 if state == POINT_IN:
                     lines.new_line()
 
         else:
             t_current = t_max
             p_current = p_max
-            state = get_state(p_current, gp_domain, handle)
+            state = get_state(p_current, gp_domain)
             if state == POINT_IN:
                 lines.new_line()
 
@@ -582,19 +564,19 @@ def project_linear(geometry not None, src_crs not None,
         GEOSContextHandle_t handle = get_geos_context_handle()
         GEOSGeometry *g_linear = geos_from_shapely(geometry)
         Interpolator interpolator
-        GEOSGeometry *g_domain
+        object g_domain
         const GEOSCoordSequence *src_coords
         unsigned int src_size, src_idx
-        const GEOSPreparedGeometry *gp_domain
+        object gp_domain
         LineAccumulator lines
         GEOSGeometry *g_multi_line_string
 
-    g_domain = geos_from_shapely(dest_projection.domain)
+    g_domain = dest_projection.domain
 
     interpolator = _interpolator(src_crs, dest_projection)
 
     src_coords = GEOSGeom_getCoordSeq_r(handle, g_linear)
-    gp_domain = GEOSPrepare_r(handle, g_domain)
+    gp_domain = sprep.prep(g_domain)
 
     GEOSCoordSeq_getSize_r(handle, src_coords, &src_size)  # check exceptions
 
@@ -603,7 +585,7 @@ def project_linear(geometry not None, src_crs not None,
         _project_segment(handle, src_coords, src_idx - 1, src_idx,
                          interpolator, gp_domain, threshold, lines);
 
-    GEOSPreparedGeom_destroy_r(handle, gp_domain)
+    del gp_domain
 
     g_multi_line_string = lines.as_geom(handle)
 
@@ -617,7 +599,7 @@ class _Testing:
     def straight_and_within(Point l_start, Point l_end,
                             double t_start, double t_end,
                             Interpolator interpolator, double threshold,
-                            domain):
+                            object domain):
         # This function is for testing/demonstration only.
         # It is not careful about freeing resources, and it short-circuits
         # optimisations that are made in the real algorithm (in exchange for
@@ -625,11 +607,10 @@ class _Testing:
 
         cdef GEOSContextHandle_t handle = get_geos_context_handle()
 
-        cdef GEOSGeometry *g_domain = geos_from_shapely(domain)
-        cdef const GEOSPreparedGeometry *gp_domain
-        gp_domain = GEOSPrepare_r(handle, g_domain)
+        cdef object gp_domain
+        gp_domain = sprep.prep(domain)
 
-        state = get_state(interpolator.project(l_start), gp_domain, handle)
+        state = get_state(interpolator.project(l_start), gp_domain)
         cdef bool p_start_inside_domain = state == POINT_IN
 
         # l_end and l_start should be un-projected.
@@ -643,7 +624,7 @@ class _Testing:
             interpolator, threshold,
             handle, gp_domain, p_start_inside_domain)
 
-        GEOSPreparedGeom_destroy_r(handle, gp_domain)
+        del gp_domain
         return valid
 
     @staticmethod

From fe84165b8c9dff96aee5881c0e88fca023aa0917 Mon Sep 17 00:00:00 2001
From: Elliott Sales de Andrade <quantum.analyst@gmail.com>
Date: Sat, 29 Dec 2018 19:06:48 -0500
Subject: [PATCH 02/10] Use shapely when grabbing projecting coordinates.

---
 lib/cartopy/trace.pyx | 24 +++++++-----------------
 1 file changed, 7 insertions(+), 17 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index 376164015..350392c0a 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -60,12 +60,6 @@ cdef GEOSContextHandle_t get_geos_context_handle():
     return handle
 
 
-cdef GEOSGeometry *geos_from_shapely(shapely_geom) except *:
-    """Get the GEOS pointer from the given shapely geometry."""
-    cdef ptr geos_geom = shapely_geom._geom
-    return <GEOSGeometry *>geos_geom
-
-
 cdef shapely_from_geos(GEOSGeometry *geom):
     """Turn the given GEOS geometry pointer into a shapely geometry."""
     return sgeom.base.geom_factory(<ptr>geom)
@@ -439,8 +433,7 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
 
 
 cdef void _project_segment(GEOSContextHandle_t handle,
-                           const GEOSCoordSequence *src_coords,
-                           unsigned int src_idx_from, unsigned int src_idx_to,
+                           tuple src_from, tuple src_to,
                            Interpolator interpolator,
                            object gp_domain,
                            double threshold, LineAccumulator lines) except *:
@@ -448,10 +441,8 @@ cdef void _project_segment(GEOSContextHandle_t handle,
     cdef double t_current, t_min, t_max
     cdef State state
 
-    GEOSCoordSeq_getX_r(handle, src_coords, src_idx_from, &p_current.x)
-    GEOSCoordSeq_getY_r(handle, src_coords, src_idx_from, &p_current.y)
-    GEOSCoordSeq_getX_r(handle, src_coords, src_idx_to, &p_end.x)
-    GEOSCoordSeq_getY_r(handle, src_coords, src_idx_to, &p_end.y)
+    p_current.x, p_current.y = src_from
+    p_end.x, p_end.y = src_to
     if DEBUG:
         print("Setting line:")
         print("   ", p_current.x, ", ", p_current.y)
@@ -562,10 +553,9 @@ def project_linear(geometry not None, src_crs not None,
     cdef:
         double threshold = dest_projection.threshold
         GEOSContextHandle_t handle = get_geos_context_handle()
-        GEOSGeometry *g_linear = geos_from_shapely(geometry)
         Interpolator interpolator
         object g_domain
-        const GEOSCoordSequence *src_coords
+        object src_coords
         unsigned int src_size, src_idx
         object gp_domain
         LineAccumulator lines
@@ -575,14 +565,14 @@ def project_linear(geometry not None, src_crs not None,
 
     interpolator = _interpolator(src_crs, dest_projection)
 
-    src_coords = GEOSGeom_getCoordSeq_r(handle, g_linear)
+    src_coords = geometry.coords
     gp_domain = sprep.prep(g_domain)
 
-    GEOSCoordSeq_getSize_r(handle, src_coords, &src_size)  # check exceptions
+    src_size = len(src_coords)  # check exceptions
 
     lines = LineAccumulator()
     for src_idx in range(1, src_size):
-        _project_segment(handle, src_coords, src_idx - 1, src_idx,
+        _project_segment(handle, src_coords[src_idx - 1], src_coords[src_idx],
                          interpolator, gp_domain, threshold, lines);
 
     del gp_domain

From e68b3c1237fc0299ee703e3db0e00ce9e07bf79a Mon Sep 17 00:00:00 2001
From: Elliott Sales de Andrade <quantum.analyst@gmail.com>
Date: Sun, 30 Dec 2018 18:28:11 -0500
Subject: [PATCH 03/10] Use shapely to create projected geometries.

---
 lib/cartopy/trace.pyx | 26 ++++++--------------------
 1 file changed, 6 insertions(+), 20 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index 350392c0a..626e1cf73 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -98,7 +98,7 @@ cdef class LineAccumulator:
         if self.lines.back().empty():
             self.add_point(point)
 
-    cdef GEOSGeometry *as_geom(self, GEOSContextHandle_t handle):
+    cdef object as_geom(self, GEOSContextHandle_t handle):
         from cython.operator cimport dereference, preincrement
 
         # self.lines.remove_if(degenerate_line) is not available in Cython.
@@ -121,24 +121,12 @@ cdef class LineAccumulator:
 
         cdef Line ilines
         cdef Point ipoints
-        cdef vector[GEOSGeometry *] geoms
-        cdef int i
-        cdef GEOSCoordSequence *coords
+        geoms = []
         for ilines in self.lines:
-            coords = GEOSCoordSeq_create_r(handle, ilines.size(), 2)
-            for i, ipoints in enumerate(ilines):
-                GEOSCoordSeq_setX_r(handle, coords, i, ipoints.x)
-                GEOSCoordSeq_setY_r(handle, coords, i, ipoints.y)
+            coords = [(ipoints.x, ipoints.y) for ipoints in ilines]
+            geoms.append(sgeom.LineString(coords))
 
-            geoms.push_back(GEOSGeom_createLineString_r(handle, coords))
-
-        cdef GEOSGeometry *geom
-        if geoms.empty():
-            geom = GEOSGeom_createEmptyCollection_r(handle,
-                                                    GEOS_MULTILINESTRING)
-        else:
-            geom = GEOSGeom_createCollection_r(handle, GEOS_MULTILINESTRING,
-                                               &geoms[0], geoms.size())
+        geom = sgeom.MultiLineString(geoms)
         return geom
 
     cdef size_t size(self):
@@ -559,7 +547,6 @@ def project_linear(geometry not None, src_crs not None,
         unsigned int src_size, src_idx
         object gp_domain
         LineAccumulator lines
-        GEOSGeometry *g_multi_line_string
 
     g_domain = dest_projection.domain
 
@@ -577,10 +564,9 @@ def project_linear(geometry not None, src_crs not None,
 
     del gp_domain
 
-    g_multi_line_string = lines.as_geom(handle)
+    multi_line_string = lines.as_geom(handle)
 
     del lines, interpolator
-    multi_line_string = shapely_from_geos(g_multi_line_string)
     return multi_line_string
 
 

From 8311eed7e3a1df7c88054f6234beb79083d65ec9 Mon Sep 17 00:00:00 2001
From: Elliott Sales de Andrade <quantum.analyst@gmail.com>
Date: Sun, 30 Dec 2018 18:45:59 -0500
Subject: [PATCH 04/10] Remove remaining GEOS bits.

---
 lib/cartopy/trace.pyx | 58 ++++++++-----------------------------------
 1 file changed, 10 insertions(+), 48 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index 626e1cf73..a27a945a8 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -18,53 +18,22 @@ from functools import lru_cache
 
 cimport cython
 from libc.math cimport HUGE_VAL, sqrt
-from libc.stdint cimport uintptr_t as ptr
+from numpy.math cimport isfinite, isnan
 from libcpp cimport bool
 from libcpp.list cimport list
-from libcpp.vector cimport vector
-from numpy.math cimport isfinite, isnan
-
 
 cdef bool DEBUG = False
 
-cdef extern from "geos_c.h":
-    ctypedef void *GEOSContextHandle_t
-    ctypedef struct GEOSGeometry:
-        pass
-    ctypedef struct GEOSCoordSequence
-    GEOSCoordSequence *GEOSCoordSeq_create_r(GEOSContextHandle_t, unsigned int, unsigned int) nogil
-    GEOSGeometry *GEOSGeom_createLineString_r(GEOSContextHandle_t, GEOSCoordSequence *) nogil
-    GEOSGeometry *GEOSGeom_createCollection_r(GEOSContextHandle_t, int, GEOSGeometry **, unsigned int) nogil
-    GEOSGeometry *GEOSGeom_createEmptyCollection_r(GEOSContextHandle_t, int) nogil
-    GEOSCoordSequence *GEOSGeom_getCoordSeq_r(GEOSContextHandle_t, GEOSGeometry *) nogil
-    int GEOSCoordSeq_getSize_r(GEOSContextHandle_t handle, const GEOSCoordSequence* s, unsigned int *size) nogil
-    int GEOSCoordSeq_getX_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double *) nogil
-    int GEOSCoordSeq_getY_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double *) nogil
-    int GEOSCoordSeq_setX_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
-    int GEOSCoordSeq_setY_r(GEOSContextHandle_t, GEOSCoordSequence *, int, double) nogil
-    cdef int GEOS_MULTILINESTRING
-
 import re
 import warnings
 
 import shapely.geometry as sgeom
 import shapely.prepared as sprep
-from shapely.geos import lgeos
-from pyproj import Geod, Transformer
+from pyproj import Geod, Transformer, proj_version_str
 from pyproj.exceptions import ProjError
 import shapely.geometry as sgeom
 
 
-cdef GEOSContextHandle_t get_geos_context_handle():
-    cdef GEOSContextHandle_t handle = GEOS_init_r()
-    return handle
-
-
-cdef shapely_from_geos(GEOSGeometry *geom):
-    """Turn the given GEOS geometry pointer into a shapely geometry."""
-    return sgeom.base.geom_factory(<ptr>geom)
-
-
 ctypedef struct Point:
     double x
     double y
@@ -98,7 +67,7 @@ cdef class LineAccumulator:
         if self.lines.back().empty():
             self.add_point(point)
 
-    cdef object as_geom(self, GEOSContextHandle_t handle):
+    cdef object as_geom(self):
         from cython.operator cimport dereference, preincrement
 
         # self.lines.remove_if(degenerate_line) is not available in Cython.
@@ -246,7 +215,6 @@ cdef State get_state(const Point &point, object gp_domain):
 cdef bool straightAndDomain(double t_start, const Point &p_start,
                             double t_end, const Point &p_end,
                             Interpolator interpolator, double threshold,
-                            GEOSContextHandle_t handle,
                             object gp_domain,
                             bool inside) except *:
     """
@@ -259,7 +227,6 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
     p_start: Projected end point.
     interpolator: Interpolator for current source line.
     threshold: Lateral tolerance in target projection coordinates.
-    handle: Thread-local context handle for GEOS.
     gp_domain: Prepared polygon of target map domain.
     inside: Whether the start point is within the map domain.
 
@@ -367,7 +334,6 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
 
 
 cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
-                 GEOSContextHandle_t handle,
                  object gp_domain, const State &state,
                  Interpolator interpolator, double threshold,
                  double &t_min, Point &p_min, double &t_max, Point &p_max) except *:
@@ -396,13 +362,13 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
             # Straight and entirely-inside-domain
             valid = straightAndDomain(t_start, p_start, t_current, p_current,
                                       interpolator, threshold,
-                                      handle, gp_domain, True)
+                                      gp_domain, True)
 
         elif state == POINT_OUT:
             # Straight and entirely-outside-domain
             valid = straightAndDomain(t_start, p_start, t_current, p_current,
                                       interpolator, threshold,
-                                      handle, gp_domain, False)
+                                      gp_domain, False)
         else:
             valid = not isfinite(p_current.x) or not isfinite(p_current.y)
 
@@ -420,8 +386,7 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
         p_current = interpolator.interpolate(t_current)
 
 
-cdef void _project_segment(GEOSContextHandle_t handle,
-                           tuple src_from, tuple src_to,
+cdef void _project_segment(tuple src_from, tuple src_to,
                            Interpolator interpolator,
                            object gp_domain,
                            double threshold, LineAccumulator lines) except *:
@@ -461,7 +426,7 @@ cdef void _project_segment(GEOSContextHandle_t handle,
             print("   ", p_current.x, ", ", p_current.y)
             print("   ", p_end.x, ", ", p_end.y)
 
-        bisect(t_current, p_current, p_end, handle, gp_domain, state,
+        bisect(t_current, p_current, p_end, gp_domain, state,
                interpolator, threshold,
                t_min, p_min, t_max, p_max)
         if DEBUG:
@@ -540,7 +505,6 @@ def project_linear(geometry not None, src_crs not None,
     """
     cdef:
         double threshold = dest_projection.threshold
-        GEOSContextHandle_t handle = get_geos_context_handle()
         Interpolator interpolator
         object g_domain
         object src_coords
@@ -559,12 +523,12 @@ def project_linear(geometry not None, src_crs not None,
 
     lines = LineAccumulator()
     for src_idx in range(1, src_size):
-        _project_segment(handle, src_coords[src_idx - 1], src_coords[src_idx],
+        _project_segment(src_coords[src_idx - 1], src_coords[src_idx],
                          interpolator, gp_domain, threshold, lines);
 
     del gp_domain
 
-    multi_line_string = lines.as_geom(handle)
+    multi_line_string = lines.as_geom()
 
     del lines, interpolator
     return multi_line_string
@@ -581,8 +545,6 @@ class _Testing:
         # optimisations that are made in the real algorithm (in exchange for
         # a convenient signature).
 
-        cdef GEOSContextHandle_t handle = get_geos_context_handle()
-
         cdef object gp_domain
         gp_domain = sprep.prep(domain)
 
@@ -598,7 +560,7 @@ class _Testing:
         valid = straightAndDomain(
             t_start, p0, t_end, p1,
             interpolator, threshold,
-            handle, gp_domain, p_start_inside_domain)
+            gp_domain, p_start_inside_domain)
 
         del gp_domain
         return valid

From 59ddff6f55f56de48f6da2170d0278be2e27c4d9 Mon Sep 17 00:00:00 2001
From: Elliott Sales de Andrade <quantum.analyst@gmail.com>
Date: Mon, 31 Dec 2018 02:37:52 -0500
Subject: [PATCH 05/10] Remove GEOS linking.

---
 INSTALL  |  4 ----
 setup.py | 43 +------------------------------------------
 2 files changed, 1 insertion(+), 46 deletions(-)

diff --git a/INSTALL b/INSTALL
index ccc5aee8b..23351841f 100644
--- a/INSTALL
+++ b/INSTALL
@@ -113,10 +113,6 @@ Further information about the required dependencies can be found here:
     Python package for 2D plotting. Python package required for any
     graphical capabilities.
 
-**GEOS** 3.7.2 or later (https://trac.osgeo.org/geos/)
-    GEOS is an API of spatial predicates and functions for processing geometry
-    written in C++.
-
 **Shapely** 1.7.1 or later (https://github.com/Toblerity/Shapely)
     Python package for the manipulation and analysis of planar geometric
     objects.
diff --git a/setup.py b/setup.py
index 783871ac1..4ddda131c 100644
--- a/setup.py
+++ b/setup.py
@@ -56,10 +56,6 @@
     raise ImportError('NumPy 1.21+ is required to install cartopy.')
 
 
-# Please keep in sync with INSTALL file.
-GEOS_MIN_VERSION = (3, 7, 2)
-
-
 def file_walk_relative(root):
     """
     Return a list of files from the top of the tree, removing
@@ -73,40 +69,6 @@ def file_walk_relative(root):
 # Dependency checks
 # =================
 
-# GEOS
-try:
-    geos_version = subprocess.check_output(['geos-config', '--version'])
-    geos_version = tuple(int(v) for v in geos_version.split(b'.')
-                         if 'dev' not in str(v))
-    geos_includes = subprocess.check_output(['geos-config', '--includes'])
-    geos_clibs = subprocess.check_output(['geos-config', '--clibs'])
-except (OSError, ValueError, subprocess.CalledProcessError):
-    warnings.warn(
-        'Unable to determine GEOS version. Ensure you have %s or later '
-        'installed, or installation may fail.' % (
-            '.'.join(str(v) for v in GEOS_MIN_VERSION), ))
-
-    geos_includes = []
-    geos_library_dirs = []
-    geos_libraries = ['geos_c']
-else:
-    if geos_version < GEOS_MIN_VERSION:
-        print('GEOS version %s is installed, but cartopy requires at least '
-              'version %s.' % ('.'.join(str(v) for v in geos_version),
-                               '.'.join(str(v) for v in GEOS_MIN_VERSION)),
-              file=sys.stderr)
-        exit(1)
-
-    geos_includes = geos_includes.decode().split()
-    geos_libraries = []
-    geos_library_dirs = []
-    for entry in geos_clibs.decode().split():
-        if entry.startswith('-L'):
-            geos_library_dirs.append(entry[2:])
-        elif entry.startswith('-l'):
-            geos_libraries.append(entry[2:])
-
-
 # Python dependencies
 extras_require = {}
 for name in (HERE / 'requirements').iterdir():
@@ -151,10 +113,7 @@ def get_config_var(name):
     Extension(
         'cartopy.trace',
         ['lib/cartopy/trace.pyx'],
-        include_dirs=([include_dir, './lib/cartopy', np.get_include()] +
-                      geos_includes),
-        libraries=geos_libraries,
-        library_dirs=[library_dir] + geos_library_dirs,
+        include_dirs=[include_dir, './lib/cartopy', np.get_include()],
         language='c++',
         **extra_extension_args),
 ]

From 0266fcd319ce08d8b2d51b38c87dc9eb39ba3776 Mon Sep 17 00:00:00 2001
From: Greg Lucas <greg.m.lucas@gmail.com>
Date: Wed, 14 Sep 2022 08:21:57 -0600
Subject: [PATCH 06/10] PERF: get source geometry coords as numpy array

This avoids __getitem__ calls from the shapely geometry and converts
to numpy array and memory views within cython.
---
 lib/cartopy/trace.pyx | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index a27a945a8..693012783 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -27,6 +27,7 @@ cdef bool DEBUG = False
 import re
 import warnings
 
+import numpy as np
 import shapely.geometry as sgeom
 import shapely.prepared as sprep
 from pyproj import Geod, Transformer, proj_version_str
@@ -386,7 +387,7 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
         p_current = interpolator.interpolate(t_current)
 
 
-cdef void _project_segment(tuple src_from, tuple src_to,
+cdef void _project_segment(double[:] src_from, double[:] src_to,
                            Interpolator interpolator,
                            object gp_domain,
                            double threshold, LineAccumulator lines) except *:
@@ -507,7 +508,7 @@ def project_linear(geometry not None, src_crs not None,
         double threshold = dest_projection.threshold
         Interpolator interpolator
         object g_domain
-        object src_coords
+        double[:, :] src_coords
         unsigned int src_size, src_idx
         object gp_domain
         LineAccumulator lines
@@ -516,14 +517,14 @@ def project_linear(geometry not None, src_crs not None,
 
     interpolator = _interpolator(src_crs, dest_projection)
 
-    src_coords = geometry.coords
+    src_coords = np.asarray(geometry.coords)
     gp_domain = sprep.prep(g_domain)
 
     src_size = len(src_coords)  # check exceptions
 
     lines = LineAccumulator()
     for src_idx in range(1, src_size):
-        _project_segment(src_coords[src_idx - 1], src_coords[src_idx],
+        _project_segment(src_coords[src_idx - 1, :], src_coords[src_idx, :],
                          interpolator, gp_domain, threshold, lines);
 
     del gp_domain

From 308e9deb162897d9f7259ede0f91a4e414f7d62d Mon Sep 17 00:00:00 2001
From: Greg Lucas <greg.m.lucas@gmail.com>
Date: Wed, 14 Sep 2022 18:09:27 -0600
Subject: [PATCH 07/10] PERF: Project source geometry right away in Cython

This transforms all points in the source geometry right away, rather
than one-by-one later on. This saves some calls to pyproj's transform()
method, which added some overhead.
---
 lib/cartopy/trace.pyx | 53 +++++++++++++++++++++++++++++++++++++++----
 1 file changed, 49 insertions(+), 4 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index 693012783..dea5d5024 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -156,6 +156,47 @@ cdef class Interpolator:
         dest_xy.y = yy * self.dest_scale
         return dest_xy
 
+    cdef double[:, :] project_points(self, double[:, :] src_xy) except *:
+        # Used for fallback to single point updates
+        cdef Point xy
+        # Make a temporary copy so we don't update the incoming memory view
+        new_src_xy = np.asarray(src_xy)*self.src_scale
+        try:
+            xx, yy = self.transformer.transform(
+                new_src_xy[:, 0],
+                new_src_xy[:, 1],
+                errcheck=True
+            )
+        except ProjError as err:
+            msg = str(err).lower()
+            if (
+                "latitude" in msg or
+                "longitude" in msg or
+                "outside of projection domain" in msg or
+                "tolerance condition error" in msg
+            ):
+                # Go back to trying to project a single point at a time
+                xx = np.empty(shape=len(src_xy))
+                yy = np.empty(shape=len(src_xy))
+                for i in range(len(src_xy)):
+                    # Update the point object's x/y coords
+                    xy.x = src_xy[i, 0]
+                    xy.y = src_xy[i, 1]
+                    xy = self.project(xy)
+                    xx[i] = xy.x
+                    yy[i] = xy.y
+            else:
+                raise
+
+        if self.to_180:
+            # Get the places where we should wrap
+            wrap_locs = (xx > 180) | (xx < -180) & (xx != HUGE_VAL)
+            # Do the wrap at those locations
+            xx[wrap_locs] = (((xx[wrap_locs] + 180) % 360) - 180)
+
+        # Destination xy [ncoords, 2]
+        return np.stack([xx, yy], axis=-1) * self.dest_scale
+
     cdef Point interpolate(self, double t) except *:
         raise NotImplementedError
 
@@ -388,6 +429,7 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
 
 
 cdef void _project_segment(double[:] src_from, double[:] src_to,
+                           double[:] dest_from, double[:] dest_to,
                            Interpolator interpolator,
                            object gp_domain,
                            double threshold, LineAccumulator lines) except *:
@@ -403,8 +445,9 @@ cdef void _project_segment(double[:] src_from, double[:] src_to,
         print("   ", p_end.x, ", ", p_end.y)
 
     interpolator.set_line(p_current, p_end)
-    p_current = interpolator.project(p_current)
-    p_end = interpolator.project(p_end)
+    # Now update the current/end with the destination (projected) coords
+    p_current.x, p_current.y = dest_from
+    p_end.x, p_end.y = dest_to
     if DEBUG:
         print("Projected as:")
         print("   ", p_current.x, ", ", p_current.y)
@@ -508,7 +551,7 @@ def project_linear(geometry not None, src_crs not None,
         double threshold = dest_projection.threshold
         Interpolator interpolator
         object g_domain
-        double[:, :] src_coords
+        double[:, :] src_coords, dest_coords
         unsigned int src_size, src_idx
         object gp_domain
         LineAccumulator lines
@@ -518,13 +561,15 @@ def project_linear(geometry not None, src_crs not None,
     interpolator = _interpolator(src_crs, dest_projection)
 
     src_coords = np.asarray(geometry.coords)
+    dest_coords = interpolator.project_points(src_coords)
     gp_domain = sprep.prep(g_domain)
 
     src_size = len(src_coords)  # check exceptions
 
     lines = LineAccumulator()
     for src_idx in range(1, src_size):
-        _project_segment(src_coords[src_idx - 1, :], src_coords[src_idx, :],
+        _project_segment(src_coords[src_idx - 1, :2], src_coords[src_idx, :2],
+                         dest_coords[src_idx - 1, :2], dest_coords[src_idx, :2],
                          interpolator, gp_domain, threshold, lines);
 
     del gp_domain

From c87cc5bed9fdc5ba5f07e5e122b501a390a7f7e4 Mon Sep 17 00:00:00 2001
From: Greg Lucas <greg.m.lucas@gmail.com>
Date: Thu, 8 Dec 2022 11:03:29 -0700
Subject: [PATCH 08/10] PERF: Avoid create/destroy shapely geometry

With Shapely 2.0 we can avoid the create/destroy point for checking
whether a point is within a geometry.
---
 lib/cartopy/trace.pyx | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index dea5d5024..4bcba61af 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -28,6 +28,7 @@ import re
 import warnings
 
 import numpy as np
+import shapely
 import shapely.geometry as sgeom
 import shapely.prepared as sprep
 from pyproj import Geod, Transformer, proj_version_str
@@ -244,10 +245,13 @@ cdef State get_state(const Point &point, object gp_domain):
     cdef State state
 
     if isfinite(point.x) and isfinite(point.y):
-        # TODO: Avoid create-destroy
-        g_point = sgeom.Point((point.x, point.y))
-        state = POINT_IN if gp_domain.covers(g_point) else POINT_OUT
-        del g_point
+        if shapely.__version__ >= "2":
+            # Shapely 2.0 doesn't need to create/destroy a point
+            state = POINT_IN if shapely.intersects_xy(gp_domain.context, point.x, point.y) else POINT_OUT
+        else:
+            g_point = sgeom.Point((point.x, point.y))
+            state = POINT_IN if gp_domain.covers(g_point) else POINT_OUT
+            del g_point
     else:
         state = POINT_NAN
     return state

From 49c6b5a105f2be7d88cdf9e75209e6bfc538b1dc Mon Sep 17 00:00:00 2001
From: Greg Lucas <greg.m.lucas@gmail.com>
Date: Tue, 27 Dec 2022 16:21:06 -0700
Subject: [PATCH 09/10] PERF: Add a fast-path geometry check for rectangular
 destinations

If we are going to a rectangular domain and all initial destination
points are within the domain, then we can avoid checks to see if
each individual line segment is within the domain.
---
 lib/cartopy/trace.pyx | 45 +++++++++++++++++++++++++++++--------------
 1 file changed, 31 insertions(+), 14 deletions(-)

diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx
index 4bcba61af..3cac0cad3 100644
--- a/lib/cartopy/trace.pyx
+++ b/lib/cartopy/trace.pyx
@@ -35,6 +35,8 @@ from pyproj import Geod, Transformer, proj_version_str
 from pyproj.exceptions import ProjError
 import shapely.geometry as sgeom
 
+import cartopy.crs as ccrs
+
 
 ctypedef struct Point:
     double x
@@ -241,9 +243,11 @@ cdef enum State:
     POINT_NAN
 
 
-cdef State get_state(const Point &point, object gp_domain):
+cdef State get_state(const Point &point, object gp_domain, bool geom_fully_inside=False):
     cdef State state
-
+    if geom_fully_inside:
+        # Fast-path return because the geometry is fully inside
+        return POINT_IN
     if isfinite(point.x) and isfinite(point.y):
         if shapely.__version__ >= "2":
             # Shapely 2.0 doesn't need to create/destroy a point
@@ -262,7 +266,8 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
                             double t_end, const Point &p_end,
                             Interpolator interpolator, double threshold,
                             object gp_domain,
-                            bool inside) except *:
+                            bool inside,
+                            bool geom_fully_inside=False) except *:
     """
     Return whether the given line segment is suitable as an
     approximation of the projection of the source line.
@@ -275,6 +280,7 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
     threshold: Lateral tolerance in target projection coordinates.
     gp_domain: Prepared polygon of target map domain.
     inside: Whether the start point is within the map domain.
+    geom_fully_inside: Whether all points are within the map domain.
 
     """
     # Straight and in-domain (de9im[7] == 'F')
@@ -361,7 +367,7 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
                     hypot = mid_dx*mid_dx + mid_dy*mid_dy
                     valid = ((separation * separation) / hypot) < 0.04
 
-        if valid:
+        if valid and not geom_fully_inside:
             # TODO: Re-use geometries, instead of create-destroy!
 
             # Create a LineString for the current end-point.
@@ -382,7 +388,8 @@ cdef bool straightAndDomain(double t_start, const Point &p_start,
 cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
                  object gp_domain, const State &state,
                  Interpolator interpolator, double threshold,
-                 double &t_min, Point &p_min, double &t_max, Point &p_max) except *:
+                 double &t_min, Point &p_min, double &t_max, Point &p_max,
+                 bool geom_fully_inside=False) except *:
     cdef double t_current
     cdef Point p_current
     cdef bool valid
@@ -408,13 +415,13 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end,
             # Straight and entirely-inside-domain
             valid = straightAndDomain(t_start, p_start, t_current, p_current,
                                       interpolator, threshold,
-                                      gp_domain, True)
+                                      gp_domain, True, geom_fully_inside=geom_fully_inside)
 
         elif state == POINT_OUT:
             # Straight and entirely-outside-domain
             valid = straightAndDomain(t_start, p_start, t_current, p_current,
                                       interpolator, threshold,
-                                      gp_domain, False)
+                                      gp_domain, False, geom_fully_inside=geom_fully_inside)
         else:
             valid = not isfinite(p_current.x) or not isfinite(p_current.y)
 
@@ -436,7 +443,8 @@ cdef void _project_segment(double[:] src_from, double[:] src_to,
                            double[:] dest_from, double[:] dest_to,
                            Interpolator interpolator,
                            object gp_domain,
-                           double threshold, LineAccumulator lines) except *:
+                           double threshold, LineAccumulator lines,
+                           bool geom_fully_inside=False) except *:
     cdef Point p_current, p_min, p_max, p_end
     cdef double t_current, t_min, t_max
     cdef State state
@@ -458,7 +466,7 @@ cdef void _project_segment(double[:] src_from, double[:] src_to,
         print("   ", p_end.x, ", ", p_end.y)
 
     t_current = 0.0
-    state = get_state(p_current, gp_domain)
+    state = get_state(p_current, gp_domain, geom_fully_inside)
 
     cdef size_t old_lines_size = lines.size()
     while t_current < 1.0 and (lines.size() - old_lines_size) < 100:
@@ -476,7 +484,7 @@ cdef void _project_segment(double[:] src_from, double[:] src_to,
 
         bisect(t_current, p_current, p_end, gp_domain, state,
                interpolator, threshold,
-               t_min, p_min, t_max, p_max)
+               t_min, p_min, t_max, p_max, geom_fully_inside=geom_fully_inside)
         if DEBUG:
             print("   => ", t_min, "to", t_max)
             print("   => (", p_min.x, ", ", p_min.y, ") to (",
@@ -491,7 +499,7 @@ cdef void _project_segment(double[:] src_from, double[:] src_to,
             else:
                 t_current = t_max
                 p_current = p_max
-                state = get_state(p_current, gp_domain)
+                state = get_state(p_current, gp_domain, geom_fully_inside)
                 if state == POINT_IN:
                     lines.new_line()
 
@@ -502,14 +510,14 @@ cdef void _project_segment(double[:] src_from, double[:] src_to,
             else:
                 t_current = t_max
                 p_current = p_max
-                state = get_state(p_current, gp_domain)
+                state = get_state(p_current, gp_domain, geom_fully_inside)
                 if state == POINT_IN:
                     lines.new_line()
 
         else:
             t_current = t_max
             p_current = p_max
-            state = get_state(p_current, gp_domain)
+            state = get_state(p_current, gp_domain, geom_fully_inside)
             if state == POINT_IN:
                 lines.new_line()
 
@@ -570,11 +578,20 @@ def project_linear(geometry not None, src_crs not None,
 
     src_size = len(src_coords)  # check exceptions
 
+    # Test the entire geometry to see if there are any domain crossings
+    # If there are none, then we can skip expensive domain checks later
+    # TODO: Handle projections other than rectangular
+    cdef bool geom_fully_inside = False
+    if isinstance(dest_projection, (ccrs._RectangularProjection, ccrs._WarpedRectangularProjection)):
+        dest_line = sgeom.LineString([(x[0], x[1]) for x in dest_coords])
+        geom_fully_inside = gp_domain.covers(dest_line)
+
     lines = LineAccumulator()
     for src_idx in range(1, src_size):
         _project_segment(src_coords[src_idx - 1, :2], src_coords[src_idx, :2],
                          dest_coords[src_idx - 1, :2], dest_coords[src_idx, :2],
-                         interpolator, gp_domain, threshold, lines);
+                         interpolator, gp_domain, threshold, lines,
+                         geom_fully_inside=geom_fully_inside);
 
     del gp_domain
 

From 2dcd21d9f0159ee300f534a1dce0577cde786431 Mon Sep 17 00:00:00 2001
From: Greg Lucas <greg.m.lucas@gmail.com>
Date: Sat, 11 Mar 2023 10:42:09 -0700
Subject: [PATCH 10/10] FIX/TST: Use nearby boundary rather than extended
 segments

We don't "project" by a factor of a two to find the boundary
intersection, rather we project along the boundary to the nearest
point. So really we just want to make sure that our cuts are "close"
to the boundary, but we don't care about the segments being extended
by a given fraction.

This is relevant when the final two coordinates are close to each
other and thus that segment is not projected very far, yet it is
quite close to the boundary.
---
 lib/cartopy/tests/test_linear_ring.py | 24 ++++++++++--------------
 1 file changed, 10 insertions(+), 14 deletions(-)

diff --git a/lib/cartopy/tests/test_linear_ring.py b/lib/cartopy/tests/test_linear_ring.py
index f610bde6c..eabdf57a6 100644
--- a/lib/cartopy/tests/test_linear_ring.py
+++ b/lib/cartopy/tests/test_linear_ring.py
@@ -23,25 +23,21 @@ def test_cuts(self):
         assert len(multi_line_string.geoms) > 1
         assert not rings
 
-        def assert_intersection_with_boundary(segment_coords):
-            # Double the length of the segment.
-            start = segment_coords[0]
-            end = segment_coords[1]
-            end = [end[i] + 2 * (end[i] - start[i]) for i in (0, 1)]
-            extended_segment = sgeom.LineString([start, end])
-            # And see if it crosses the boundary.
-            intersection = extended_segment.intersection(projection.boundary)
-            assert not intersection.is_empty, 'Bad topology near boundary'
-
-        # Each line resulting from the split should start and end with a
-        # segment that crosses the boundary when extended to double length.
+        def assert_close_to_boundary(xy):
+            # Are we close to the boundary, which we are considering within
+            # a fraction of the x domain limits
+            limit = (projection.x_limits[1] - projection.x_limits[0]) / 1e4
+            assert sgeom.Point(*xy).distance(projection.boundary) < limit, \
+                'Bad topology near boundary'
+
+        # Each line resulting from the split should be close to the boundary.
         # (This is important when considering polygon rings which need to be
         # attached to the boundary.)
         for line_string in multi_line_string.geoms:
             coords = list(line_string.coords)
             assert len(coords) >= 2
-            assert_intersection_with_boundary(coords[1::-1])
-            assert_intersection_with_boundary(coords[-2:])
+            assert_close_to_boundary(coords[0])
+            assert_close_to_boundary(coords[-1])
 
     def test_out_of_bounds(self):
         # Check that a ring that is completely out of the map boundary