diff --git a/.appveyor.yml b/.appveyor.yml index f5c11f703..cfb916a53 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -2,7 +2,7 @@ environment: matrix: - PYTHON_VERSION: "3.9" CONDA_INSTALL_LOCN: "C:\\Miniconda3-x64" - PACKAGES: "cython numpy matplotlib-base proj=7 pykdtree scipy" + PACKAGES: "cython numpy matplotlib-base pyproj pykdtree scipy" install: # Install miniconda @@ -15,7 +15,7 @@ install: - conda config --add channels conda-forge - conda config --add channels conda-forge/label/testing - set ENV_NAME=test-environment - - set PACKAGES=%PACKAGES% flufl.lock owslib pep8 pillow pyepsg pyshp pytest + - set PACKAGES=%PACKAGES% flufl.lock owslib pep8 pillow pyshp pytest - set PACKAGES=%PACKAGES% requests setuptools_scm setuptools_scm_git_archive - set PACKAGES=%PACKAGES% shapely - conda create -n %ENV_NAME% python=%PYTHON_VERSION% %PACKAGES% @@ -37,6 +37,7 @@ build_script: test_script: - set MPLBACKEND=Agg + - set PYPROJ_GLOBAL_CONTEXT=ON - pytest --pyargs cartopy artifacts: diff --git a/.circleci/config.yml b/.circleci/config.yml index 8dbc8e6a6..3b626d79b 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -36,8 +36,7 @@ deps-run: &deps-install numpy \ owslib \ pillow \ - 'proj<8' \ - pyepsg \ + pyproj \ pykdtree \ pyshp \ requests \ @@ -59,6 +58,7 @@ doc-run: &doc-build name: Build documentation command: | source activate test-environment + PYPROJ_GLOBAL_CONTEXT=ON make -C docs html diff --git a/.github/workflows/ci-testing.yml b/.github/workflows/ci-testing.yml index 976b23be2..2c207c4ee 100644 --- a/.github/workflows/ci-testing.yml +++ b/.github/workflows/ci-testing.yml @@ -10,7 +10,7 @@ jobs: fail-fast: false matrix: os: [ubuntu-latest, macos-latest] - python-version: [3.6, 3.9] + python-version: [3.7, 3.8, 3.9] defaults: run: shell: bash -l {0} @@ -26,16 +26,16 @@ jobs: - name: Minimum packages # Only run on macos for now # Conda's linux packages don't grab the testing label of matplotlib causing failures due to freetype differences - if: matrix.python-version == '3.6' && matrix.os == 'macos-latest' + if: matrix.python-version == '3.7' && matrix.os == 'macos-latest' id: minimum-packages run: | - echo "PACKAGES=cython=0.28.5 matplotlib=2.2.2 numpy=1.16 owslib=0.17 proj4=5.2.0 scipy=1.2.0" >> $GITHUB_ENV + echo "PACKAGES=cython=0.28.5 matplotlib=3.1 numpy=1.18 owslib=0.17 pyproj=3.0 proj=8.0 scipy=1.2.0" >> $GITHUB_ENV echo "CFLAGS=-stdlib=libc++" >> $GITHUB_ENV - name: Latest packages if: steps.minimum-packages.conclusion == 'skipped' run: | - echo "PACKAGES=cython fiona matplotlib-base numpy proj<8 pykdtree scipy" >> $GITHUB_ENV + echo "PACKAGES=cython fiona matplotlib-base numpy pyproj pykdtree scipy" >> $GITHUB_ENV - name: Coverage packages id: coverage @@ -48,7 +48,7 @@ jobs: - name: Install dependencies run: | - PACKAGES="$PACKAGES flufl.lock owslib pep8 pillow pyepsg pyshp pytest" + PACKAGES="$PACKAGES flufl.lock owslib pep8 pillow pyshp pytest" PACKAGES="$PACKAGES pytest-xdist requests setuptools_scm" PACKAGES="$PACKAGES setuptools_scm_git_archive shapely" conda install $PACKAGES @@ -70,7 +70,7 @@ jobs: # Check that the downloader tool at least knows where to get the data from (but don't actually download it) python tools/cartopy_feature_download.py gshhs physical --dry-run CARTOPY_GIT_DIR=$PWD - pytest -n 4 --doctest-modules --pyargs cartopy ${EXTRA_TEST_ARGS} + PYPROJ_GLOBAL_CONTEXT=ON pytest -n 4 --doctest-modules --pyargs cartopy ${EXTRA_TEST_ARGS} - name: Coveralls if: steps.coverage.conclusion == 'success' diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index 1e3d96c15..fd09d9024 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -24,8 +24,8 @@ jobs: - name: Install dependencies run: | - PACKAGES="cython fiona matplotlib-base numpy proj<8 pykdtree scipy" - PACKAGES="$PACKAGES flufl.lock owslib pep8 pillow pyepsg pyshp pytest" + PACKAGES="cython fiona matplotlib-base numpy pyproj pykdtree scipy" + PACKAGES="$PACKAGES flufl.lock owslib pep8 pillow pyshp pytest" PACKAGES="$PACKAGES pytest-xdist requests setuptools_scm" PACKAGES="$PACKAGES setuptools_scm_git_archive shapely" conda install $PACKAGES diff --git a/.gitignore b/.gitignore index 6b5d66c28..42f2e407d 100644 --- a/.gitignore +++ b/.gitignore @@ -53,6 +53,7 @@ benchmarks/envs/ *.swp .ipynb_checkpoints/ .idea/ +.vscode/ # Operating system files \.DS_Store diff --git a/INSTALL b/INSTALL index 7d83e05c0..42d81b44c 100644 --- a/INSTALL +++ b/INSTALL @@ -90,8 +90,10 @@ Further information about the required dependencies can be found here: **pyshp** 2.0 or later (https://pypi.python.org/pypi/pyshp) Pure Python read/write support for ESRI Shapefile format. -**PROJ** 4.9.0 or later (https://proj4.org/) - Cartographic Projections library. +**PROJ** 8.0.0 or later (https://proj.org/). + +**pyproj** 3.0.0 or later (https://github.com/pyproj4/pyproj/) + Python interface to PROJ (cartographic projections and coordinate transformations library). Optional Dependencies ~~~~~~~~~~~~~~~~~~~~~ @@ -109,9 +111,6 @@ to install these optional dependencies. **Pillow** 1.7.8 or later (https://pypi.python.org/pypi/Pillow/2.3.0) A popular fork of PythonImagingLibrary. -**pyepsg** 0.4.0 or later (https://github.com/rhattersley/pyepsg) - A simple Python interface to https://epsg.io - **pykdtree** 1.2.2 or later (https://github.com/storpipfugl/pykdtree) A fast kd-tree implementation that is used for faster warping of images than SciPy. diff --git a/examples/lines_and_polygons/effects_of_the_ellipse.py b/examples/lines_and_polygons/effects_of_the_ellipse.py index dac76685a..f1404c582 100644 --- a/examples/lines_and_polygons/effects_of_the_ellipse.py +++ b/examples/lines_and_polygons/effects_of_the_ellipse.py @@ -46,10 +46,8 @@ def transform_fn(x, y, z=None): def main(): # Define the two coordinate systems with different ellipses. - wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84', - ellipse='WGS84')) - sphere = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84', - ellipse='sphere')) + wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(ellipse='WGS84')) + sphere = ccrs.PlateCarree(globe=ccrs.Globe(ellipse='sphere')) # Define the coordinate system of the data we have from Natural Earth and # acquire the 1:10m physical coastline shapefile. diff --git a/examples/miscellanea/eccentric_ellipse.py b/examples/miscellanea/eccentric_ellipse.py index b6c2d05ed..084ce6996 100644 --- a/examples/miscellanea/eccentric_ellipse.py +++ b/examples/miscellanea/eccentric_ellipse.py @@ -50,8 +50,8 @@ def vesta_image(): img_globe = ccrs.Globe(semimajor_axis=285000., semiminor_axis=229000., ellipse=None) img_proj = ccrs.PlateCarree(globe=img_globe) - img_extent = (-895353.906273091, 895353.906273091, - 447676.9531365455, -447676.9531365455) + img_extent = (-180, 180, + -90, 90) return img, img_globe, img_proj, img_extent @@ -60,7 +60,7 @@ def main(): projection = ccrs.Geostationary(globe=globe) fig = plt.figure() - ax = fig.add_subplot(1, 1, 1, projection=projection) + ax = fig.add_subplot(1, 1, 1, projection=projection, frameon=False) ax.imshow(img, transform=crs, extent=extent) fig.text(.075, .012, "Image credit: NASA/JPL-Caltech/UCLA/MPS/DLR/IDA/PSI", bbox={'facecolor': 'w', 'edgecolor': 'k'}) diff --git a/lib/cartopy/_crs.pxd b/lib/cartopy/_crs.pxd deleted file mode 100644 index 5d106393f..000000000 --- a/lib/cartopy/_crs.pxd +++ /dev/null @@ -1,20 +0,0 @@ -# Copyright Cartopy Contributors -# -# This file is part of Cartopy and is released under the LGPL license. -# See COPYING and COPYING.LESSER in the root of the repository for full -# licensing details. - -from ._proj4 cimport projPJ - - -cdef class CRS: - """ - Defines a Coordinate Reference System using proj. - - """ - - cdef projPJ proj4 - cdef readonly proj4_init - cdef proj4_params - - cpdef is_geodetic(self) diff --git a/lib/cartopy/_crs.pyx b/lib/cartopy/_crs.pyx deleted file mode 100644 index a8b1eb8a1..000000000 --- a/lib/cartopy/_crs.pyx +++ /dev/null @@ -1,620 +0,0 @@ -# Copyright Cartopy Contributors -# -# This file is part of Cartopy and is released under the LGPL license. -# See COPYING and COPYING.LESSER in the root of the repository for full -# licensing details. -# -# cython: embedsignature=True - -""" -This module defines the core CRS class which can interface with Proj. -The CRS class is the base-class for all projections defined in :mod:`cartopy.crs`. - -""" - -from collections import OrderedDict -import re -import warnings - -import numpy as np -cimport numpy as np - -from cython.operator cimport dereference as deref - - -from ._proj4 cimport (pj_init_plus, pj_free, pj_transform, pj_is_latlong, - pj_strerrno, pj_get_errno_ref, pj_get_release, - DEG_TO_RAD, RAD_TO_DEG) - - -cdef double NAN = float('nan') - - -PROJ4_RELEASE = pj_get_release().decode('utf-8') -_match = re.search(r"\d+\.\d+.\d+", PROJ4_RELEASE) -if _match is not None: - PROJ4_VERSION = tuple(int(v) for v in _match.group().split('.')) -else: - PROJ4_VERSION = () - -WGS84_SEMIMAJOR_AXIS = 6378137.0 -WGS84_SEMIMINOR_AXIS = 6356752.3142 - - -class Proj4Error(Exception): - """ - Raised when there has been an exception calling proj.4. - - Add a ``status`` attribute to the exception which has the - proj.4 error reference. - - """ - def __init__(self): - cdef int status - status = deref(pj_get_errno_ref()) - msg = 'Error from proj: {}'.format(pj_strerrno(status)) - self.status = status - Exception.__init__(self, msg) - - -def _safe_pj_transform_611(CRS src_crs not None, CRS tgt_crs not None, - int npts, int offset, - np.ndarray[np.double_t] x not None, - np.ndarray[np.double_t] y not None, - np.ndarray[np.double_t] z): - """ - Workaround bug in Proj 6.1.1+ with +to_meter on +proj=ob_tran. - - See https://github.com/OSGeo/proj#1782. - """ - cdef int status - - lonlat = ('latlon', 'latlong', 'lonlat', 'longlat') - - if (src_crs.proj4_params.get('proj', '') == 'ob_tran' and - src_crs.proj4_params.get('o_proj', '') in lonlat and - 'to_meter' in src_crs.proj4_params): - x *= src_crs.proj4_params['to_meter'] - y *= src_crs.proj4_params['to_meter'] - - if z is not None: - status = pj_transform(src_crs.proj4, tgt_crs.proj4, npts, offset, - &x[0], &y[0], &z[0]) - else: - status = pj_transform(src_crs.proj4, tgt_crs.proj4, npts, offset, - &x[0], &y[0], NULL) - - if (tgt_crs.proj4_params.get('proj', '') == 'ob_tran' and - tgt_crs.proj4_params.get('o_proj', '') in lonlat and - 'to_meter' in tgt_crs.proj4_params): - x /= tgt_crs.proj4_params['to_meter'] - y /= tgt_crs.proj4_params['to_meter'] - - return status - - -def _safe_pj_transform_pre_611(CRS src_crs not None, CRS tgt_crs not None, - int npts, int offset, - np.ndarray[np.double_t] x not None, - np.ndarray[np.double_t] y not None, - np.ndarray[np.double_t] z): - if z is not None: - return pj_transform(src_crs.proj4, tgt_crs.proj4, npts, offset, - &x[0], &y[0], &z[0]) - else: - return pj_transform(src_crs.proj4, tgt_crs.proj4, npts, offset, - &x[0], &y[0], NULL) - - -if (6, 1, 1) <= PROJ4_VERSION < (6, 3, 0): - _safe_pj_transform = _safe_pj_transform_611 -else: - _safe_pj_transform = _safe_pj_transform_pre_611 - - -class Globe(object): - """ - Define an ellipsoid and, optionally, how to relate it to the real world. - - """ - def __init__(self, datum=None, ellipse='WGS84', - semimajor_axis=None, semiminor_axis=None, - flattening=None, inverse_flattening=None, - towgs84=None, nadgrids=None): - """ - Parameters - ---------- - datum - Proj "datum" definition. Defaults to None. - ellipse - Proj "ellps" definition. Defaults to 'WGS84'. - semimajor_axis - Semimajor axis of the spheroid / ellipsoid. Defaults to None. - semiminor_axis - Semiminor axis of the ellipsoid. Defaults to None. - flattening - Flattening of the ellipsoid. Defaults to None. - inverse_flattening - Inverse flattening of the ellipsoid. Defaults to None. - towgs84 - Passed through to the Proj definition. Defaults to None. - nadgrids - Passed through to the Proj definition. Defaults to None. - - """ - self.datum = datum - self.ellipse = ellipse - self.semimajor_axis = semimajor_axis - self.semiminor_axis = semiminor_axis - self.flattening = flattening - self.inverse_flattening = inverse_flattening - self.towgs84 = towgs84 - self.nadgrids = nadgrids - - def to_proj4_params(self): - """ - Create an OrderedDict of key value pairs which represents this globe - in terms of proj params. - - """ - proj4_params = (['datum', self.datum], ['ellps', self.ellipse], - ['a', self.semimajor_axis], ['b', self.semiminor_axis], - ['f', self.flattening], ['rf', self.inverse_flattening], - ['towgs84', self.towgs84], ['nadgrids', self.nadgrids]) - return OrderedDict((k, v) for k, v in proj4_params if v is not None) - - -cdef class CRS: - """ - Define a Coordinate Reference System using proj. - - """ - - #: Whether this projection can handle ellipses. - _handles_ellipses = True - - def __cinit__(self): - self.proj4 = NULL - - def __dealloc__(self): - if self.proj4 != NULL: - pj_free(self.proj4) - - def __init__(self, proj4_params, globe=None): - """ - Parameters - ---------- - proj4_params: iterable of key-value pairs - The proj4 parameters required to define the - desired CRS. The parameters should not describe - the desired elliptic model, instead create an - appropriate Globe instance. The ``proj4_params`` - parameters will override any parameters that the - Globe defines. - globe: :class:`~cartopy.crs.Globe` instance, optional - If omitted, the default Globe instance will be created. - See :class:`~cartopy.crs.Globe` for details. - - """ - if globe is None: - if self._handles_ellipses: - globe = Globe() - else: - globe = Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS, - ellipse=None) - if not self._handles_ellipses: - a = globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS - b = globe.semiminor_axis or a - if a != b or globe.ellipse is not None: - warnings.warn('The "{}" projection does not handle elliptical ' - 'globes.'.format(self.__class__.__name__)) - self.globe = globe - self.proj4_params = self.globe.to_proj4_params() - self.proj4_params.update(proj4_params) - - init_items = [] - for k, v in self.proj4_params.items(): - if v is not None: - if isinstance(v, float): - init_items.append('+{}={:.16}'.format(k, v)) - elif isinstance(v, np.float32): - init_items.append('+{}={:.8}'.format(k, v)) - else: - init_items.append('+{}={}'.format(k, v)) - else: - init_items.append('+{}'.format(k)) - self.proj4_init = ' '.join(init_items) + ' +no_defs' - proj4_init_bytes = self.proj4_init.encode('utf-8') - if self.proj4 != NULL: - pj_free(self.proj4) - self.proj4 = pj_init_plus(proj4_init_bytes) - if not self.proj4: - raise Proj4Error() - - # Cython uses this method instead of the normal rich comparisons. - def __richcmp__(self, other, op): - # We're only interested in: - # == -> 2 - # != -> 3 - result = NotImplemented - if isinstance(other, CRS): - if op == 2: - result = self.proj4_init == other.proj4_init - elif op == 3: - result = self.proj4_init != other.proj4_init - return result - - def __hash__(self): - """Hash the CRS based on its proj4_init string.""" - return hash(self.proj4_init) - - def __reduce__(self): - """ - Implement the __reduce__ API so that unpickling produces a stateless - instance of this class (e.g. an empty tuple). The state will then be - added via __getstate__ and __setstate__. - - We are forced to this approach because a CRS does not store - the constructor keyword arguments in its state. - - """ - return self.__class__, (), self.__getstate__() - - def __getstate__(self): - """Return the full state of this instance for reconstruction - in ``__setstate__``. - - """ - state = self.__dict__.copy() - # Remove the proj4 instance and the proj4_init string, which can - # be re-created (in __setstate__) from the other arguments. - state.pop('proj4', None) - state.pop('proj4_init', None) - state['proj4_params'] = self.proj4_params - return state - - def __setstate__(self, state): - """ - Take the dictionary created by ``__getstate__`` and passes it - through to this implementation's __init__ method. - - """ - # Strip out the key state items for a CRS instance. - CRS_state = {key: state.pop(key) for key in ['proj4_params', 'globe']} - # Put everything else directly into the dict of the instance. - self.__dict__.update(state) - # Call the init of this class to ensure that the projection is - # properly initialised with proj4. - CRS.__init__(self, **CRS_state) - - # TODO - #def __str__ - #def _geod(self): # to return the pyproj.Geod - - def _as_mpl_transform(self, axes=None): - """ - Cast this CRS instance into a :class:`matplotlib.axes.Axes` using - the Matplotlib ``_as_mpl_transform`` interface. - - """ - # lazy import mpl.geoaxes (and therefore matplotlib) as mpl - # is only an optional dependency - import cartopy.mpl.geoaxes as geoaxes - if not isinstance(axes, geoaxes.GeoAxes): - raise ValueError('Axes should be an instance of GeoAxes, got %s' % type(axes)) - return geoaxes.InterProjectionTransform(self, axes.projection) + axes.transData - - property proj4_params: - def __get__(self): - return dict(self.proj4_params) - - def as_geocentric(self): - """ - Return a new Geocentric CRS with the same ellipse/datum as this - CRS. - - """ - return Geocentric(self.globe) - - def as_geodetic(self): - """ - Return a new Geodetic CRS with the same ellipse/datum as this - CRS. - - """ - return Geodetic(self.globe) - - cpdef is_geodetic(self): - return bool(pj_is_latlong(self.proj4)) - - def transform_point(self, double x, double y, CRS src_crs not None, trap=True): - """ - transform_point(x, y, src_crs) - - Transform the given float64 coordinate pair, in the given source - coordinate system (``src_crs``), to this coordinate system. - - Parameters - ---------- - x - the x coordinate, in ``src_crs`` coordinates, to transform - y - the y coordinate, in ``src_crs`` coordinates, to transform - src_crs - instance of :class:`CRS` that represents the coordinate - system of ``x`` and ``y``. - trap - Whether proj errors for "latitude or longitude exceeded limits" and - "tolerance condition error" should be trapped. - - Returns - ------- - (x, y) in this coordinate system - - """ - cdef: - np.ndarray[np.double_t, ndim=2] result - - result = self.transform_points(src_crs, np.array([x]), np.array([y])) - return result[0, 0], result[0, 1] - - - def transform_points(self, CRS src_crs not None, - np.ndarray x not None, - np.ndarray y not None, - np.ndarray z=None, trap=True): - """ - transform_points(src_crs, x, y[, z]) - - Transform the given coordinates, in the given source - coordinate system (``src_crs``), to this coordinate system. - - Parameters - ---------- - src_crs - instance of :class:`CRS` that represents the - coordinate system of ``x``, ``y`` and ``z``. - x - the x coordinates (array), in ``src_crs`` coordinates, - to transform. May be 1 or 2 dimensional. - y - the y coordinates (array), in ``src_crs`` coordinates, - to transform. Its shape must match that of x. - z: optional - the z coordinates (array), in ``src_crs`` coordinates, to - transform. Defaults to None. - If supplied, its shape must match that of x. - trap - Whether proj errors for "latitude or longitude exceeded limits" and - "tolerance condition error" should be trapped. - - Returns - ------- - Array of shape ``x.shape + (3, )`` in this coordinate system. - - """ - cdef: - np.ndarray[np.double_t, ndim=2] result - int status - - result_shape = tuple(x.shape[i] for i in range(x.ndim)) + (3, ) - - if z is None: - if x.ndim > 2 or y.ndim > 2: - raise ValueError('x and y arrays must be 1 or 2 dimensional') - elif x.ndim != 1 or y.ndim != 1: - x, y = x.flatten(), y.flatten() - - if x.shape[0] != y.shape[0]: - raise ValueError('x and y arrays must have the same length') - else: - if x.ndim > 2 or y.ndim > 2 or z.ndim > 2: - raise ValueError('x, y and z arrays must be 1 or 2 ' - 'dimensional') - elif x.ndim != 1 or y.ndim != 1 or z.ndim != 1: - x, y, z = x.flatten(), y.flatten(), z.flatten() - - if not x.shape[0] == y.shape[0] == z.shape[0]: - raise ValueError('x, y, and z arrays must have the same ' - 'length') - - npts = x.shape[0] - - result = np.empty([npts, 3], dtype=np.double) - if src_crs.is_geodetic(): - result[:, 0] = np.deg2rad(x) - result[:, 1] = np.deg2rad(y) - else: - result[:, 0] = x - result[:, 1] = y - # if a z has been given, put it in the result array which will be - # transformed in-place - if z is None: - result[:, 2] = 0 - else: - result[:, 2] = z - - # call proj. The result array is modified in place. This is only - # safe if npts is not 0. - if npts == 0: - return result - status = _safe_pj_transform(src_crs, self, npts, 3, - result[:, 0], result[:, 1], - result[:, 2]) - - if trap and status in (-14, -20): - # -14 => "latitude or longitude exceeded limits" - # -20 => "tolerance condition error" - result[:] = np.nan - elif trap and status != 0: - raise Proj4Error() - - if self.is_geodetic(): - result[:, :2] = np.rad2deg(result[:, :2]) - - if len(result_shape) > 2: - return result.reshape(result_shape) - - return result - - def transform_vectors(self, src_proj, x, y, u, v): - """ - transform_vectors(src_proj, x, y, u, v) - - Transform the given vector components, with coordinates in the - given source coordinate system (``src_proj``), to this coordinate - system. The vector components must be given relative to the - source projection's coordinate reference system (grid eastward and - grid northward). - - Parameters - ---------- - src_proj - The :class:`CRS.Projection` that represents the coordinate system - the vectors are defined in. - x - The x coordinates of the vectors in the source projection. - y - The y coordinates of the vectors in the source projection. - u - The grid-eastward components of the vectors. - v - The grid-northward components of the vectors. - - Note - ---- - x, y, u and v may be 1 or 2 dimensional, but must all have matching - shapes. - - Returns - ------- - ut, vt: The transformed vector components. - - Note - ---- - The algorithm used to transform vectors is an approximation - rather than an exact transform, but the accuracy should be - good enough for visualization purposes. - - """ - if not (x.shape == y.shape == u.shape == v.shape): - raise ValueError('x, y, u and v arrays must be the same shape') - if x.ndim not in (1, 2): - raise ValueError('x, y, u and v must be 1 or 2 dimensional') - # Transform the coordinates to the target projection. - proj_xyz = self.transform_points(src_proj, x, y) - target_x, target_y = proj_xyz[..., 0], proj_xyz[..., 1] - # Rotate the input vectors to the projection. - # - # 1: Find the magnitude and direction of the input vectors. - vector_magnitudes = (u**2 + v**2)**0.5 - vector_angles = np.arctan2(v, u) - # 2: Find a point in the direction of the original vector that is - # a small distance away from the base point of the vector (near - # the poles the point may have to be in the opposite direction - # to be valid). - factor = 360000. - delta = (src_proj.x_limits[1] - src_proj.x_limits[0]) / factor - x_perturbations = delta * np.cos(vector_angles) - y_perturbations = delta * np.sin(vector_angles) - # 3: Handle points that are invalid. These come from picking a new - # point that is outside the domain of the CRS. The first step is - # to apply the native transform to the input coordinates to make - # sure they are in the appropriate range. Then detect all the - # coordinates where the perturbation takes the point out of the - # valid x-domain and fix them. After that do the same for points - # that are outside the valid y-domain, which may reintroduce some - # points outside of the valid x-domain - proj_xyz = src_proj.transform_points(src_proj, x, y) - source_x, source_y = proj_xyz[..., 0], proj_xyz[..., 1] - # Detect all the coordinates where the perturbation takes the point - # outside of the valid x-domain, and reverse the direction of the - # perturbation to fix this. - eps = 1e-9 - invalid_x = np.logical_or( - source_x + x_perturbations < src_proj.x_limits[0]-eps, - source_x + x_perturbations > src_proj.x_limits[1]+eps) - if invalid_x.any(): - x_perturbations[invalid_x] *= -1 - y_perturbations[invalid_x] *= -1 - # Do the same for coordinates where the perturbation takes the point - # outside of the valid y-domain. This may reintroduce some points - # that will be outside the x-domain when the perturbation is - # applied. - invalid_y = np.logical_or( - source_y + y_perturbations < src_proj.y_limits[0]-eps, - source_y + y_perturbations > src_proj.y_limits[1]+eps) - if invalid_y.any(): - x_perturbations[invalid_y] *= -1 - y_perturbations[invalid_y] *= -1 - # Keep track of the points where the perturbation direction was - # reversed. - reversed_vectors = np.logical_xor(invalid_x, invalid_y) - # See if there were any points where we cannot reverse the direction - # of the perturbation to get the perturbed point within the valid - # domain of the projection, and issue a warning if there are. - problem_points = np.logical_or( - source_x + x_perturbations < src_proj.x_limits[0]-eps, - source_x + x_perturbations > src_proj.x_limits[1]+eps) - if problem_points.any(): - warnings.warn('Some vectors at source domain corners ' - 'may not have been transformed correctly') - # 4: Transform this set of points to the projection coordinates and - # find the angle between the base point and the perturbed point - # in the projection coordinates (reversing the direction at any - # points where the original was reversed in step 3). - proj_xyz = self.transform_points(src_proj, - source_x + x_perturbations, - source_y + y_perturbations) - target_x_perturbed = proj_xyz[..., 0] - target_y_perturbed = proj_xyz[..., 1] - projected_angles = np.arctan2(target_y_perturbed - target_y, - target_x_perturbed - target_x) - if reversed_vectors.any(): - projected_angles[reversed_vectors] += np.pi - # 5: Form the projected vector components, preserving the magnitude - # of the original vectors. - projected_u = vector_magnitudes * np.cos(projected_angles) - projected_v = vector_magnitudes * np.sin(projected_angles) - return projected_u, projected_v - - -class Geodetic(CRS): - """ - Define a latitude/longitude coordinate system with spherical topology, - geographical distance and coordinates are measured in degrees. - - """ - def __init__(self, globe=None): - """ - Parameters - ---------- - globe: A :class:`cartopy.crs.Globe`, optional - Defaults to a "WGS84" datum. - - """ - proj4_params = [('proj', 'lonlat')] - globe = globe or Globe(datum='WGS84') - super(Geodetic, self).__init__(proj4_params, globe) - - # XXX Implement fwd such as Basemap's Geod. Would be used in the tissot example. - # Could come from https://geographiclib.sourceforge.io - - -class Geocentric(CRS): - """ - Define a Geocentric coordinate system, where x, y, z are Cartesian - coordinates from the center of the Earth. - - """ - def __init__(self, globe=None): - """ - Parameters - ---------- - globe: A :class:`cartopy.crs.Globe`, optional - Defaults to a "WGS84" datum. - - """ - proj4_params = [('proj', 'geocent')] - globe = globe or Globe(datum='WGS84') - super(Geocentric, self).__init__(proj4_params, globe) diff --git a/lib/cartopy/_epsg.py b/lib/cartopy/_epsg.py index acbef771b..4e664025c 100644 --- a/lib/cartopy/_epsg.py +++ b/lib/cartopy/_epsg.py @@ -7,82 +7,20 @@ Provide support for converting EPSG codes to Projection instances. """ - -import numpy as np -import shapely.geometry as sgeom - import cartopy.crs as ccrs - - -_GLOBE_PARAMS = {'datum': 'datum', - 'ellps': 'ellipse', - 'a': 'semimajor_axis', - 'b': 'semiminor_axis', - 'f': 'flattening', - 'rf': 'inverse_flattening', - 'towgs84': 'towgs84', - 'nadgrids': 'nadgrids'} +from pyproj.crs import CRS as _CRS class _EPSGProjection(ccrs.Projection): def __init__(self, code): - import pyepsg - projection = pyepsg.get(code) - if not (isinstance(projection, pyepsg.ProjectedCRS) or - isinstance(projection, pyepsg.CompoundCRS)): + crs = _CRS.from_epsg(code) + if not crs.is_projected: raise ValueError('EPSG code does not define a projection') + if not crs.area_of_use: + raise ValueError("Area of use not defined.") self.epsg_code = code - - proj4_str = projection.as_proj4().strip() - terms = [term.strip('+').split('=') for term in proj4_str.split(' ')] - globe_terms = filter(lambda term: term[0] in _GLOBE_PARAMS, terms) - globe = ccrs.Globe(**{_GLOBE_PARAMS[name]: value for name, value in - globe_terms}) - other_terms = [] - for term in terms: - if term[0] not in _GLOBE_PARAMS: - if len(term) == 1: - other_terms.append([term[0], None]) - else: - other_terms.append(term) - super().__init__(other_terms, globe) - - # Convert lat/lon bounds to projected bounds. - # GML defines gmd:EX_GeographicBoundingBox as: - # Geographic area of the entire dataset referenced to WGS 84 - # NB. We can't use a polygon transform at this stage because - # that relies on the existence of the map boundary... the very - # thing we're trying to work out! ;-) - x0, x1, y0, y1 = projection.domain_of_validity() - geodetic = ccrs.Geodetic() - lons = np.array([x0, x0, x1, x1]) - lats = np.array([y0, y1, y1, y0]) - points = self.transform_points(geodetic, lons, lats) - x = points[:, 0] - y = points[:, 1] - self.bounds = (x.min(), x.max(), y.min(), y.max()) + super().__init__(crs.to_wkt()) def __repr__(self): return '_EPSGProjection({})'.format(self.epsg_code) - - @property - def boundary(self): - x0, x1, y0, y1 = self.bounds - return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0), - (x0, y0)]) - - @property - def x_limits(self): - x0, x1, y0, y1 = self.bounds - return (x0, x1) - - @property - def y_limits(self): - x0, x1, y0, y1 = self.bounds - return (y0, y1) - - @property - def threshold(self): - x0, x1, y0, y1 = self.bounds - return min(x1 - x0, y1 - y0) / 100. diff --git a/lib/cartopy/_proj4.pxd b/lib/cartopy/_proj4.pxd deleted file mode 100644 index 336d15edf..000000000 --- a/lib/cartopy/_proj4.pxd +++ /dev/null @@ -1,28 +0,0 @@ -# Copyright Cartopy Contributors -# -# This file is part of Cartopy and is released under the LGPL license. -# See COPYING and COPYING.LESSER in the root of the repository for full -# licensing details. - -""" -This file declares the Proj API, version 4. - -""" - - -cdef extern from "proj_api.h": - ctypedef void *projPJ - ctypedef struct projLP: - double u - double v - - projPJ pj_init_plus(char *) nogil - void pj_free(projPJ) nogil - void pj_get_spheroid_defn(projPJ, double *, double *) nogil - int pj_transform(projPJ, projPJ, long, int, double *, double *, double *) nogil - int pj_is_latlong(projPJ) nogil - char *pj_strerrno(int) nogil - int *pj_get_errno_ref() nogil - char *pj_get_release() nogil - cdef double DEG_TO_RAD - cdef double RAD_TO_DEG diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 36f0181e8..ede0f4107 100644 --- a/lib/cartopy/crs.py +++ b/lib/cartopy/crs.py @@ -10,23 +10,619 @@ """ -from abc import ABCMeta, abstractproperty +from abc import ABCMeta +from collections import OrderedDict import io +import json import math import warnings import numpy as np import shapely.geometry as sgeom +from pyproj import Transformer +from pyproj.exceptions import ProjError from shapely.prepared import prep -from cartopy._crs import (CRS, Geodetic, Globe, PROJ4_VERSION, - WGS84_SEMIMAJOR_AXIS, WGS84_SEMIMINOR_AXIS) -from cartopy._crs import Geocentric # noqa: F401 (flake8 = unused import) import cartopy.trace +try: + # https://github.com/pyproj4/pyproj/pull/912 + from pyproj.crs import CustomConstructorCRS as _CRS +except ImportError: + from pyproj import CRS as _CRS + __document_these__ = ['CRS', 'Geocentric', 'Geodetic', 'Globe'] +PROJ_VERSION = cartopy.trace.PROJ_VERSION +WGS84_SEMIMAJOR_AXIS = 6378137.0 +WGS84_SEMIMINOR_AXIS = 6356752.3142 + + +def _safe_pj_transform_611(src_crs, tgt_crs, x, y, z=None, trap=True): + """ + Workaround bug in Proj 6.1.1+ with +to_meter on +proj=ob_tran. + + See https://github.com/OSGeo/proj#1782. + """ + lonlat = ('latlon', 'latlong', 'lonlat', 'longlat') + + if ( + src_crs.proj4_params.get('proj', '') == 'ob_tran' and + src_crs.proj4_params.get('o_proj', '') in lonlat and + 'to_meter' in src_crs.proj4_params + ): + x *= src_crs.proj4_params['to_meter'] + y *= src_crs.proj4_params['to_meter'] + + transformer = Transformer.from_crs(src_crs, tgt_crs, always_xy=True) + transformed_coords = transformer.transform(x, y, z, errcheck=trap) + if z is None: + xx, yy = transformed_coords + zz = 0 + else: + xx, yy, zz = transformed_coords + + if ( + tgt_crs.proj4_params.get('proj', '') == 'ob_tran' and + tgt_crs.proj4_params.get('o_proj', '') in lonlat and + 'to_meter' in tgt_crs.proj4_params + ): + xx /= tgt_crs.proj4_params['to_meter'] + yy /= tgt_crs.proj4_params['to_meter'] + + return xx, yy, zz + + +def _safe_pj_transform_pre_611(src_crs, tgt_crs, x, y, z=None, trap=True): + transformer = Transformer.from_crs(src_crs, tgt_crs, always_xy=True) + transformed_coords = transformer.transform(x, y, z, errcheck=trap) + if z is None: + xx, yy = transformed_coords + zz = 0 + else: + xx, yy, zz = transformed_coords + return xx, yy, zz + + +if (6, 1, 1) <= PROJ_VERSION < (6, 3, 0): + _safe_pj_transform = _safe_pj_transform_611 +else: + _safe_pj_transform = _safe_pj_transform_pre_611 + + +class Globe(object): + """ + Define an ellipsoid and, optionally, how to relate it to the real world. + + """ + def __init__(self, datum=None, ellipse='WGS84', + semimajor_axis=None, semiminor_axis=None, + flattening=None, inverse_flattening=None, + towgs84=None, nadgrids=None): + """ + Parameters + ---------- + datum + Proj "datum" definition. Defaults to None. + ellipse + Proj "ellps" definition. Defaults to 'WGS84'. + semimajor_axis + Semimajor axis of the spheroid / ellipsoid. Defaults to None. + semiminor_axis + Semiminor axis of the ellipsoid. Defaults to None. + flattening + Flattening of the ellipsoid. Defaults to None. + inverse_flattening + Inverse flattening of the ellipsoid. Defaults to None. + towgs84 + Passed through to the Proj definition. Defaults to None. + nadgrids + Passed through to the Proj definition. Defaults to None. + + """ + self.datum = datum + self.ellipse = ellipse + self.semimajor_axis = semimajor_axis + self.semiminor_axis = semiminor_axis + self.flattening = flattening + self.inverse_flattening = inverse_flattening + self.towgs84 = towgs84 + self.nadgrids = nadgrids + + def to_proj4_params(self): + """ + Create an OrderedDict of key value pairs which represents this globe + in terms of proj params. + + """ + proj4_params = ( + ['datum', self.datum], + ['ellps', self.ellipse], + ['a', self.semimajor_axis], + ['b', self.semiminor_axis], + ['f', self.flattening], + ['rf', self.inverse_flattening], + ['towgs84', self.towgs84], + ['nadgrids', self.nadgrids] + ) + return OrderedDict((k, v) for k, v in proj4_params if v is not None) + + +class CRS(_CRS): + """ + Define a Coordinate Reference System using proj. + """ + + #: Whether this projection can handle ellipses. + _handles_ellipses = True + + def __init__(self, proj4_params, globe=None): + """ + Parameters + ---------- + proj4_params: iterable of key-value pairs + The proj4 parameters required to define the + desired CRS. The parameters should not describe + the desired elliptic model, instead create an + appropriate Globe instance. The ``proj4_params`` + parameters will override any parameters that the + Globe defines. + globe: :class:`~cartopy.crs.Globe` instance, optional + If omitted, the default Globe instance will be created. + See :class:`~cartopy.crs.Globe` for details. + + """ + # for compatibility with pyproj.CRS and rasterio.crs.CRS + try: + proj4_params = proj4_params.to_wkt() + except AttributeError: + pass + # handle PROJ JSON + if ( + isinstance(proj4_params, dict) and + "proj" not in proj4_params and + "init" not in proj4_params + ): + proj4_params = json.dumps(proj4_params) + + if globe is not None and isinstance(proj4_params, str): + raise ValueError("Cannot have 'globe' with string params.") + if globe is None and not isinstance(proj4_params, str): + if self._handles_ellipses: + globe = Globe() + else: + globe = Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS, + ellipse=None) + if globe is not None and not self._handles_ellipses: + a = globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS + b = globe.semiminor_axis or a + if a != b or globe.ellipse is not None: + warnings.warn('The "{}" projection does not handle elliptical ' + 'globes.'.format(self.__class__.__name__)) + self.globe = globe + if isinstance(proj4_params, str): + self._proj4_params = {} + self.proj4_init = proj4_params + else: + self._proj4_params = self.globe.to_proj4_params() + self._proj4_params.update(proj4_params) + + init_items = [] + for k, v in self._proj4_params.items(): + if v is not None: + if isinstance(v, float): + init_items.append('+{}={:.16}'.format(k, v)) + elif isinstance(v, np.float32): + init_items.append('+{}={:.8}'.format(k, v)) + else: + init_items.append('+{}={}'.format(k, v)) + else: + init_items.append('+{}'.format(k)) + self.proj4_init = ' '.join(init_items) + ' +no_defs' + super().__init__(self.proj4_init) + + def __eq__(self, other): + if isinstance(other, CRS) and self.proj4_init == other.proj4_init: + # Fast path Cartopy's CRS + return True + # For everything else, we let pyproj handle the comparison + return super().__eq__(other) + + def __hash__(self): + """Hash the CRS based on its proj4_init string.""" + return hash(self.proj4_init) + + def __reduce__(self): + """ + Implement the __reduce__ API so that unpickling produces a stateless + instance of this class (e.g. an empty tuple). The state will then be + added via __getstate__ and __setstate__. + We are forced to this approach because a CRS does not store + the constructor keyword arguments in its state. + """ + return self.__class__, (), self.__getstate__() + + def __getstate__(self): + """Return the full state of this instance for reconstruction + in ``__setstate__``. + """ + state = self.__dict__.copy() + # remove pyproj specific attrs + state.pop('srs', None) + state.pop('_local', None) + # Remove the proj4 instance and the proj4_init string, which can + # be re-created (in __setstate__) from the other arguments. + state.pop('proj4', None) + state.pop('proj4_init', None) + state['proj4_params'] = self.proj4_params + return state + + def __setstate__(self, state): + """ + Take the dictionary created by ``__getstate__`` and passes it + through to this implementation's __init__ method. + """ + # Strip out the key state items for a CRS instance. + CRS_state = {key: state.pop(key) for key in ['proj4_params', 'globe']} + # Put everything else directly into the dict of the instance. + self.__dict__.update(state) + # Call the init of this class to ensure that the projection is + # properly initialised with proj4. + CRS.__init__(self, **CRS_state) + + def _as_mpl_transform(self, axes=None): + """ + Cast this CRS instance into a :class:`matplotlib.axes.Axes` using + the Matplotlib ``_as_mpl_transform`` interface. + + """ + # lazy import mpl.geoaxes (and therefore matplotlib) as mpl + # is only an optional dependency + import cartopy.mpl.geoaxes as geoaxes + if not isinstance(axes, geoaxes.GeoAxes): + raise ValueError( + 'Axes should be an instance of GeoAxes, got %s' % type(axes) + ) + return ( + geoaxes.InterProjectionTransform(self, axes.projection) + + axes.transData + ) + + @property + def proj4_params(self): + return dict(self._proj4_params) + + def as_geocentric(self): + """ + Return a new Geocentric CRS with the same ellipse/datum as this + CRS. + + """ + return CRS( + { + "$schema": ( + "https://proj.org/schemas/v0.2/projjson.schema.json" + ), + "type": "GeodeticCRS", + "name": "unknown", + "datum": self.datum.to_json_dict(), + "coordinate_system": { + "subtype": "Cartesian", + "axis": [ + { + "name": "Geocentric X", + "abbreviation": "X", + "direction": "geocentricX", + "unit": "metre" + }, + { + "name": "Geocentric Y", + "abbreviation": "Y", + "direction": "geocentricY", + "unit": "metre" + }, + { + "name": "Geocentric Z", + "abbreviation": "Z", + "direction": "geocentricZ", + "unit": "metre" + } + ] + } + } + ) + + def as_geodetic(self): + """ + Return a new Geodetic CRS with the same ellipse/datum as this + CRS. + + """ + return CRS(self.geodetic_crs.srs) + + def is_geodetic(self): + return self.is_geographic + + def transform_point(self, x, y, src_crs, trap=True): + """ + transform_point(x, y, src_crs) + + Transform the given float64 coordinate pair, in the given source + coordinate system (``src_crs``), to this coordinate system. + + Parameters + ---------- + x + the x coordinate, in ``src_crs`` coordinates, to transform + y + the y coordinate, in ``src_crs`` coordinates, to transform + src_crs + instance of :class:`CRS` that represents the coordinate + system of ``x`` and ``y``. + trap + Whether proj errors for "latitude or longitude exceeded limits" and + "tolerance condition error" should be trapped. + + Returns + ------- + (x, y) in this coordinate system + + """ + result = self.transform_points( + src_crs, np.array([x]), np.array([y]), trap=trap, + ).reshape((1, 3)) + return result[0, 0], result[0, 1] + + def transform_points(self, src_crs, x, y, z=None, trap=False): + """ + transform_points(src_crs, x, y[, z]) + + Transform the given coordinates, in the given source + coordinate system (``src_crs``), to this coordinate system. + + Parameters + ---------- + src_crs + instance of :class:`CRS` that represents the + coordinate system of ``x``, ``y`` and ``z``. + x + the x coordinates (array), in ``src_crs`` coordinates, + to transform. May be 1 or 2 dimensional. + y + the y coordinates (array), in ``src_crs`` coordinates, + to transform. Its shape must match that of x. + z: optional + the z coordinates (array), in ``src_crs`` coordinates, to + transform. Defaults to None. + If supplied, its shape must match that of x. + trap + Whether proj errors for "latitude or longitude exceeded limits" and + "tolerance condition error" should be trapped. + + Returns + ------- + Array of shape ``x.shape + (3, )`` in this coordinate system. + + """ + result_shape = tuple(x.shape[i] for i in range(x.ndim)) + (3, ) + + if z is None: + if x.ndim > 2 or y.ndim > 2: + raise ValueError('x and y arrays must be 1 or 2 dimensional') + elif x.ndim != 1 or y.ndim != 1: + x, y = x.flatten(), y.flatten() + + if x.shape[0] != y.shape[0]: + raise ValueError('x and y arrays must have the same length') + else: + if x.ndim > 2 or y.ndim > 2 or z.ndim > 2: + raise ValueError('x, y and z arrays must be 1 or 2 ' + 'dimensional') + elif x.ndim != 1 or y.ndim != 1 or z.ndim != 1: + x, y, z = x.flatten(), y.flatten(), z.flatten() + + if not x.shape[0] == y.shape[0] == z.shape[0]: + raise ValueError('x, y, and z arrays must have the same ' + 'length') + + npts = x.shape[0] + + result = np.empty([npts, 3], dtype=np.double) + if npts: + if self == src_crs and ( + isinstance(src_crs, _CylindricalProjection) or + self.is_geodetic()): + # convert from [0,360] to [-180,180] + x = np.array(x, copy=True) + to_180 = x > 180 + x[to_180] = (((x[to_180] + 180) % 360) - 180) + try: + result[:, 0], result[:, 1], result[:, 2] = \ + _safe_pj_transform(src_crs, self, x, y, z, trap=trap) + 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 + ): + result[:] = np.nan + else: + raise + + if not trap: + result[np.isinf(result)] = np.nan + + if len(result_shape) > 2: + return result.reshape(result_shape) + + return result + + def transform_vectors(self, src_proj, x, y, u, v): + """ + transform_vectors(src_proj, x, y, u, v) + + Transform the given vector components, with coordinates in the + given source coordinate system (``src_proj``), to this coordinate + system. The vector components must be given relative to the + source projection's coordinate reference system (grid eastward and + grid northward). + + Parameters + ---------- + src_proj + The :class:`CRS.Projection` that represents the coordinate system + the vectors are defined in. + x + The x coordinates of the vectors in the source projection. + y + The y coordinates of the vectors in the source projection. + u + The grid-eastward components of the vectors. + v + The grid-northward components of the vectors. + + Note + ---- + x, y, u and v may be 1 or 2 dimensional, but must all have matching + shapes. + + Returns + ------- + ut, vt: The transformed vector components. + + Note + ---- + The algorithm used to transform vectors is an approximation + rather than an exact transform, but the accuracy should be + good enough for visualization purposes. + + """ + if not (x.shape == y.shape == u.shape == v.shape): + raise ValueError('x, y, u and v arrays must be the same shape') + if x.ndim not in (1, 2): + raise ValueError('x, y, u and v must be 1 or 2 dimensional') + # Transform the coordinates to the target projection. + proj_xyz = self.transform_points(src_proj, x, y) + target_x, target_y = proj_xyz[..., 0], proj_xyz[..., 1] + # Rotate the input vectors to the projection. + # + # 1: Find the magnitude and direction of the input vectors. + vector_magnitudes = (u**2 + v**2)**0.5 + vector_angles = np.arctan2(v, u) + # 2: Find a point in the direction of the original vector that is + # a small distance away from the base point of the vector (near + # the poles the point may have to be in the opposite direction + # to be valid). + factor = 360000. + delta = (src_proj.x_limits[1] - src_proj.x_limits[0]) / factor + x_perturbations = delta * np.cos(vector_angles) + y_perturbations = delta * np.sin(vector_angles) + # 3: Handle points that are invalid. These come from picking a new + # point that is outside the domain of the CRS. The first step is + # to apply the native transform to the input coordinates to make + # sure they are in the appropriate range. Then detect all the + # coordinates where the perturbation takes the point out of the + # valid x-domain and fix them. After that do the same for points + # that are outside the valid y-domain, which may reintroduce some + # points outside of the valid x-domain + proj_xyz = src_proj.transform_points(src_proj, x, y) + source_x, source_y = proj_xyz[..., 0], proj_xyz[..., 1] + # Detect all the coordinates where the perturbation takes the point + # outside of the valid x-domain, and reverse the direction of the + # perturbation to fix this. + eps = 1e-9 + invalid_x = np.logical_or( + source_x + x_perturbations < src_proj.x_limits[0]-eps, + source_x + x_perturbations > src_proj.x_limits[1]+eps) + if invalid_x.any(): + x_perturbations[invalid_x] *= -1 + y_perturbations[invalid_x] *= -1 + # Do the same for coordinates where the perturbation takes the point + # outside of the valid y-domain. This may reintroduce some points + # that will be outside the x-domain when the perturbation is + # applied. + invalid_y = np.logical_or( + source_y + y_perturbations < src_proj.y_limits[0]-eps, + source_y + y_perturbations > src_proj.y_limits[1]+eps) + if invalid_y.any(): + x_perturbations[invalid_y] *= -1 + y_perturbations[invalid_y] *= -1 + # Keep track of the points where the perturbation direction was + # reversed. + reversed_vectors = np.logical_xor(invalid_x, invalid_y) + # See if there were any points where we cannot reverse the direction + # of the perturbation to get the perturbed point within the valid + # domain of the projection, and issue a warning if there are. + problem_points = np.logical_or( + source_x + x_perturbations < src_proj.x_limits[0]-eps, + source_x + x_perturbations > src_proj.x_limits[1]+eps) + if problem_points.any(): + warnings.warn('Some vectors at source domain corners ' + 'may not have been transformed correctly') + # 4: Transform this set of points to the projection coordinates and + # find the angle between the base point and the perturbed point + # in the projection coordinates (reversing the direction at any + # points where the original was reversed in step 3). + proj_xyz = self.transform_points(src_proj, + source_x + x_perturbations, + source_y + y_perturbations) + target_x_perturbed = proj_xyz[..., 0] + target_y_perturbed = proj_xyz[..., 1] + projected_angles = np.arctan2(target_y_perturbed - target_y, + target_x_perturbed - target_x) + if reversed_vectors.any(): + projected_angles[reversed_vectors] += np.pi + # 5: Form the projected vector components, preserving the magnitude + # of the original vectors. + projected_u = vector_magnitudes * np.cos(projected_angles) + projected_v = vector_magnitudes * np.sin(projected_angles) + return projected_u, projected_v + + +class Geodetic(CRS): + """ + Define a latitude/longitude coordinate system with spherical topology, + geographical distance and coordinates are measured in degrees. + + """ + def __init__(self, globe=None): + """ + Parameters + ---------- + globe: A :class:`cartopy.crs.Globe`, optional + Defaults to a "WGS84" datum. + + """ + proj4_params = [('proj', 'lonlat')] + globe = globe or Globe(datum='WGS84') + super(Geodetic, self).__init__(proj4_params, globe) + + # XXX Implement fwd such as Basemap's Geod. + # Would be used in the tissot example. + # Could come from https://geographiclib.sourceforge.io + + +class Geocentric(CRS): + """ + Define a Geocentric coordinate system, where x, y, z are Cartesian + coordinates from the center of the Earth. + + """ + def __init__(self, globe=None): + """ + Parameters + ---------- + globe: A :class:`cartopy.crs.Globe`, optional + Defaults to a "WGS84" datum. + + """ + proj4_params = [('proj', 'geocent')] + globe = globe or Globe(datum='WGS84') + super(Geocentric, self).__init__(proj4_params, globe) + class RotatedGeodetic(CRS): """ @@ -59,12 +655,14 @@ def __init__(self, pole_longitude, pole_latitude, A :class:`cartopy.crs.Globe`. Defaults to a "WGS84" datum. """ + globe = globe or Globe(datum='WGS84') proj4_params = [('proj', 'ob_tran'), ('o_proj', 'latlon'), ('o_lon_p', central_rotated_longitude), ('o_lat_p', pole_latitude), ('lon_0', 180 + pole_longitude), - ('to_meter', math.radians(1))] - globe = globe or Globe(datum='WGS84') + ('to_meter', math.radians(1) * ( + globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS))] + super().__init__(proj4_params, globe=globe) @@ -85,17 +683,49 @@ class Projection(CRS, metaclass=ABCMeta): 'MultiPolygon': '_project_multipolygon', } - @abstractproperty + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.bounds = None + if self.area_of_use: + # Convert lat/lon bounds to projected bounds. + # Geographic area of the entire dataset referenced to WGS 84 + # NB. We can't use a polygon transform at this stage because + # that relies on the existence of the map boundary... the very + # thing we're trying to work out! ;-) + x0 = self.area_of_use.west + x1 = self.area_of_use.east + y0 = self.area_of_use.south + y1 = self.area_of_use.north + lons = np.array([x0, x0, x1, x1]) + lats = np.array([y0, y1, y1, y0]) + points = self.transform_points(self.as_geodetic(), lons, lats) + x = points[:, 0] + y = points[:, 1] + self.bounds = (x.min(), x.max(), y.min(), y.max()) + x0, x1, y0, y1 = self.bounds + self.threshold = min(x1 - x0, y1 - y0) / 100. + + @property def boundary(self): - pass + if self.bounds is None: + raise NotImplementedError + x0, x1, y0, y1 = self.bounds + return sgeom.LineString([(x0, y0), (x0, y1), (x1, y1), (x1, y0), + (x0, y0)]) - @abstractproperty + @property def x_limits(self): - pass + if self.bounds is None: + raise NotImplementedError + x0, x1, y0, y1 = self.bounds + return (x0, x1) - @abstractproperty + @property def y_limits(self): - pass + if self.bounds is None: + raise NotImplementedError + x0, x1, y0, y1 = self.bounds + return (y0, y1) @property def threshold(self): @@ -131,6 +761,9 @@ def domain(self): domain = self._domain = sgeom.Polygon(self.boundary) return domain + def is_geodetic(self): + return False + def _determine_longitude_bounds(self, central_longitude): # In new proj, using exact limits will wrap-around, so subtract a # small epsilon: @@ -166,7 +799,7 @@ def _repr_html_(self): # "Rewind" the buffer to the start and return it as an svg string. buf.seek(0) svg = buf.read() - return '{}
{}
'.format(svg, escape(repr(self))) + return '{}
{}
'.format(svg, escape(object.__repr__(self))) def _as_mpl_axes(self): import cartopy.mpl.geoaxes as geoaxes @@ -700,12 +1333,13 @@ def _ellipse_boundary(semimajor=2, semiminor=1, easting=0, northing=0, n=201): class PlateCarree(_CylindricalProjection): def __init__(self, central_longitude=0.0, globe=None): - proj4_params = [('proj', 'eqc'), ('lon_0', central_longitude)] - if globe is None: - globe = Globe(semimajor_axis=math.degrees(1)) - a_rad = math.radians(globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) - x_max = a_rad * 180 - y_max = a_rad * 90 + globe = globe or Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS) + proj4_params = [('proj', 'eqc'), ('lon_0', central_longitude), + ('to_meter', math.radians(1) * ( + globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS)), + ('vto_meter', 1)] + x_max = 180 + y_max = 90 # Set the threshold around 0.5 if the x max is 180. self.threshold = x_max / 360 super().__init__(proj4_params, x_max, y_max, globe=globe) @@ -720,7 +1354,7 @@ def _bbox_and_offset(self, other_plate_carree): >>> src = ccrs.PlateCarree(central_longitude=10) >>> bboxes, offset = ccrs.PlateCarree()._bbox_and_offset(src) >>> print(bboxes) - [[-180.0, -170.0], [-170.0, 180.0]] + [[-180, -170.0], [-170.0, 180]] >>> print(offset) 10.0 @@ -769,7 +1403,6 @@ def quick_vertices_transform(self, vertices, src_crs): potential = (self_params == src_params and self.y_limits[0] <= ys.min() and self.y_limits[1] >= ys.max()) - if potential: mod = np.diff(src_crs.x_limits)[0] bboxes, proj_offset = self._bbox_and_offset(src_crs) @@ -833,7 +1466,7 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, ('lat_0', central_latitude), ('k', scale_factor), ('x_0', false_easting), ('y_0', false_northing), ('units', 'm')] - if PROJ4_VERSION < (6, 0, 0): + if PROJ_VERSION < (6, 0, 0): if not approx: proj4_params[0] = ('proj', 'etmerc') else: @@ -1046,7 +1679,7 @@ def __init__(self, central_longitude=0.0, # Calculate limits. minlon, maxlon = self._determine_longitude_bounds(central_longitude) - limits = self.transform_points(Geodetic(), + limits = self.transform_points(self.as_geodetic(), np.array([minlon, maxlon]), np.array([min_latitude, max_latitude])) self._x_limits = tuple(limits[..., 0]) @@ -1096,9 +1729,11 @@ def y_limits(self): class LambertCylindrical(_RectangularProjection): - def __init__(self, central_longitude=0.0): - proj4_params = [('proj', 'cea'), ('lon_0', central_longitude)] - globe = Globe(semimajor_axis=math.degrees(1)) + def __init__(self, central_longitude=0.0, globe=None): + globe = globe or Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS) + proj4_params = [('proj', 'cea'), ('lon_0', central_longitude), + ('to_meter', math.radians(1) * ( + globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS))] super().__init__(proj4_params, 180, math.degrees(1), globe=globe) @@ -1266,14 +1901,14 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, super().__init__(proj4_params, globe=globe) - a = float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) + a = float(self.ellipsoid.semi_major_metre or WGS84_SEMIMAJOR_AXIS) # Find the antipode, and shift it a small amount in latitude to # approximate the extent of the projection: lon = central_longitude + 180 sign = np.sign(central_latitude) or 1 lat = -central_latitude + sign * 0.01 - x, max_y = self.transform_point(lon, lat, PlateCarree()) + x, max_y = self.transform_point(lon, lat, PlateCarree(globe=globe)) coords = _ellipse_boundary(a * 1.9999, max_y - false_northing, false_easting, false_northing, 61) @@ -1302,10 +1937,9 @@ class Miller(_RectangularProjection): def __init__(self, central_longitude=0.0, globe=None): if globe is None: - globe = Globe(semimajor_axis=math.degrees(1), ellipse=None) + globe = Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS, ellipse=None) - # TODO: Let the globe return the semimajor axis always. - a = float(globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) + a = globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS proj4_params = [('proj', 'mill'), ('lon_0', central_longitude)] # See Snyder, 1987. Eqs (11-1) and (11-2) substituting maximums of @@ -1347,12 +1981,13 @@ def __init__(self, pole_longitude=0.0, pole_latitude=90.0, datum. """ - + globe = globe or Globe(semimajor_axis=WGS84_SEMIMAJOR_AXIS) proj4_params = [('proj', 'ob_tran'), ('o_proj', 'latlon'), ('o_lon_p', central_rotated_longitude), ('o_lat_p', pole_latitude), ('lon_0', 180 + pole_longitude), - ('to_meter', math.radians(1))] + ('to_meter', math.radians(1) * ( + globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS))] super().__init__(proj4_params, 180, 90, globe=globe) @@ -1389,8 +2024,8 @@ def __init__(self, central_latitude=0.0, central_longitude=0.0, # incorrect transformation with lon_0=0 (see # https://github.com/OSGeo/proj.4/issues/194). if central_latitude == 0: - if PROJ4_VERSION != (): - if PROJ4_VERSION < (5, 0, 0): + if PROJ_VERSION != (): + if PROJ_VERSION < (5, 0, 0): warnings.warn( 'The Stereographic projection in Proj older than ' '5.0.0 incorrectly transforms points when ' @@ -1427,8 +2062,8 @@ def __init__(self, central_latitude=0.0, central_longitude=0.0, super().__init__(proj4_params, globe=globe) # TODO: Let the globe return the semimajor axis always. - a = float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) - b = float(self.globe.semiminor_axis or WGS84_SEMIMINOR_AXIS) + a = float(self.ellipsoid.semi_major_metre or WGS84_SEMIMAJOR_AXIS) + b = float(self.ellipsoid.semi_minor_metre or WGS84_SEMIMINOR_AXIS) # Note: The magic number has been picked to maintain consistent # behaviour with a wgs84 globe. There is no guarantee that the scaling @@ -1482,8 +2117,8 @@ class Orthographic(Projection): def __init__(self, central_longitude=0.0, central_latitude=0.0, globe=None): - if PROJ4_VERSION != (): - if (5, 0, 0) <= PROJ4_VERSION < (5, 1, 0): + if PROJ_VERSION != (): + if (5, 0, 0) <= PROJ_VERSION < (5, 1, 0): warnings.warn( 'The Orthographic projection in the v5.0.x series of Proj ' 'incorrectly transforms points. Use this projection with ' @@ -1500,7 +2135,7 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, super().__init__(proj4_params, globe=globe) # TODO: Let the globe return the semimajor axis always. - a = float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) + a = float(self.ellipsoid.semi_major_metre or WGS84_SEMIMAJOR_AXIS) # To stabilise the projection of geometries, we reduce the boundary by # a tiny fraction at the cost of the extreme edges. @@ -1717,10 +2352,10 @@ def __init__(self, central_longitude=0, false_easting=None, If omitted, a default globe is created. """ - if PROJ4_VERSION < (5, 2, 0): + if PROJ_VERSION < (5, 2, 0): raise ValueError('The EqualEarth projection requires Proj version ' '5.2.0, but you are using {}.' - .format('.'.join(str(v) for v in PROJ4_VERSION))) + .format('.'.join(str(v) for v in PROJ_VERSION))) proj_params = [('proj', 'eqearth'), ('lon_0', central_longitude)] super().__init__(proj_params, central_longitude, @@ -1806,8 +2441,8 @@ def __init__(self, central_longitude=0, globe=None, # Warn when using Robinson with proj 4.8 due to discontinuity at # 40 deg N introduced by incomplete fix to issue #113 (see # https://github.com/OSGeo/proj.4/issues/113). - if PROJ4_VERSION != (): - if (4, 8) <= PROJ4_VERSION < (4, 9): + if PROJ_VERSION != (): + if (4, 8) <= PROJ_VERSION < (4, 9): warnings.warn('The Robinson projection in the v4.8.x series ' 'of Proj contains a discontinuity at ' '40 deg latitude. Use this projection with ' @@ -1826,7 +2461,7 @@ def __init__(self, central_longitude=0, globe=None, globe=globe) self.threshold = 1e4 - def transform_point(self, x, y, src_crs): + def transform_point(self, x, y, src_crs, trap=True): """ Capture and handle any input NaNs, else invoke parent function, :meth:`_WarpedRectangularProjection.transform_point`. @@ -1843,10 +2478,10 @@ def transform_point(self, x, y, src_crs): if np.isnan(x) or np.isnan(y): result = (np.nan, np.nan) else: - result = super().transform_point(x, y, src_crs) + result = super().transform_point(x, y, src_crs, trap=trap) return result - def transform_points(self, src_crs, x, y, z=None): + def transform_points(self, src_crs, x, y, z=None, trap=False): """ Capture and handle NaNs in input points -- else as parent function, :meth:`_WarpedRectangularProjection.transform_points`. @@ -1871,7 +2506,7 @@ def transform_points(self, src_crs, x, y, z=None): y[input_point_nans] = 0.0 if z is not None: z[input_point_nans] = 0.0 - result = super().transform_points(src_crs, x, y, z) + result = super().transform_points(src_crs, x, y, z, trap=trap) if handle_nans: # Result always has shape (N, 3). # Blank out each (whole) point where we had a NaN in the input. @@ -1910,8 +2545,8 @@ def __init__(self, central_longitude=0, globe=None, emphasis='land'): super().__init__(proj4_params, globe=globe) elif emphasis == 'ocean': - if PROJ4_VERSION < (7, 1, 0): - _proj_ver = '.'.join(str(v) for v in PROJ4_VERSION) + if PROJ_VERSION < (7, 1, 0): + _proj_ver = '.'.join(str(v) for v in PROJ_VERSION) raise ValueError('The Interrupted Goode Homolosine ocean ' 'projection requires Proj version 7.1.0, ' 'but you are using ' + _proj_ver) @@ -2075,7 +2710,7 @@ def __init__(self, central_longitude=0.0, satellite_height=35785831, sweep_axis=sweep_axis) # TODO: Let the globe return the semimajor axis always. - a = float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) + a = float(self.ellipsoid.semi_major_metre or WGS84_SEMIMAJOR_AXIS) h = float(satellite_height) # These are only exact for a spherical Earth, owing to assuming a is @@ -2140,7 +2775,7 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, globe=globe) # TODO: Let the globe return the semimajor axis always. - a = self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS + a = self.ellipsoid.semi_major_metre or WGS84_SEMIMAJOR_AXIS h = float(satellite_height) max_x = a * np.sqrt(h / (2 * a + h)) @@ -2264,8 +2899,8 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, # Warn when using Azimuthal Equidistant with proj < 4.9.2 due to # incorrect transformation past 90 deg distance (see # https://github.com/OSGeo/proj.4/issues/246). - if PROJ4_VERSION != (): - if PROJ4_VERSION < (4, 9, 2): + if PROJ_VERSION != (): + if PROJ_VERSION < (4, 9, 2): warnings.warn('The Azimuthal Equidistant projection in Proj ' 'older than 4.9.2 incorrectly transforms points ' 'farther than 90 deg from the origin. Use this ' @@ -2283,8 +2918,8 @@ def __init__(self, central_longitude=0.0, central_latitude=0.0, super().__init__(proj4_params, globe=globe) # TODO: Let the globe return the semimajor axis always. - a = float(self.globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS) - b = float(self.globe.semiminor_axis or a) + a = float(self.ellipsoid.semi_major_metre or WGS84_SEMIMAJOR_AXIS) + b = float(self.ellipsoid.semi_minor_metre or a) coords = _ellipse_boundary(a * np.pi, b * np.pi, false_easting, false_northing, 61) @@ -2508,8 +3143,7 @@ def epsg(code): Note ---- - The conversion is performed by querying https://epsg.io/ so a - live internet connection is required. + The conversion is performed by pyproj.CRS. """ import cartopy._epsg diff --git a/lib/cartopy/img_transform.py b/lib/cartopy/img_transform.py index ae2047006..c1214e530 100644 --- a/lib/cartopy/img_transform.py +++ b/lib/cartopy/img_transform.py @@ -252,21 +252,14 @@ def regrid(array, source_x_coords, source_y_coords, source_cs, target_proj, The data array regridded in the target projection. """ - - # n.b. source_cs is actually a projection (the coord system of the - # source coordinates), but not necessarily the native projection of - # the source array (i.e. you can provide a warped image with lat lon - # coordinates). - - # XXX NB. target_x and target_y must currently be rectangular (i.e. - # be a 2d np array) - geo_cent = source_cs.as_geocentric() - xyz = geo_cent.transform_points(source_cs, - source_x_coords.flatten(), - source_y_coords.flatten()) - target_xyz = geo_cent.transform_points(target_proj, - target_x_points.flatten(), - target_y_points.flatten()) + # Stack our original xyz array, this will also wrap coords when necessary + xyz = source_cs.transform_points(source_cs, + source_x_coords.flatten(), + source_y_coords.flatten()) + # Transform the target points into the source projection + target_xyz = source_cs.transform_points(target_proj, + target_x_points.flatten(), + target_y_points.flatten()) if _is_pykdtree: kdtree = pykdtree.kdtree.KDTree(xyz) @@ -297,14 +290,11 @@ def regrid(array, source_x_coords, source_y_coords, source_cs, target_proj, # Do double transform to clip points that do not map back and forth # to the same point to within a fixed fractional offset. - # XXX THIS ONLY NEEDS TO BE DONE FOR (PSEUDO-)CYLINDRICAL PROJECTIONS - # (OR ANY OTHERS WHICH HAVE THE CONCEPT OF WRAPPING) - source_desired_xyz = source_cs.transform_points(target_proj, - target_x_points.flatten(), - target_y_points.flatten()) + # NOTE: This only needs to be done for (pseudo-)cylindrical projections, + # or any others which have the concept of wrapping back_to_target_xyz = target_proj.transform_points(source_cs, - source_desired_xyz[:, 0], - source_desired_xyz[:, 1]) + target_xyz[:, 0], + target_xyz[:, 1]) back_to_target_x = back_to_target_xyz[:, 0].reshape(desired_ny, desired_nx) back_to_target_y = back_to_target_xyz[:, 1].reshape(desired_ny, @@ -327,10 +317,10 @@ def regrid(array, source_x_coords, source_y_coords, source_cs, target_proj, # Transform the target points to the source projection and mask any points # that fall outside the original source domain. if mask_extrapolated: - target_in_source_xyz = source_cs.transform_points( - target_proj, target_x_points, target_y_points) - target_in_source_x = target_in_source_xyz[..., 0] - target_in_source_y = target_in_source_xyz[..., 1] + target_in_source_x = target_xyz[:, 0].reshape(desired_ny, + desired_nx) + target_in_source_y = target_xyz[:, 1].reshape(desired_ny, + desired_nx) bounds = _determine_bounds(source_x_coords, source_y_coords, source_cs) diff --git a/lib/cartopy/mpl/geoaxes.py b/lib/cartopy/mpl/geoaxes.py index 7081f7626..ae2e12f39 100644 --- a/lib/cartopy/mpl/geoaxes.py +++ b/lib/cartopy/mpl/geoaxes.py @@ -581,7 +581,9 @@ def format_coord(self, x, y): A string formatted for the Matplotlib GUI status bar. """ - lon, lat = ccrs.Geodetic().transform_point(x, y, self.projection) + lon, lat = self.projection.as_geodetic().transform_point( + x, y, self.projection, + ) ns = 'N' if lat >= 0.0 else 'S' ew = 'E' if lon >= 0.0 else 'W' @@ -796,7 +798,7 @@ def _get_extent_geom(self, crs=None): warnings.warn('Approximating coordinate system {!r} with a ' 'RotatedPole projection.'.format(crs)) elif hasattr(crs, 'is_geodetic') and crs.is_geodetic(): - proj = ccrs.PlateCarree(crs.globe) + proj = ccrs.PlateCarree(globe=crs.globe) warnings.warn('Approximating coordinate system {!r} with the ' 'PlateCarree projection.'.format(crs)) else: diff --git a/lib/cartopy/tests/crs/test_azimuthal_equidistant.py b/lib/cartopy/tests/crs/test_azimuthal_equidistant.py index d7ec80ffd..f8c4ce96f 100644 --- a/lib/cartopy/tests/crs/test_azimuthal_equidistant.py +++ b/lib/cartopy/tests/crs/test_azimuthal_equidistant.py @@ -21,7 +21,7 @@ def test_default(self): assert_almost_equal(np.array(aeqd.x_limits), [-20037508.34278924, 20037508.34278924], decimal=6) assert_almost_equal(np.array(aeqd.y_limits), - [-20037508.34278924, 20037508.34278924], decimal=6) + [-19970326.371123, 19970326.371123], decimal=6) def test_eccentric_globe(self): globe = ccrs.Globe(semimajor_axis=1000, semiminor_axis=500, @@ -47,7 +47,7 @@ def test_eastings(self): assert_almost_equal(np.array(aeqd_offset.x_limits), [-20036274.34278924, 20038742.34278924], decimal=6) assert_almost_equal(np.array(aeqd_offset.y_limits), - [-20041829.34278924, 20033187.34278924], decimal=6) + [-19974647.371123, 19966005.371123], decimal=6) def test_grid(self): # USGS Professional Paper 1395, pp 196--197, Table 30 @@ -156,10 +156,7 @@ def test_ellipsoid_polar_transform(self): assert_almost_equal(np.array(aeqd.x_limits), [-20038296.88254529, 20038296.88254529], decimal=6) assert_almost_equal(np.array(aeqd.y_limits), - # TODO: This is wrong. Globe.semiminor_axis does - # not account for flattening. - # [-19970827.86969727, 19970827.86969727] - [-20038296.88254529, 20038296.88254529], decimal=6) + [-19970827.86969727, 19970827.86969727], decimal=6) result = aeqd.transform_point(5.0, 80.0, geodetic) @@ -186,10 +183,7 @@ def test_ellipsoid_guam_transform(self): assert_almost_equal(np.array(aeqd.x_limits), [-19987726.36931940, 20087726.36931940], decimal=6) assert_almost_equal(np.array(aeqd.y_limits), - # TODO: This is wrong. Globe.semiminor_axis does - # not account for flattening. - # [-19919796.94787477, 20019796.94787477] - [-19987726.36931940, 20087726.36931940], decimal=6) + [-19919796.94787477, 20019796.94787477], decimal=6) pt_lat = 13 + (20 + 20.53846 / 60) / 60 pt_lon = 144 + (38 + 7.19265 / 60) / 60 @@ -220,10 +214,7 @@ def test_ellipsoid_micronesia_transform(self): assert_almost_equal(np.array(aeqd.x_limits), [-20009068.84931940, 20066383.88931940], decimal=6) assert_almost_equal(np.array(aeqd.y_limits), - # TODO: This is wrong. Globe.semiminor_axis does - # not account for flattening. - # [-19902596.95787477, 20036996.93787477] - [-19970526.37931940, 20104926.35931940], decimal=6) + [-19902596.95787477, 20036996.93787477], decimal=6) pt_lat = 15 + (14 + 47.4930 / 60) / 60 pt_lon = 145 + (47 + 34.9080 / 60) / 60 diff --git a/lib/cartopy/tests/crs/test_equal_earth.py b/lib/cartopy/tests/crs/test_equal_earth.py index 49db608c2..652c6a570 100644 --- a/lib/cartopy/tests/crs/test_equal_earth.py +++ b/lib/cartopy/tests/crs/test_equal_earth.py @@ -16,7 +16,7 @@ from .helpers import check_proj_params -pytestmark = pytest.mark.skipif(ccrs.PROJ4_VERSION < (5, 2, 0), +pytestmark = pytest.mark.skipif(ccrs.PROJ_VERSION < (5, 2, 0), reason='Proj is too old.') diff --git a/lib/cartopy/tests/crs/test_interrupted_goode_homolosine.py b/lib/cartopy/tests/crs/test_interrupted_goode_homolosine.py index 1f9d8452c..65743ad29 100644 --- a/lib/cartopy/tests/crs/test_interrupted_goode_homolosine.py +++ b/lib/cartopy/tests/crs/test_interrupted_goode_homolosine.py @@ -18,7 +18,7 @@ @pytest.mark.parametrize("emphasis", ["land", "ocean"]) def test_default(emphasis): - if emphasis == "ocean" and ccrs.PROJ4_VERSION < (7, 1, 0): + if emphasis == "ocean" and ccrs.PROJ_VERSION < (7, 1, 0): pytest.skip() igh = ccrs.InterruptedGoodeHomolosine(emphasis=emphasis) other_args = {"ellps=WGS84", "lon_0=0"} @@ -36,7 +36,7 @@ def test_default(emphasis): @pytest.mark.parametrize("emphasis", ["land", "ocean"]) def test_eccentric_globe(emphasis): - if emphasis == "ocean" and ccrs.PROJ4_VERSION < (7, 1, 0): + if emphasis == "ocean" and ccrs.PROJ_VERSION < (7, 1, 0): pytest.skip() globe = ccrs.Globe(semimajor_axis=1000, semiminor_axis=500, ellipse=None) igh = ccrs.InterruptedGoodeHomolosine(globe=globe, emphasis=emphasis) @@ -55,7 +55,7 @@ def test_eccentric_globe(emphasis): [("land", -10.0), ("land", 10.0), ("ocean", -10.0), ("ocean", 10.0)], ) def test_central_longitude(emphasis, lon): - if emphasis == "ocean" and ccrs.PROJ4_VERSION < (7, 1, 0): + if emphasis == "ocean" and ccrs.PROJ_VERSION < (7, 1, 0): pytest.skip() igh = ccrs.InterruptedGoodeHomolosine( central_longitude=lon, emphasis=emphasis diff --git a/lib/cartopy/tests/crs/test_lambert_conformal.py b/lib/cartopy/tests/crs/test_lambert_conformal.py index 209fca06f..bf7d55d04 100644 --- a/lib/cartopy/tests/crs/test_lambert_conformal.py +++ b/lib/cartopy/tests/crs/test_lambert_conformal.py @@ -69,17 +69,17 @@ def test_too_many_parallel(self): def test_single_spole(self): s_pole_crs = ccrs.LambertConformal(standard_parallels=[-1.]) assert_array_almost_equal(s_pole_crs.x_limits, - (-19840440, 19840440.), + (-19939660, 19939660), decimal=0) assert_array_almost_equal(s_pole_crs.y_limits, - (-370239953, -8191953), + (-735590302, -8183795), decimal=0) def test_single_npole(self): n_pole_crs = ccrs.LambertConformal(standard_parallels=[1.]) assert_array_almost_equal(n_pole_crs.x_limits, - (-20222156, 20222156), + (-20130569, 20130569), decimal=0) assert_array_almost_equal(n_pole_crs.y_limits, - (-8164817, 360848720), + (-8170229, 726200683), decimal=0) diff --git a/lib/cartopy/tests/crs/test_miller.py b/lib/cartopy/tests/crs/test_miller.py index 95f425f8d..264b5c6e8 100644 --- a/lib/cartopy/tests/crs/test_miller.py +++ b/lib/cartopy/tests/crs/test_miller.py @@ -18,13 +18,13 @@ def test_default(): mill = ccrs.Miller() - other_args = {'a=57.29577951308232', 'lon_0=0.0'} + other_args = {'a=6378137.0', 'lon_0=0.0'} check_proj_params('mill', mill, other_args) assert_almost_equal(np.array(mill.x_limits), - [-180, 180]) + [-20037508.3427892, 20037508.3427892]) assert_almost_equal(np.array(mill.y_limits), - [-131.9758172, 131.9758172]) + [-14691480.7691731, 14691480.7691731]) def test_sphere_globe(): @@ -75,13 +75,13 @@ def test_eccentric_globe(): @pytest.mark.parametrize('lon', [-10.0, 10.0]) def test_central_longitude(lon): mill = ccrs.Miller(central_longitude=lon) - other_args = {'a=57.29577951308232', 'lon_0={}'.format(lon)} + other_args = {'a=6378137.0', 'lon_0={}'.format(lon)} check_proj_params('mill', mill, other_args) assert_almost_equal(np.array(mill.x_limits), - [-180, 180]) + [-20037508.3427892, 20037508.3427892]) assert_almost_equal(np.array(mill.y_limits), - [-131.9758172, 131.9758172]) + [-14691480.7691731, 14691480.7691731]) def test_grid(): diff --git a/lib/cartopy/tests/crs/test_orthographic.py b/lib/cartopy/tests/crs/test_orthographic.py index d8708f398..9ce42be47 100644 --- a/lib/cartopy/tests/crs/test_orthographic.py +++ b/lib/cartopy/tests/crs/test_orthographic.py @@ -44,7 +44,7 @@ def test_ellipse_globe(): match='does not handle elliptical globes.') as w: ortho = ccrs.Orthographic(globe=globe) assert len(w) == (2 - if (5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0) + if (5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0) else 1) other_args = {'ellps=WGS84', 'lon_0=0.0', 'lat_0=0.0'} @@ -62,7 +62,7 @@ def test_eccentric_globe(): match='does not handle elliptical globes.') as w: ortho = ccrs.Orthographic(globe=globe) assert len(w) == (2 - if (5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0) + if (5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0) else 1) other_args = {'a=1000', 'b=500', 'lon_0=0.0', 'lat_0=0.0'} diff --git a/lib/cartopy/tests/crs/test_robinson.py b/lib/cartopy/tests/crs/test_robinson.py index 6569cdbbd..255ec3722 100644 --- a/lib/cartopy/tests/crs/test_robinson.py +++ b/lib/cartopy/tests/crs/test_robinson.py @@ -20,13 +20,13 @@ _CRS_ROB = ccrs.Robinson() # Increase tolerance if using older proj releases -if ccrs.PROJ4_VERSION >= (6, 3, 1): +if ccrs.PROJ_VERSION >= (6, 3, 1): _TRANSFORM_TOL = 7 -elif ccrs.PROJ4_VERSION >= (4, 9): +elif ccrs.PROJ_VERSION >= (4, 9): _TRANSFORM_TOL = 0 else: _TRANSFORM_TOL = -1 -_LIMIT_TOL = -1 # if ccrs.PROJ4_VERSION < (5, 2, 0) else 7 +_LIMIT_TOL = -1 # if ccrs.PROJ_VERSION < (5, 2, 0) else 7 def test_default(): diff --git a/lib/cartopy/tests/crs/test_rotated_geodetic.py b/lib/cartopy/tests/crs/test_rotated_geodetic.py index f5206c4fc..229d45519 100644 --- a/lib/cartopy/tests/crs/test_rotated_geodetic.py +++ b/lib/cartopy/tests/crs/test_rotated_geodetic.py @@ -12,7 +12,8 @@ from .helpers import check_proj_params -common_other_args = {'o_proj=latlon', 'to_meter=0.0174532925199433'} +common_other_args = {'o_proj=latlon', 'to_meter=111319.4907932736', + 'a=6378137.0'} def test_default(): diff --git a/lib/cartopy/tests/crs/test_rotated_pole.py b/lib/cartopy/tests/crs/test_rotated_pole.py index 118ce8246..c59dc612b 100644 --- a/lib/cartopy/tests/crs/test_rotated_pole.py +++ b/lib/cartopy/tests/crs/test_rotated_pole.py @@ -12,7 +12,7 @@ from .helpers import check_proj_params -common_other_args = {'o_proj=latlon', 'to_meter=0.0174532925199433'} +common_other_args = {'o_proj=latlon', 'to_meter=111319.4907932736'} def test_default(): diff --git a/lib/cartopy/tests/crs/test_stereographic.py b/lib/cartopy/tests/crs/test_stereographic.py index c9a041a17..2644198e4 100644 --- a/lib/cartopy/tests/crs/test_stereographic.py +++ b/lib/cartopy/tests/crs/test_stereographic.py @@ -19,9 +19,9 @@ def test_default(self): check_proj_params('stere', stereo, other_args) assert_almost_equal(np.array(stereo.x_limits), - [-5e7, 5e7], decimal=4) + [-5e7, 5e7], decimal=3) assert_almost_equal(np.array(stereo.y_limits), - [-5e7, 5e7], decimal=4) + [-5e7, 5e7], decimal=3) def test_eccentric_globe(self): globe = ccrs.Globe(semimajor_axis=1000, semiminor_axis=500, diff --git a/lib/cartopy/tests/crs/test_transverse_mercator.py b/lib/cartopy/tests/crs/test_transverse_mercator.py index 3b35d038d..aef320cf3 100644 --- a/lib/cartopy/tests/crs/test_transverse_mercator.py +++ b/lib/cartopy/tests/crs/test_transverse_mercator.py @@ -94,7 +94,7 @@ def test_default(self, approx): res = proj.transform_point(*self.point_a, src_crs=self.src_crs) np.testing.assert_array_almost_equal( res, (275614.26762651594, 386984.206429612), - decimal=0 if ccrs.PROJ4_VERSION < (5, 0, 0) else 6) + decimal=0 if ccrs.PROJ_VERSION < (5, 0, 0) else 6) def test_nan(self): proj = ccrs.OSNI(approx=True) diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_crs/mercator_squashed.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_crs/mercator_squashed.png deleted file mode 100644 index 24fa4d45a..000000000 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_crs/mercator_squashed.png and /dev/null differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label.png index 7bfac3203..b7e4518c2 100644 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label.png and b/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label.png differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label_3.4.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label_3.4.png deleted file mode 100644 index ec07cbf4e..000000000 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_examples/contour_label_3.4.png and /dev/null differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_gridliner/gridliner_labels_tight.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_gridliner/gridliner_labels_tight.png index c9ccf345c..39899fefb 100644 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_gridliner/gridliner_labels_tight.png and b/lib/cartopy/tests/mpl/baseline_images/mpl/test_gridliner/gridliner_labels_tight.png differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/global_contour_wrap_mpl_pre_3.0.0.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/global_contour_wrap_mpl_pre_3.0.0.png deleted file mode 100644 index 855b37450..000000000 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/global_contour_wrap_mpl_pre_3.0.0.png and /dev/null differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/multiple_projections4.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/multiple_projections4.png deleted file mode 100644 index 733be3c42..000000000 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/multiple_projections4.png and /dev/null differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_3.2.0.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot.png similarity index 100% rename from lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_3.2.0.png rename to lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot.png diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_2.2.2.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_2.2.2.png deleted file mode 100644 index 6a1b490c0..000000000 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_2.2.2.png and /dev/null differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_3.0.0.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_3.0.0.png deleted file mode 100644 index 33a0836ae..000000000 Binary files a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/streamplot_mpl_3.0.0.png and /dev/null differ diff --git a/lib/cartopy/tests/mpl/test_crs.py b/lib/cartopy/tests/mpl/test_crs.py index 4e3ea0cd7..583c6e03b 100644 --- a/lib/cartopy/tests/mpl/test_crs.py +++ b/lib/cartopy/tests/mpl/test_crs.py @@ -21,7 +21,7 @@ def test_igh_land(): ax.gridlines() -@pytest.mark.skipif(ccrs.PROJ4_VERSION < (7, 1, 0), reason="Proj is too old.") +@pytest.mark.skipif(ccrs.PROJ_VERSION < (7, 1, 0), reason="Proj is too old.") @pytest.mark.natural_earth @ImageTesting(["igh_ocean"]) def test_igh_ocean(): @@ -44,17 +44,6 @@ def test_lambert_south(): ax.gridlines() -@pytest.mark.natural_earth -@ImageTesting(['mercator_squashed']) -def test_mercator_squashed(): - globe = ccrs.Globe(semimajor_axis=10000, semiminor_axis=9000, - ellipse=None) - crs = ccrs.Mercator(globe=globe, min_latitude=-40, max_latitude=40) - ax = plt.axes(projection=crs) - ax.coastlines() - ax.gridlines() - - @pytest.mark.natural_earth @cleanup def test_repr_html(): diff --git a/lib/cartopy/tests/mpl/test_examples.py b/lib/cartopy/tests/mpl/test_examples.py index f2f07e1ee..2daf5a0ea 100644 --- a/lib/cartopy/tests/mpl/test_examples.py +++ b/lib/cartopy/tests/mpl/test_examples.py @@ -47,11 +47,9 @@ def test_global_map(): ax.plot([-0.08, 132], [51.53, 43.17], transform=ccrs.Geodetic()) -contour_image = 'contour_label' if MPL_VERSION < '3.4' else 'contour_label_3.4' - - @pytest.mark.natural_earth -@ExampleImageTesting([contour_image], tolerance=0) +@ExampleImageTesting(['contour_label'], + tolerance=9.9 if MPL_VERSION < "3.2" else 0) def test_contour_label(): from cartopy.tests.mpl.test_caching import sample_data fig = plt.figure() diff --git a/lib/cartopy/tests/mpl/test_gridliner.py b/lib/cartopy/tests/mpl/test_gridliner.py index 550c934d0..4c8f125f9 100644 --- a/lib/cartopy/tests/mpl/test_gridliner.py +++ b/lib/cartopy/tests/mpl/test_gridliner.py @@ -54,7 +54,7 @@ @pytest.mark.natural_earth @ImageTesting(['gridliner1'], # Robinson projection is slightly better in Proj 6+. - tolerance=0.7 if ccrs.PROJ4_VERSION >= (6, 0, 0) else 0.5) + tolerance=0.7 if ccrs.PROJ_VERSION >= (6, 0, 0) else 0.5) def test_gridliner(): ny, nx = 2, 4 @@ -129,9 +129,10 @@ def test_gridliner_specified_lines(): if MPL_VERSION < "3": TOL = 15 else: - TOL = 0.5 + TOL = 3.1 grid_label_tol = grid_label_inline_tol = grid_label_inline_usa_tol = TOL -grid_label_inline_tol += 1.1 +grid_label_tol += 0.8 +grid_label_inline_tol += 1.6 grid_label_image = 'gridliner_labels' grid_label_inline_image = 'gridliner_labels_inline' grid_label_inline_usa_image = 'gridliner_labels_inline_usa' @@ -217,7 +218,7 @@ def test_grid_labels(): @pytest.mark.skipif(geos_version == (3, 9, 0), reason="GEOS intersection bug") @pytest.mark.natural_earth @ImageTesting(['gridliner_labels_tight'], - tolerance=grid_label_tol if ccrs.PROJ4_VERSION < (7, 1, 0) + tolerance=grid_label_tol if ccrs.PROJ_VERSION < (7, 1, 0) else 4) def test_grid_labels_tight(): # Ensure tight layout accounts for gridlines @@ -268,7 +269,7 @@ def test_grid_labels_inline(): else: kwargs = {} ax = plt.subplot(7, 4, i, projection=proj(**kwargs)) - if (ccrs.PROJ4_VERSION[:2] == (5, 0) and + if (ccrs.PROJ_VERSION[:2] == (5, 0) and proj in (ccrs.Orthographic, ccrs.AlbersEqualArea, ccrs.Geostationary, ccrs.NearsidePerspective)): # Above projections are broken, so skip labels. @@ -303,7 +304,7 @@ def test_grid_labels_inline_usa(): except Exception: pass ax.set_title(proj, y=1.075) - if (ccrs.PROJ4_VERSION[:2] == (5, 0) and + if (ccrs.PROJ_VERSION[:2] == (5, 0) and proj in (ccrs.Orthographic, ccrs.AlbersEqualArea, ccrs.Geostationary, ccrs.NearsidePerspective)): # Above projections are broken, so skip labels. diff --git a/lib/cartopy/tests/mpl/test_images.py b/lib/cartopy/tests/mpl/test_images.py index 2b5676b19..0ba8ad229 100644 --- a/lib/cartopy/tests/mpl/test_images.py +++ b/lib/cartopy/tests/mpl/test_images.py @@ -35,7 +35,7 @@ # care that it is putting images onto the map which are roughly correct. @pytest.mark.natural_earth @pytest.mark.network -@pytest.mark.xfail(ccrs.PROJ4_VERSION == (5, 0, 0), +@pytest.mark.xfail(ccrs.PROJ_VERSION == (5, 0, 0), reason='Proj returns slightly different bounds.', strict=True) @ImageTesting(['web_tiles'], tolerance=5.91) @@ -74,7 +74,7 @@ def test_web_tiles(): @pytest.mark.natural_earth @pytest.mark.network -@pytest.mark.xfail(ccrs.PROJ4_VERSION == (5, 0, 0), +@pytest.mark.xfail(ccrs.PROJ_VERSION == (5, 0, 0), reason='Proj returns slightly different bounds.', strict=True) @ImageTesting(['image_merge'], tolerance=0.01) @@ -103,7 +103,7 @@ def test_image_merge(): plt.imshow(img, origin=origin, extent=extent, alpha=0.5) -@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0), +@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0), reason='Proj Orthographic projection is buggy.', strict=True) @ImageTesting(['imshow_natural_earth_ortho'], tolerance=0.7) @@ -170,7 +170,7 @@ def test_imshow_rgb(): plt.close() -@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0), +@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0), reason='Proj Orthographic projection is buggy.', strict=True) @ImageTesting(['imshow_natural_earth_ortho'], tolerance=0.7) @@ -179,7 +179,7 @@ def test_stock_img(): ax.stock_img() -@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0), +@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0), reason='Proj Orthographic projection is buggy.', strict=True) @ImageTesting(['imshow_natural_earth_ortho'], tolerance=0.7) @@ -191,7 +191,7 @@ def test_pil_Image(): extent=[-180, 180, -90, 90]) -@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0), +@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0), reason='Proj Orthographic projection is buggy.', strict=True) @ImageTesting(['imshow_natural_earth_ortho']) diff --git a/lib/cartopy/tests/mpl/test_mpl_integration.py b/lib/cartopy/tests/mpl/test_mpl_integration.py index 16753d402..2669d139f 100644 --- a/lib/cartopy/tests/mpl/test_mpl_integration.py +++ b/lib/cartopy/tests/mpl/test_mpl_integration.py @@ -4,7 +4,6 @@ # See COPYING and COPYING.LESSER in the root of the repository for full # licensing details. -import math import re import warnings @@ -17,52 +16,24 @@ from cartopy.tests.mpl import MPL_VERSION, ImageTesting -_ROB_TOL = 0.5 if ccrs.PROJ4_VERSION < (4, 9) else 0.111 -_CONTOUR_STYLE = _STREAMPLOT_STYLE = 'classic' -_CONTOUR_TOL = 0.5 -if MPL_VERSION >= '3.0.0': - _CONTOUR_IMAGE = 'global_contour_wrap' - _CONTOUR_STYLE = 'mpl20' - if MPL_VERSION < '3.2.0': - _CONTOUR_TOL = 0.74 - if MPL_VERSION >= '3.2.0': - _STREAMPLOT_IMAGE = 'streamplot_mpl_3.2.0' - else: - _STREAMPLOT_IMAGE = 'streamplot_mpl_3.0.0' - # Should have been the case for anything but _1.4.3, but we don't want to - # regenerate those images again. - _STREAMPLOT_STYLE = 'mpl20' -else: - _CONTOUR_IMAGE = 'global_contour_wrap_mpl_pre_3.0.0' - _STREAMPLOT_IMAGE = 'streamplot_mpl_2.2.2' - - -@pytest.mark.natural_earth -@ImageTesting([_CONTOUR_IMAGE], style=_CONTOUR_STYLE, tolerance=_CONTOUR_TOL) +@pytest.mark.natural_earth +@ImageTesting(['global_contour_wrap'], style='mpl20') def test_global_contour_wrap_new_transform(): ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines() x, y = np.meshgrid(np.linspace(0, 360), np.linspace(-90, 90)) data = np.sin(np.sqrt(x ** 2 + y ** 2)) plt.contour(x, y, data, transform=ccrs.PlateCarree()) - if MPL_VERSION <= '3.0.0': - # Rather than updating image test for old version - # Remove when only MPL > 3 is required - ax.set_extent((-176.3265306122449, 176.3265306122449, -90.0, 90.0)) @pytest.mark.natural_earth -@ImageTesting([_CONTOUR_IMAGE], style=_CONTOUR_STYLE, tolerance=_CONTOUR_TOL) +@ImageTesting(['global_contour_wrap'], style='mpl20') def test_global_contour_wrap_no_transform(): ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines() x, y = np.meshgrid(np.linspace(0, 360), np.linspace(-90, 90)) data = np.sin(np.sqrt(x ** 2 + y ** 2)) plt.contour(x, y, data) - if MPL_VERSION <= '3.0.0': - # Rather than updating image test for old version - # Remove when only MPL > 3 is required - ax.set_extent((-176.3265306122449, 176.3265306122449, -90.0, 90.0)) @pytest.mark.natural_earth @@ -130,7 +101,7 @@ def test_global_scatter_wrap_no_transform(): @pytest.mark.natural_earth @ImageTesting(['global_hexbin_wrap'], - tolerance=2 if MPL_VERSION < '3' else 0.5) + tolerance=2 if MPL_VERSION < '3.2' else 0.5) def test_global_hexbin_wrap(): ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines(zorder=2) @@ -147,7 +118,7 @@ def test_global_hexbin_wrap(): @pytest.mark.natural_earth @ImageTesting(['global_hexbin_wrap'], - tolerance=2 if MPL_VERSION < '3' else 0.5) + tolerance=2 if MPL_VERSION < '3.2' else 0.5) def test_global_hexbin_wrap_transform(): ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines(zorder=2) @@ -191,9 +162,8 @@ def test_simple_global(): @pytest.mark.filterwarnings("ignore:Unable to determine extent") @pytest.mark.natural_earth -@ImageTesting(['multiple_projections4' if ccrs.PROJ4_VERSION < (5, 0, 0) - else 'multiple_projections5'], - tolerance=0.81) +@ImageTesting(['multiple_projections5'], + tolerance=2.05) def test_multiple_projections(): projections = [ccrs.PlateCarree(), @@ -201,9 +171,7 @@ def test_multiple_projections(): ccrs.RotatedPole(pole_latitude=45, pole_longitude=180), ccrs.OSGB(approx=True), ccrs.TransverseMercator(approx=True), - ccrs.Mercator( - globe=ccrs.Globe(semimajor_axis=math.degrees(1)), - min_latitude=-85., max_latitude=85.), + ccrs.Mercator(min_latitude=-85., max_latitude=85.), ccrs.LambertCylindrical(), ccrs.Miller(), ccrs.Gnomonic(), @@ -232,15 +200,13 @@ def test_multiple_projections(): ax.coastlines(resolution="110m") plt.plot(-0.08, 51.53, 'o', transform=ccrs.PlateCarree()) - plt.plot([-0.08, 132], [51.53, 43.17], color='red', transform=ccrs.PlateCarree()) - plt.plot([-0.08, 132], [51.53, 43.17], color='blue', - transform=ccrs.Geodetic()) + transform=prj.as_geodetic()) -@pytest.mark.skipif(ccrs.PROJ4_VERSION < (5, 2, 0), +@pytest.mark.skipif(ccrs.PROJ_VERSION < (5, 2, 0), reason='Proj is too old.') @pytest.mark.natural_earth @ImageTesting(['multiple_projections520']) @@ -287,7 +253,7 @@ def test_cursor_values(): @pytest.mark.natural_earth -@ImageTesting(['natural_earth_interface'], tolerance=_ROB_TOL) +@ImageTesting(['natural_earth_interface'], tolerance=0.111) def test_axes_natural_earth_interface(): rob = ccrs.Robinson() @@ -374,8 +340,8 @@ def test_pcolormesh_get_array_with_mask(): 'Data supplied does not match data retrieved in unwrapped case' -tolerance = 1.61 -if (5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0): +tolerance = 1.87 +if (5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0): tolerance += 0.8 @@ -411,7 +377,7 @@ def test_pcolormesh_global_with_wrap2(): tolerance = 1.39 -if (5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0): +if (5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0): tolerance += 1.4 @@ -747,7 +713,7 @@ def test_quiver_regrid(): @pytest.mark.natural_earth -@ImageTesting(['quiver_regrid_with_extent'], tolerance=0.53) +@ImageTesting(['quiver_regrid_with_extent'], tolerance=0.54) def test_quiver_regrid_with_extent(): x = np.arange(-60, 42.5, 2.5) y = np.arange(30, 72.5, 2.5) @@ -858,7 +824,8 @@ def test_barbs_1d_transformed(): @pytest.mark.natural_earth -@ImageTesting([_STREAMPLOT_IMAGE], style=_STREAMPLOT_STYLE, tolerance=0.54) +@ImageTesting(['streamplot'], style='mpl20', + tolerance=42 if MPL_VERSION < "3.2" else 0.54) def test_streamplot(): x = np.arange(-60, 42.5, 2.5) y = np.arange(30, 72.5, 2.5) diff --git a/lib/cartopy/tests/mpl/test_pseudo_color.py b/lib/cartopy/tests/mpl/test_pseudo_color.py index b08436566..762bd6ec7 100644 --- a/lib/cartopy/tests/mpl/test_pseudo_color.py +++ b/lib/cartopy/tests/mpl/test_pseudo_color.py @@ -14,13 +14,13 @@ def test_pcolormesh_partially_masked(): - data = np.ma.masked_all((30, 40)) + data = np.ma.masked_all((40, 30)) data[0:100] = 10 # Check that a partially masked data array does trigger a pcolor call. with mock.patch('cartopy.mpl.geoaxes.GeoAxes.pcolor') as pcolor: ax = plt.axes(projection=ccrs.PlateCarree()) - ax.pcolormesh(np.linspace(-90, 90, 40), np.linspace(0, 360, 30), data) + ax.pcolormesh(np.linspace(0, 360, 30), np.linspace(-90, 90, 40), data) assert pcolor.call_count == 1, ("pcolor should have been called " "exactly once.") plt.close() diff --git a/lib/cartopy/tests/mpl/test_ticks.py b/lib/cartopy/tests/mpl/test_ticks.py index ad0f372fa..d79d55f5c 100644 --- a/lib/cartopy/tests/mpl/test_ticks.py +++ b/lib/cartopy/tests/mpl/test_ticks.py @@ -3,40 +3,14 @@ # This file is part of Cartopy and is released under the LGPL license. # See COPYING and COPYING.LESSER in the root of the repository for full # licensing details. - -import math - import matplotlib.pyplot as plt -import matplotlib.ticker import pytest import cartopy.crs as ccrs +from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter from cartopy.tests.mpl import ImageTesting -def _format_lat(val, i): - if val > 0: - return '%.0fN' % val - elif val < 0: - return '%.0fS' % abs(val) - else: - return '0' - - -def _format_lon(val, i): - # Apply periodic boundary conditions, with an almost equal test on 180 lon. - while val > 180: - val -= 360 - while val < -180: - val += 360 - if abs(abs(val) - 180.) <= 1e-06 or val == 0: - return '%.0f' % abs(val) - elif val > 0: - return '%.0fE' % val - elif val < 0: - return '%.0fW' % abs(val) - - ticks_tolerance = 7 @@ -46,8 +20,7 @@ def _format_lon(val, i): def test_set_xticks_no_transform(): ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines('110m') - ax.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lon)) - ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lat)) + ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol='')) ax.set_xticks([-180, -90, 0, 90, 180]) ax.set_xticks([-135, -45, 45, 135], minor=True) @@ -58,11 +31,9 @@ def test_set_xticks_no_transform(): def test_set_xticks_cylindrical(): ax = plt.axes(projection=ccrs.Mercator( min_latitude=-85., - max_latitude=85., - globe=ccrs.Globe(semimajor_axis=math.degrees(1)))) + max_latitude=85.)) ax.coastlines('110m') - ax.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lon)) - ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lat)) + ax.xaxis.set_major_formatter(LongitudeFormatter(degree_symbol='')) ax.set_xticks([-180, -90, 0, 90, 180], crs=ccrs.PlateCarree()) ax.set_xticks([-135, -45, 45, 135], minor=True, crs=ccrs.PlateCarree()) @@ -82,8 +53,7 @@ def test_set_xticks_non_cylindrical(): def test_set_yticks_no_transform(): ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines('110m') - ax.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lon)) - ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lat)) + ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol='')) ax.set_yticks([-60, -30, 0, 30, 60]) ax.set_yticks([-75, -45, 15, 45, 75], minor=True) @@ -94,11 +64,9 @@ def test_set_yticks_no_transform(): def test_set_yticks_cylindrical(): ax = plt.axes(projection=ccrs.Mercator( min_latitude=-85., - max_latitude=85., - globe=ccrs.Globe(semimajor_axis=math.degrees(1)))) + max_latitude=85.)) ax.coastlines('110m') - ax.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lon)) - ax.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(_format_lat)) + ax.yaxis.set_major_formatter(LatitudeFormatter(degree_symbol='')) ax.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree()) ax.set_yticks([-75, -45, 15, 45, 75], minor=True, crs=ccrs.PlateCarree()) @@ -117,8 +85,7 @@ def test_set_yticks_non_cylindrical(): def test_set_xyticks(): fig = plt.figure(figsize=(10, 10)) projections = (ccrs.PlateCarree(), - ccrs.Mercator(globe=ccrs.Globe( - semimajor_axis=math.degrees(1))), + ccrs.Mercator(), ccrs.TransverseMercator(approx=False)) x = -3.275024 y = 50.753998 diff --git a/lib/cartopy/tests/mpl/test_web_services.py b/lib/cartopy/tests/mpl/test_web_services.py index 92525b594..aff0f0180 100644 --- a/lib/cartopy/tests/mpl/test_web_services.py +++ b/lib/cartopy/tests/mpl/test_web_services.py @@ -37,7 +37,7 @@ def test_wms_tight_layout(): @pytest.mark.network -@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ4_VERSION < (5, 1, 0), +@pytest.mark.xfail((5, 0, 0) <= ccrs.PROJ_VERSION < (5, 1, 0), reason='Proj Orthographic projection is buggy.', strict=True) @pytest.mark.skipif(not _OWSLIB_AVAILABLE, reason='OWSLib is unavailable.') diff --git a/lib/cartopy/tests/test_crs.py b/lib/cartopy/tests/test_crs.py index d93636475..8871816e9 100644 --- a/lib/cartopy/tests/test_crs.py +++ b/lib/cartopy/tests/test_crs.py @@ -11,10 +11,7 @@ import numpy as np from numpy.testing import assert_almost_equal, assert_array_equal from numpy.testing import assert_array_almost_equal as assert_arr_almost_eq -try: - import pyepsg -except ImportError: - pyepsg = None +import pyproj import pytest import shapely.geometry as sgeom @@ -53,7 +50,7 @@ def _check_osgb(self, osgb): # results obtained by streetmap.co.uk. lat, lon = np.array([50.462023, -3.478831], dtype=np.double) - east, north = np.array([295131, 63511], dtype=np.double) + east, north = np.array([295132.1, 63512.6], dtype=np.double) # note the handling of precision here... assert_arr_almost_eq(np.array(osgb.transform_point(lon, lat, ll)), @@ -75,20 +72,23 @@ def _check_osgb(self, osgb): def test_osgb(self, approx): self._check_osgb(ccrs.OSGB(approx=approx)) - @pytest.mark.network - @pytest.mark.skipif(pyepsg is None, reason='requires pyepsg') def test_epsg(self): uk = ccrs.epsg(27700) assert uk.epsg_code == 27700 - assert_almost_equal(uk.x_limits, (-118365.7406176, 751581.5647514), - decimal=3) - assert_almost_equal(uk.y_limits, (-5268.1704980, 1272227.7987656), - decimal=2) - assert_almost_equal(uk.threshold, 8699.47, decimal=2) + if ccrs.PROJ_VERSION >= (8, 0, 0): + assert_almost_equal(uk.x_limits, (-104009.357, 688806.007), + decimal=3) + assert_almost_equal(uk.y_limits, (-8908.37, 1256558.45), + decimal=2) + assert_almost_equal(uk.threshold, 7928.15, decimal=2) + else: + assert_almost_equal(uk.x_limits, (-118397.001, 751441.779), + decimal=3) + assert_almost_equal(uk.y_limits, (-5192.07, 1272149.35), + decimal=2) + assert_almost_equal(uk.threshold, 8698.39, decimal=2) self._check_osgb(uk) - @pytest.mark.network - @pytest.mark.skipif(pyepsg is None, reason='requires pyepsg') def test_epsg_compound_crs(self): projection = ccrs.epsg(5973) assert projection.epsg_code == 5973 @@ -180,8 +180,12 @@ def test_globe(self): rugby_moll = ccrs.Mollweide(globe=rugby_globe) footy_moll = ccrs.Mollweide(globe=footy_globe) - rugby_pt = rugby_moll.transform_point(10, 10, ccrs.Geodetic()) - footy_pt = footy_moll.transform_point(10, 10, ccrs.Geodetic()) + rugby_pt = rugby_moll.transform_point( + 10, 10, rugby_moll.as_geodetic(), + ) + footy_pt = footy_moll.transform_point( + 10, 10, footy_moll.as_geodetic(), + ) assert_arr_almost_eq(rugby_pt, (1400915, 1741319), decimal=0) assert_arr_almost_eq(footy_pt, (155657, 193479), decimal=0) @@ -294,6 +298,11 @@ def test_transform_points_outside_domain(): crs = ccrs.Orthographic() result = crs.transform_points(ccrs.PlateCarree(), np.array([-120]), np.array([80])) + assert np.all(np.isnan(result[..., :2])) + assert result[..., -1] == 0 + result = crs.transform_points(ccrs.PlateCarree(), + np.array([-120]), np.array([80]), + trap=True) assert np.all(np.isnan(result)) # A length-2 array of the same transform produces "inf" rather # than nan due to PROJ never returning nan itself. @@ -305,3 +314,16 @@ def test_transform_points_outside_domain(): # the same as the transform_points call with a length-1 array result = crs.transform_point(-120, 80, ccrs.PlateCarree()) assert np.all(np.isnan(result)) + + +def test_projection__from_string(): + crs = ccrs.Projection("NAD83 / Pennsylvania South") + assert crs.as_geocentric().datum.name == "North American Datum 1983" + assert_almost_equal( + crs.bounds, + [361633.1351868, 859794.6690229, 45575.5693199, 209415.9845754], + ) + + +def test_crs__from_pyproj_crs(): + assert ccrs.CRS(pyproj.CRS("EPSG:4326")) == "EPSG:4326" diff --git a/lib/cartopy/tests/test_img_tiles.py b/lib/cartopy/tests/test_img_tiles.py index 7c7448ac4..3f3cc77e1 100644 --- a/lib/cartopy/tests/test_img_tiles.py +++ b/lib/cartopy/tests/test_img_tiles.py @@ -31,7 +31,7 @@ (8, 9, 4): (0, 2504688.542848654, -5009377.085697312, -2504688.542848654), } -if ccrs.PROJ4_VERSION == (5, 0, 0): +if ccrs.PROJ_VERSION == (5, 0, 0): KNOWN_EXTENTS = { (0, 0, 0): (-20037508.342789244, 20037508.342789244, -19994827.892149, 19994827.892149), @@ -158,7 +158,7 @@ def test_image_for_domain(): ll_extent = ccrs.Geodetic().transform_points(gt.crs, np.array(extent[:2]), np.array(extent[2:])) - if ccrs.PROJ4_VERSION == (5, 0, 0): + if ccrs.PROJ_VERSION == (5, 0, 0): assert_arr_almost(ll_extent[:, :2], [[-11.25, 49.033955], [11.25, 61.687101]]) diff --git a/lib/cartopy/tests/test_img_transform.py b/lib/cartopy/tests/test_img_transform.py index 586441655..4ddc0804c 100644 --- a/lib/cartopy/tests/test_img_transform.py +++ b/lib/cartopy/tests/test_img_transform.py @@ -43,7 +43,7 @@ def test_mesh_projection_extent(xmin, xmax, ymin, ymax): assert_array_equal(np.diff(target_y, axis=0), (ymax - ymin) / ny) -def test_griding_data_std_range(): +def test_gridding_data_std_range(): # Data which exists inside the standard projection bounds i.e. # [-180, 180]. target_prj = ccrs.PlateCarree() @@ -76,7 +76,7 @@ def test_griding_data_std_range(): assert_array_equal(expected_mask, image.mask) -def test_griding_data_outside_projection(): +def test_gridding_data_outside_projection(): # Data which exists outside the standard projection e.g. [0, 360] rather # than [-180, 180]. target_prj = ccrs.PlateCarree() @@ -94,10 +94,10 @@ def test_griding_data_outside_projection(): # The expected image. n.b. on a map the data is reversed in the y axis. expected = np.array( - [[3, 3, 3, 3, 3, 3, 3, 3], - [3, 3, 3, 3, 3, 1, 2, 2], - [2, 2, 3, 1, 1, 1, 1, 2], - [1, 1, 1, 1, 1, 1, 1, 1]], dtype=np.float64) + [[3, 3, 3, 3, 3, 2, 2, 2], + [3, 3, 3, 3, 1, 1, 2, 2], + [3, 3, 3, 3, 1, 1, 1, 2], + [3, 3, 3, 1, 1, 1, 1, 1]], dtype=np.float64) expected_mask = np.array( [[True, True, True, True, True, True, True, True], diff --git a/lib/cartopy/tests/test_linear_ring.py b/lib/cartopy/tests/test_linear_ring.py index f0f4bf582..707a07702 100644 --- a/lib/cartopy/tests/test_linear_ring.py +++ b/lib/cartopy/tests/test_linear_ring.py @@ -175,4 +175,4 @@ def test_at_boundary(self): # Test area of smallest Polygon that contains all the points in the # geometry. - assert round(abs(mlinestr.convex_hull.area - 2347.75623076), 7) == 0 + assert round(abs(mlinestr.convex_hull.area - 2347.7562), 4) == 0 diff --git a/lib/cartopy/tests/test_polygon.py b/lib/cartopy/tests/test_polygon.py index 63bc61c6e..ebf2557d0 100644 --- a/lib/cartopy/tests/test_polygon.py +++ b/lib/cartopy/tests/test_polygon.py @@ -289,7 +289,7 @@ def setup_class(self): pole_latitude=37.5) polygon = sgeom.Polygon([ (177.5, -57.38460319), - (180.0, -57.445077), + (180.1, -57.445077), (175.0, -57.19913331), ]) self.multi_polygon = projection.project_geometry(polygon) diff --git a/lib/cartopy/tests/test_vector_transform.py b/lib/cartopy/tests/test_vector_transform.py index f9a71294f..61956cd8a 100644 --- a/lib/cartopy/tests/test_vector_transform.py +++ b/lib/cartopy/tests/test_vector_transform.py @@ -143,10 +143,10 @@ def test_with_transform(self): expected_y_grid = np.array([[5., 5., 5., 5., 5.], [7.5, 7.5, 7.5, 7.5, 7.5], [10., 10., 10., 10., 10]]) - expected_u_grid = np.array([[np.nan, np.nan, 3, np.nan, np.nan], + expected_u_grid = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan], [np.nan, 2.3838, 3.5025, 2.6152, np.nan], [2, 3.0043, 4, 2.9022, 2]]) - expected_v_grid = np.array([[np.nan, np.nan, 0.3, np.nan, np.nan], + expected_v_grid = np.array([[np.nan, np.nan, np.nan, np.nan, np.nan], [np.nan, 2.6527, 2.1904, 2.4192, np.nan], [5.5, 4.6483, 4, 4.47, 5.5]]) diff --git a/lib/cartopy/trace.pyx b/lib/cartopy/trace.pyx index 7dcba5332..e741e3343 100644 --- a/lib/cartopy/trace.pyx +++ b/lib/cartopy/trace.pyx @@ -13,6 +13,7 @@ In general, this should never be called manually, instead leaving the processing to be done by the :class:`cartopy.crs.Projection` subclasses. """ from __future__ import print_function +from functools import lru_cache cimport cython from libc.math cimport HUGE_VAL, sqrt @@ -48,18 +49,29 @@ cdef extern from "geos_c.h": void GEOSPreparedGeom_destroy_r(GEOSContextHandle_t handle, const GEOSPreparedGeometry* g) nogil cdef int GEOS_MULTILINESTRING -from cartopy._crs cimport CRS -from cartopy._crs import PROJ4_VERSION -from ._proj4 cimport (projPJ, projLP, pj_get_spheroid_defn, pj_transform, - pj_strerrno, DEG_TO_RAD) from .geodesic._geodesic cimport (geod_geodesic, geod_geodesicline, geod_init, geod_geninverse, geod_lineinit, geod_genposition, GEOD_ARCMODE, GEOD_LATITUDE, 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.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(): @@ -160,73 +172,81 @@ cdef class LineAccumulator: cdef class Interpolator: cdef Point start cdef Point end - cdef projPJ src_proj - cdef projPJ dest_proj + cdef readonly transformer cdef double src_scale cdef double dest_scale + cdef bint to_180 def __cinit__(self): self.src_scale = 1 self.dest_scale = 1 + self.to_180 = False - cdef void init(self, projPJ src_proj, projPJ dest_proj): - self.src_proj = src_proj - self.dest_proj = dest_proj + cdef void init(self, src_crs, dest_crs) except *: + self.transformer = Transformer.from_crs(src_crs, dest_crs, always_xy=True) + self.to_180 = ( + self.transformer.name == "noop" and + src_crs.__class__.__name__ in ("PlateCarree", "RotatedPole") + ) cdef void set_line(self, const Point &start, const Point &end): self.start = start self.end = end - cdef Point interpolate(self, double t): - raise NotImplementedError + cdef Point project(self, const Point &src_xy) except *: + cdef Point dest_xy + + try: + xx, yy = self.transformer.transform( + src_xy.x * self.src_scale, + src_xy.y * self.src_scale, + 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 + ): + xx = HUGE_VAL + yy = HUGE_VAL + else: + raise + + if self.to_180 and xx > 180 and xx != HUGE_VAL: + xx = (((xx + 180) % 360) - 180) + + dest_xy.x = xx * self.dest_scale + dest_xy.y = yy * self.dest_scale + return dest_xy - cdef Point project(self, const Point &point): + cdef Point interpolate(self, double t) except *: raise NotImplementedError cdef class CartesianInterpolator(Interpolator): - cdef Point interpolate(self, double t): + cdef Point interpolate(self, double t) except *: cdef Point xy xy.x = self.start.x + (self.end.x - self.start.x) * t xy.y = self.start.y + (self.end.y - self.start.y) * t return self.project(xy) - cdef Point project(self, const Point &src_xy): - cdef Point dest_xy - cdef projLP xy - - xy.u = src_xy.x * self.src_scale - xy.v = src_xy.y * self.src_scale - - cdef int status = pj_transform(self.src_proj, self.dest_proj, - 1, 1, &xy.u, &xy.v, NULL) - if status in (-14, -20): - # -14 => "latitude or longitude exceeded limits" - # -20 => "tolerance condition error" - xy.u = xy.v = HUGE_VAL - elif status != 0: - raise Exception('pj_transform failed: %d\n%s' % ( - status, - pj_strerrno(status))) - - dest_xy.x = xy.u * self.dest_scale - dest_xy.y = xy.v * self.dest_scale - return dest_xy - cdef class SphericalInterpolator(Interpolator): cdef geod_geodesic geod cdef geod_geodesicline geod_line cdef double a13 - cdef void init(self, projPJ src_proj, projPJ dest_proj): - self.src_proj = src_proj - self.dest_proj = dest_proj + cdef void init(self, src_crs, dest_crs) except *: + self.transformer = Transformer.from_crs(src_crs, dest_crs, always_xy=True) - cdef double major_axis - cdef double eccentricity_squared - pj_get_spheroid_defn(self.src_proj, &major_axis, &eccentricity_squared) - geod_init(&self.geod, major_axis, 1 - sqrt(1 - eccentricity_squared)) + cdef double major_axis = src_crs.ellipsoid.semi_major_metre + 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) cdef void set_line(self, const Point &start, const Point &end): cdef double azi1 @@ -236,7 +256,7 @@ cdef class SphericalInterpolator(Interpolator): geod_lineinit(&self.geod_line, &self.geod, start.y, start.x, azi1, GEOD_LATITUDE | GEOD_LONGITUDE); - cdef Point interpolate(self, double t): + cdef Point interpolate(self, double t) except *: cdef Point lonlat geod_genposition(&self.geod_line, GEOD_ARCMODE, self.a13 * t, @@ -245,28 +265,6 @@ cdef class SphericalInterpolator(Interpolator): return self.project(lonlat) - cdef Point project(self, const Point &lonlat): - cdef Point xy - cdef projLP dest - - dest.u = (lonlat.x * DEG_TO_RAD) * self.src_scale - dest.v = (lonlat.y * DEG_TO_RAD) * self.src_scale - - cdef int status = pj_transform(self.src_proj, self.dest_proj, - 1, 1, &dest.u, &dest.v, NULL) - if status in (-14, -20): - # -14 => "latitude or longitude exceeded limits" - # -20 => "tolerance condition error" - dest.u = dest.v = HUGE_VAL - elif status != 0: - raise Exception('pj_transform failed: %d\n%s' % ( - status, - pj_strerrno(status))) - - xy.x = dest.u * self.dest_scale - xy.y = dest.v * self.dest_scale - return xy - cdef enum State: POINT_IN = 1, @@ -301,7 +299,7 @@ cdef bool straightAndDomain(double t_start, const Point &p_start, Interpolator interpolator, double threshold, GEOSContextHandle_t handle, const GEOSPreparedGeometry *gp_domain, - bool inside): + bool inside) except *: """ Return whether the given line segment is suitable as an approximation of the projection of the source line. @@ -428,7 +426,7 @@ cdef void bisect(double t_start, const Point &p_start, const Point &p_end, GEOSContextHandle_t handle, const GEOSPreparedGeometry *gp_domain, const State &state, Interpolator interpolator, double threshold, - double &t_min, Point &p_min, double &t_max, Point &p_max): + double &t_min, Point &p_min, double &t_max, Point &p_max) except *: cdef double t_current cdef Point p_current cdef bool valid @@ -483,7 +481,7 @@ cdef void _project_segment(GEOSContextHandle_t handle, unsigned int src_idx_from, unsigned int src_idx_to, Interpolator interpolator, const GEOSPreparedGeometry *gp_domain, - double threshold, LineAccumulator lines): + double threshold, LineAccumulator lines) except *: cdef Point p_current, p_min, p_max, p_end cdef double t_current, t_min, t_max cdef State state @@ -562,7 +560,8 @@ cdef void _project_segment(GEOSContextHandle_t handle, lines.new_line() -cdef _interpolator(CRS src_crs, CRS dest_projection): +@lru_cache(maxsize=4) +def _interpolator(src_crs, dest_projection): # Get an Interpolator from the given CRS and projection. # Callers must hold a reference to these systems for the lifetime # of the interpolator. If they get garbage-collected while interpolator @@ -573,8 +572,8 @@ cdef _interpolator(CRS src_crs, CRS dest_projection): interpolator = SphericalInterpolator() else: interpolator = CartesianInterpolator() - interpolator.init(src_crs.proj4, (dest_projection).proj4) - if (6, 1, 1) <= PROJ4_VERSION < (6, 3, 0): + interpolator.init(src_crs, dest_projection) + if (6, 1, 1) <= PROJ_VERSION < (6, 3, 0): # Workaround bug in Proj 6.1.1+ with +to_meter on +proj=ob_tran. # See https://github.com/OSGeo/proj#1782. lonlat = ('latlon', 'latlong', 'lonlat', 'longlat') @@ -589,7 +588,7 @@ cdef _interpolator(CRS src_crs, CRS dest_projection): return interpolator -def project_linear(geometry not None, CRS src_crs not None, +def project_linear(geometry not None, src_crs not None, dest_projection not None): """ Project a geometry from one projection to another. diff --git a/requirements/default.txt b/requirements/default.txt index 16ed4e653..aab0c4ceb 100644 --- a/requirements/default.txt +++ b/requirements/default.txt @@ -1,3 +1,4 @@ numpy>=1.13.3 shapely>=1.5.6 pyshp>=2 +pyproj>=3.0.0 diff --git a/requirements/epsg.txt b/requirements/epsg.txt deleted file mode 100644 index cb1db080b..000000000 --- a/requirements/epsg.txt +++ /dev/null @@ -1 +0,0 @@ -pyepsg>=0.4.0 diff --git a/setup.py b/setup.py index b6a37d700..bde2cfec5 100644 --- a/setup.py +++ b/setup.py @@ -280,10 +280,6 @@ def get_config_var(name): cython_coverage_enabled = os.environ.get('CYTHON_COVERAGE', None) -if proj_version >= (6, 0, 0): - extra_extension_args["define_macros"].append( - ('ACCEPT_USE_OF_DEPRECATED_PROJ_API_H', '1') - ) if cython_coverage_enabled: extra_extension_args["define_macros"].append( ('CYTHON_TRACE_NOGIL', '1') @@ -299,13 +295,6 @@ def get_config_var(name): library_dirs=[library_dir] + proj_library_dirs + geos_library_dirs, language='c++', **extra_extension_args), - Extension( - 'cartopy._crs', - ['lib/cartopy/_crs.pyx'], - include_dirs=[include_dir, np.get_include()] + proj_includes, - libraries=proj_libraries, - library_dirs=[library_dir] + proj_library_dirs, - **extra_extension_args), # Requires proj v4.9 Extension( 'cartopy.geodesic._geodesic',