Merge pull request #2069 from greglucas/pyproj-geodesic-trace

Pyproj geodesic trace

.github/workflows/ci-testing.yml

@@ -30,12 +30,12 @@ jobs:
         if: matrix.python-version == '3.8' && matrix.os == 'ubuntu-latest'
         id: minimum-packages
         run: |
-          echo "PACKAGES=cython=0.29.15 matplotlib-base=3.2.1 numpy=1.19 owslib=0.19.1 pyproj=3.0 proj=8.0 scipy=1.4.0 shapely=1.6.4" >> $GITHUB_ENV
+          echo "PACKAGES=cython=0.29.15 matplotlib-base=3.2.1 numpy=1.19 owslib=0.19.1 pyproj=3.0 scipy=1.4.0 shapely=1.6.4" >> $GITHUB_ENV
 
       - name: Latest packages
         if: steps.minimum-packages.conclusion == 'skipped'
         run: |
-          echo "PACKAGES=cython fiona matplotlib-base numpy pyproj 'proj>=8' pykdtree scipy shapely" >> $GITHUB_ENV
+          echo "PACKAGES=cython fiona matplotlib-base numpy pyproj pykdtree scipy shapely" >> $GITHUB_ENV
 
       - name: Coverage packages
         id: coverage

INSTALL

@@ -70,7 +70,7 @@ _______________
 
 On Debian, you can install the required system libraries using the system package manager::
 
-    sudo apt -y install libproj-dev libgeos-dev
+    sudo apt -y install libgeos-dev
 
 Then you can proceed to install cartopy using a python package manager::
 
@@ -81,7 +81,7 @@ _____
 
 For macOS, the required dependencies can be installed in the following way::
 
-    brew install proj geos
+    brew install geos
     pip3 install --upgrade pyshp
     # shapely needs to be built from source to link to geos. If it is already
     # installed, uninstall it by: pip3 uninstall shapely
@@ -124,9 +124,6 @@ Further information about the required dependencies can be found here:
 **pyshp** 2.1 or later (https://pypi.python.org/pypi/pyshp)
     Pure Python read/write support for ESRI Shapefile format.
 
-**PROJ** 8.0.0 or later (https://proj.org/)
-    Coordinate transformation library.
-
 **pyproj** 3.0.0 or later (https://github.com/pyproj4/pyproj/)
     Python interface to PROJ (cartographic projections and coordinate transformations library).
 

environment.yml

@@ -14,8 +14,6 @@ dependencies:
   - geos>=3.7.2
   - pyshp>=2.1
   - pyproj>=3.0.0
-  # PROJ is still required for now
-  - proj>=8
   # The testing label has the proper version of freetype included
   - conda-forge/label/testing::matplotlib-base>=3.1
 

lib/cartopy/crs.py

@@ -35,7 +35,6 @@
 
 __document_these__ = ['CRS', 'Geocentric', 'Geodetic', 'Globe']
 
-PROJ_VERSION = cartopy.trace.PROJ_VERSION
 WGS84_SEMIMAJOR_AXIS = 6378137.0
 WGS84_SEMIMINOR_AXIS = 6356752.3142
 

lib/cartopy/trace.pyx

@@ -49,47 +49,15 @@ cdef extern from "geos_c.h":
     void GEOSPreparedGeom_destroy_r(GEOSContextHandle_t handle, const GEOSPreparedGeometry* g) nogil
     cdef int GEOS_MULTILINESTRING
 
-cdef extern from "geodesic.h":
-    # External imports of Proj4.9 functions
-    cdef struct geod_geodesic:
-        pass
-    cdef struct geod_geodesicline:
-        pass
-
-    void geod_init(geod_geodesic*, double, double) nogil
-    double geod_geninverse(geod_geodesic*, double, double, double, double,
-                           double*, double*, double*, double*, double*,
-                           double*, double*) nogil
-    void geod_lineinit(geod_geodesicline*, geod_geodesic*, double, double,
-                       double, int) nogil
-    void geod_genposition(geod_geodesicline*, int, double, double*,
-                          double*, double*, double*, double*, double*,
-                          double*, double*) nogil
-
-    cdef int GEOD_ARCMODE
-    cdef int GEOD_LATITUDE
-    cdef int GEOD_LONGITUDE
-
 import re
 import warnings
 
 import shapely.geometry as sgeom
 from shapely.geos import lgeos
-from pyproj import Transformer, proj_version_str
+from pyproj import Geod, Transformer
 from pyproj.exceptions import ProjError
 
 
-_match = re.search(r"\d+\.\d+.\d+", proj_version_str)
-if _match is not None:
-    PROJ_VERSION = tuple(int(v) for v in _match.group().split('.'))
-    if PROJ_VERSION < (8, 0, 0):
-        warnings.warn(
-            "PROJ 8+ is required. Current version: {}".format(proj_version_str)
-        )
-else:
-    PROJ_VERSION = ()
-
-
 cdef GEOSContextHandle_t get_geos_context_handle():
     cdef ptr handle = lgeos.geos_handle
     return <GEOSContextHandle_t>handle
@@ -251,9 +219,9 @@ cdef class CartesianInterpolator(Interpolator):
 
 
 cdef class SphericalInterpolator(Interpolator):
-    cdef geod_geodesic geod
-    cdef geod_geodesicline geod_line
-    cdef double a13
+    cdef object geod
+    cdef double azim
+    cdef double s12
 
     cdef void init(self, src_crs, dest_crs) except *:
         self.transformer = Transformer.from_crs(src_crs, dest_crs, always_xy=True)
@@ -262,23 +230,16 @@ cdef class SphericalInterpolator(Interpolator):
         cdef double flattening = 0
         if src_crs.ellipsoid.inverse_flattening > 0:
             flattening = 1 / src_crs.ellipsoid.inverse_flattening
-        geod_init(&self.geod, major_axis, flattening)
+        self.geod = Geod(a=major_axis, f=flattening)
 
     cdef void set_line(self, const Point &start, const Point &end):
-        cdef double azi1
-        self.a13 = geod_geninverse(&self.geod,
-                                   start.y, start.x, end.y, end.x,
-                                   NULL, &azi1, NULL, NULL, NULL, NULL, NULL)
-        geod_lineinit(&self.geod_line, &self.geod, start.y, start.x, azi1,
-                      GEOD_LATITUDE | GEOD_LONGITUDE);
+        Interpolator.set_line(self, start, end)
+        self.azim, _, self.s12 = self.geod.inv(start.x, start.y, end.x, end.y)
 
     cdef Point interpolate(self, double t) except *:
         cdef Point lonlat
 
-        geod_genposition(&self.geod_line, GEOD_ARCMODE, self.a13 * t,
-                         &lonlat.y, &lonlat.x, NULL, NULL, NULL, NULL, NULL,
-                         NULL)
-
+        lonlat.x, lonlat.y, _ = self.geod.fwd(self.start.x, self.start.y, self.azim, self.s12 * t)
         return self.project(lonlat)
 
 

setup.py

@@ -61,8 +61,6 @@
 
 # Please keep in sync with INSTALL file.
 GEOS_MIN_VERSION = (3, 7, 2)
-PROJ_MIN_VERSION = (8, 0, 0)
-PROJ_MIN_VERSION_STRING = '.'.join(str(v) for v in PROJ_MIN_VERSION)
 
 
 def file_walk_relative(top, remove=''):
@@ -115,99 +113,6 @@ def file_walk_relative(top, remove=''):
             geos_libraries.append(entry[2:])
 
 
-# Proj
-def find_proj_version_by_program(conda=None):
-    proj = shutil.which('proj')
-    if proj is None:
-        print(
-            f'Proj {PROJ_MIN_VERSION_STRING} or later must be installed.',
-            file=sys.stderr)
-        exit(1)
-
-    if conda is not None and conda not in proj:
-        print(
-            f'Proj {PROJ_MIN_VERSION_STRING} or later must be installed in'
-            f' Conda environment "{conda}".',
-            file=sys.stderr)
-        exit(1)
-
-    try:
-        proj_version = subprocess.check_output([proj],
-                                               stderr=subprocess.STDOUT)
-        proj_version = proj_version.split()[1].split(b'.')
-        proj_version = tuple(int(v.strip(b',')) for v in proj_version)
-    except (OSError, IndexError, ValueError, subprocess.CalledProcessError):
-        warnings.warn(
-            f'Unable to determine Proj version. Ensure you have '
-            f'{PROJ_MIN_VERSION_STRING} or later installed, or installation '
-            f'may fail.'
-        )
-        proj_version = (0, 0, 0)
-
-    return proj_version
-
-
-conda = os.getenv('CONDA_DEFAULT_ENV')
-if conda is not None and conda in sys.prefix:
-    # Conda does not provide pkg-config compatibility, but the search paths
-    # should be set up so that nothing extra is required. We'll still check
-    # the version, though.
-    proj_version = find_proj_version_by_program(conda)
-    if proj_version < PROJ_MIN_VERSION:
-        proj_version_string = '.'.join(str(v) for v in proj_version)
-        print(
-            f'Proj version {proj_version_string} is installed, but cartopy'
-            f' requires at least version {PROJ_MIN_VERSION_STRING}',
-            file=sys.stderr)
-        exit(1)
-
-    proj_includes = []
-    proj_libraries = ["proj"]
-    proj_library_dirs = []
-
-else:
-    try:
-        proj_version = subprocess.check_output(['pkg-config', '--modversion',
-                                                'proj'],
-                                               stderr=subprocess.STDOUT)
-        proj_version = tuple(int(v) for v in proj_version.split(b'.'))
-        proj_includes = subprocess.check_output(['pkg-config', '--cflags',
-                                                 'proj'])
-        proj_clibs = subprocess.check_output(['pkg-config', '--libs', 'proj'])
-    except (OSError, ValueError, subprocess.CalledProcessError):
-        proj_version = find_proj_version_by_program()
-        if proj_version < PROJ_MIN_VERSION:
-            print(
-                'Proj version %s is installed, but cartopy requires at least '
-                'version %s.' % ('.'.join(str(v) for v in proj_version),
-                                 '.'.join(str(v) for v in PROJ_MIN_VERSION)),
-                file=sys.stderr)
-            exit(1)
-
-        proj_includes = []
-        proj_libraries = ["proj"]
-        proj_library_dirs = []
-    else:
-        if proj_version < PROJ_MIN_VERSION:
-            print(
-                'Proj version %s is installed, but cartopy requires at least '
-                'version %s.' % ('.'.join(str(v) for v in proj_version),
-                                 '.'.join(str(v) for v in PROJ_MIN_VERSION)),
-                file=sys.stderr)
-            exit(1)
-
-        proj_includes = [
-            proj_include[2:] if proj_include.startswith('-I') else
-            proj_include for proj_include in proj_includes.decode().split()]
-
-        proj_libraries = []
-        proj_library_dirs = []
-        for entry in proj_clibs.decode().split():
-            if entry.startswith('-L'):
-                proj_library_dirs.append(entry[2:])
-            elif entry.startswith('-l'):
-                proj_libraries.append(entry[2:])
-
 # Python dependencies
 extras_require = {}
 for name in os.listdir(os.path.join(HERE, 'requirements')):
@@ -253,9 +158,9 @@ def get_config_var(name):
         'cartopy.trace',
         ['lib/cartopy/trace.pyx'],
         include_dirs=([include_dir, './lib/cartopy', np.get_include()] +
-                      proj_includes + geos_includes),
-        libraries=proj_libraries + geos_libraries,
-        library_dirs=[library_dir] + proj_library_dirs + geos_library_dirs,
+                      geos_includes),
+        libraries=geos_libraries,
+        library_dirs=[library_dir] + geos_library_dirs,
         language='c++',
         **extra_extension_args),
 ]
@@ -336,8 +241,6 @@ def decythonize(extensions, **_ignore):
                   ['io/srtm.npz']},
 
     scripts=['tools/cartopy_feature_download.py'],
-
-    # requires proj headers
     ext_modules=extensions,
     cmdclass=cmdclass,
     python_requires='>=' + '.'.join(str(n) for n in PYTHON_MIN_VERSION),