ENH: Added Proj.get_factors() (#548)

Author: snowman2

docs/api/proj.rst

@@ -18,3 +18,10 @@ pyproj.Proj
     :inherited-members:
     :special-members: __init__, __call__
     :show-inheritance:
+
+
+pyproj.proj.Factors
+-------------------
+
+.. autoclass:: pyproj.proj.Factors
+    :members:

docs/history.rst

@@ -1,6 +1,10 @@
 Change Log
 ==========
 
+2.6.0
+~~~~~
+* ENH: Added Proj.get_factors() (issue #503) 
+
 2.5.0
 ~~~~~
 * Wheels contain PROJ version is 6.3.1

pyproj/_proj.pxd

@@ -1,6 +1,6 @@
 include "proj.pxi"
 
-cdef class Proj:
+cdef class _Proj:
     cdef PJ * projobj
     cdef PJ_CONTEXT* context
     cdef PJ_PROJ_INFO projobj_info

pyproj/_proj.pyx

@@ -1,21 +1,76 @@
 include "base.pxi"
 
+import copy
 import warnings
+from collections import namedtuple
 
 cimport cython
 
 from pyproj._datadir cimport pyproj_context_initialize
 from pyproj.compat import cstrencode, pystrdecode
 from pyproj.exceptions import ProjError
 
-# # version number string for PROJ
+# version number string for PROJ
 proj_version_str = "{0}.{1}.{2}".format(
     PROJ_VERSION_MAJOR,
     PROJ_VERSION_MINOR,
     PROJ_VERSION_PATCH
 )
 
-cdef class Proj:
+
+Factors = namedtuple("Factors",
+    [
+        "meridional_scale",
+        "parallel_scale",
+        "areal_scale",
+        "angular_distortion",
+        "meridian_parallel_angle",
+        "meridian_convergence",
+        "tissot_semimajor",
+        "tissot_semiminor",
+        "dx_dlam",
+        "dx_dphi",
+        "dy_dlam",
+        "dy_dphi",
+    ]
+)
+Factors.__doc__ = """
+.. versionadded:: 2.6.0
+
+These are the scaling and angular distortion factors.
+
+See `PJ_FACTORS documentation <https://proj.org/development/reference/datatypes.html?highlight=pj_factors#c.PJ_FACTORS>`__
+
+Parameters
+----------
+meridional_scale: List[float]
+     Meridional scale at coordinate.
+parallel_scale: List[float]
+    Parallel scale at coordinate.
+areal_scale: List[float]
+    Areal scale factor at coordinate.
+angular_distortion: List[float]
+    Angular distortion at coordinate.
+meridian_parallel_angle: List[float]
+    Meridian/parallel angle at coordinate.
+meridian_convergence: List[float]
+    Meridian convergence at coordinate. Sometimes also described as *grid declination*.
+tissot_semimajor: List[float]
+    Maximum scale factor.
+tissot_semiminor: List[float]
+    Minimum scale factor.
+dx_dlam: List[float]
+    Partial derivative of coordinate.
+dx_dphi: List[float]
+    Partial derivative of coordinate.
+dy_dlam: List[float]
+    Partial derivative of coordinate.
+dy_dphi: List[float]
+    Partial derivative of coordinate.
+"""
+
+
+cdef class _Proj:
     def __cinit__(self):
         self.projobj = NULL
         self.context = NULL
@@ -187,19 +242,134 @@ cdef class Proj:
                     ybuff.data[iii] = projlonlatout.uv.v
         ProjError.clear()
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _get_factors(self, longitude, latitude, bint radians, bint errcheck):
+        """
+        Calculates the projection factors PJ_FACTORS
+
+        Equivalent to `proj -S` command line.
+        """
+        cdef PyBuffWriteManager lonbuff = PyBuffWriteManager(longitude)
+        cdef PyBuffWriteManager latbuff = PyBuffWriteManager(longitude)
+
+        if not lonbuff.len or not (lonbuff.len == latbuff.len):
+            raise ProjError('longitude and latitude must be same size')
+
+        # prepare the factors output
+        meridional_scale = copy.copy(longitude)
+        parallel_scale = copy.copy(longitude)
+        areal_scale = copy.copy(longitude)
+        angular_distortion = copy.copy(longitude)
+        meridian_parallel_angle = copy.copy(longitude)
+        meridian_convergence = copy.copy(longitude)
+        tissot_semimajor = copy.copy(longitude)
+        tissot_semiminor = copy.copy(longitude)
+        dx_dlam = copy.copy(longitude)
+        dx_dphi = copy.copy(longitude)
+        dy_dlam = copy.copy(longitude)
+        dy_dphi = copy.copy(longitude)
+        cdef PyBuffWriteManager meridional_scale_buff = PyBuffWriteManager(meridional_scale)
+        cdef PyBuffWriteManager parallel_scale_buff = PyBuffWriteManager(parallel_scale)
+        cdef PyBuffWriteManager areal_scale_buff = PyBuffWriteManager(areal_scale)
+        cdef PyBuffWriteManager angular_distortion_buff = PyBuffWriteManager(angular_distortion)
+        cdef PyBuffWriteManager meridian_parallel_angle_buff = PyBuffWriteManager(meridian_parallel_angle)
+        cdef PyBuffWriteManager meridian_convergence_buff = PyBuffWriteManager(meridian_convergence)
+        cdef PyBuffWriteManager tissot_semimajor_buff = PyBuffWriteManager(tissot_semimajor)
+        cdef PyBuffWriteManager tissot_semiminor_buff = PyBuffWriteManager(tissot_semiminor)
+        cdef PyBuffWriteManager dx_dlam_buff = PyBuffWriteManager(dx_dlam)
+        cdef PyBuffWriteManager dx_dphi_buff = PyBuffWriteManager(dx_dphi)
+        cdef PyBuffWriteManager dy_dlam_buff = PyBuffWriteManager(dy_dlam)
+        cdef PyBuffWriteManager dy_dphi_buff = PyBuffWriteManager(dy_dphi)
+
+        # calculate the factors
+        cdef PJ_COORD pj_coord = proj_coord(0, 0, 0, HUGE_VAL)
+        cdef PJ_FACTORS pj_factors
+        cdef int errno = 0
+        cdef bint invalid_coord = 0
+        cdef Py_ssize_t iii
+        with nogil:
+            for iii in range(lonbuff.len):
+                pj_coord.uv.u = lonbuff.data[iii]
+                pj_coord.uv.v = latbuff.data[iii]
+                if not radians:
+                    pj_coord.uv.u *= _DG2RAD
+                    pj_coord.uv.v *= _DG2RAD
+
+                # set both to HUGE_VAL if inf or nan
+                proj_errno_reset(self.projobj)
+                if pj_coord.uv.v == HUGE_VAL \
+                        or pj_coord.uv.v != pj_coord.uv.v \
+                        or pj_coord.uv.u == HUGE_VAL \
+                        or pj_coord.uv.u != pj_coord.uv.u:
+                    invalid_coord = True
+                else:
+                    invalid_coord = False
+                    pj_factors = proj_factors(self.projobj, pj_coord)
+
+                errno = proj_errno(self.projobj)
+                if errcheck and errno:
+                    with gil:
+                        raise ProjError("proj error: {}".format(
+                            pystrdecode(proj_errno_string(errno))))
+
+                if errno or invalid_coord:
+                    meridional_scale_buff.data[iii] = HUGE_VAL
+                    parallel_scale_buff.data[iii] = HUGE_VAL
+                    areal_scale_buff.data[iii] = HUGE_VAL
+                    angular_distortion_buff.data[iii] = HUGE_VAL
+                    meridian_parallel_angle_buff.data[iii] = HUGE_VAL
+                    meridian_convergence_buff.data[iii] = HUGE_VAL
+                    tissot_semimajor_buff.data[iii] = HUGE_VAL
+                    tissot_semiminor_buff.data[iii] = HUGE_VAL
+                    dx_dlam_buff.data[iii] = HUGE_VAL
+                    dx_dphi_buff.data[iii] = HUGE_VAL
+                    dy_dlam_buff.data[iii] = HUGE_VAL
+                    dy_dphi_buff.data[iii] = HUGE_VAL
+                else:
+                    meridional_scale_buff.data[iii] = pj_factors.meridional_scale
+                    parallel_scale_buff.data[iii] = pj_factors.parallel_scale
+                    areal_scale_buff.data[iii] = pj_factors.areal_scale
+                    angular_distortion_buff.data[iii] = pj_factors.angular_distortion * _RAD2DG
+                    meridian_parallel_angle_buff.data[iii] = pj_factors.meridian_parallel_angle * _RAD2DG
+                    meridian_convergence_buff.data[iii] = pj_factors.meridian_convergence * _RAD2DG
+                    tissot_semimajor_buff.data[iii] = pj_factors.tissot_semimajor
+                    tissot_semiminor_buff.data[iii] = pj_factors.tissot_semiminor
+                    dx_dlam_buff.data[iii] = pj_factors.dx_dlam
+                    dx_dphi_buff.data[iii] = pj_factors.dx_dphi
+                    dy_dlam_buff.data[iii] = pj_factors.dy_dlam
+                    dy_dphi_buff.data[iii] = pj_factors.dy_dphi
+
+        ProjError.clear()
+
+        return Factors(
+            meridional_scale=meridional_scale,
+            parallel_scale=parallel_scale,
+            areal_scale=areal_scale,
+            angular_distortion=angular_distortion,
+            meridian_parallel_angle=meridian_parallel_angle,
+            meridian_convergence=meridian_convergence,
+            tissot_semimajor=tissot_semimajor,
+            tissot_semiminor=tissot_semiminor,
+            dx_dlam=dx_dlam,
+            dx_dphi=dx_dphi,
+            dy_dlam=dy_dlam,
+            dy_dphi=dy_dphi,
+        )
+
     def __repr__(self):
         return "Proj('{srs}', preserve_units=True)".format(srs=self.srs)
 
-    def _is_exact_same(self, Proj other):
+    def _is_exact_same(self, _Proj other):
         return proj_is_equivalent_to(
             self.projobj, other.projobj, PJ_COMP_STRICT) == 1
 
-    def _is_equivalent(self, Proj other):
+    def _is_equivalent(self, _Proj other):
         return proj_is_equivalent_to(
             self.projobj, other.projobj, PJ_COMP_EQUIVALENT) == 1
 
     def is_exact_same(self, other):
         """Compares Proj objects to see if they are exactly the same."""
-        if not isinstance(other, Proj):
+        if not isinstance(other, _Proj):
             return False
         return self._is_exact_same(other)

pyproj/proj.pxi

@@ -434,3 +434,19 @@ cdef extern from "proj.h":
         PROJ_GRID_AVAILABILITY_USED_FOR_SORTING
         PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID
         PROJ_GRID_AVAILABILITY_IGNORED
+
+    ctypedef struct PJ_FACTORS:
+        double meridional_scale
+        double parallel_scale
+        double areal_scale
+        double angular_distortion
+        double meridian_parallel_angle
+        double meridian_convergence
+        double tissot_semimajor
+        double tissot_semiminor
+        double dx_dlam
+        double dx_dphi
+        double dy_dlam
+        double dy_dphi
+
+    PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) nogil

pyproj/proj.py

@@ -37,19 +37,16 @@
 import re
 import warnings
 
-from pyproj import _proj
 from pyproj._list import get_proj_operations_map
+from pyproj._proj import Factors, _Proj, proj_version_str  # noqa: F401
 from pyproj.compat import cstrencode, pystrdecode
 from pyproj.crs import CRS
 from pyproj.utils import _convertback, _copytobuffer
 
-# import numpy as np
-proj_version_str = _proj.proj_version_str
-
 pj_list = get_proj_operations_map()
 
 
-class Proj(_proj.Proj):
+class Proj(_Proj):
     """
     Performs cartographic transformations (converts from
     longitude,latitude to native map projection x,y coordinates and
@@ -204,6 +201,73 @@ def __call__(self, *args, **kw):
         outy = _convertback(yisfloat, yislist, xistuple, iny)
         return outx, outy
 
+    def get_factors(
+        self, longitude, latitude, radians=False, errcheck=False,
+    ):
+        """
+        .. versionadded:: 2.6.0
+
+        Calculate various cartographic properties, such as scale factors, angular
+        distortion and meridian convergence. Depending on the underlying projection
+        values will be calculated either numerically (default) or analytically.
+
+        The function also calculates the partial derivatives of the given
+        coordinate.
+
+        Parameters
+        ----------
+        longitude: scalar or array (numpy or python)
+            Input longitude coordinate(s).
+        latitude: scalar or array (numpy or python)
+            Input latitude coordinate(s).
+        radians: boolean, optional
+            If True, will expect input data to be in radians.
+            Default is False (degrees).
+        errcheck: boolean, optional (default False)
+            If True an exception is raised if the errors are found in the process.
+            By default errcheck=False and ``inf`` is returned.
+
+        Returns
+        -------
+        Factors
+        """
+        # process inputs, making copies that support buffer API.
+        inx, xisfloat, xislist, xistuple = _copytobuffer(longitude)
+        iny, yisfloat, yislist, yistuple = _copytobuffer(latitude)
+
+        # calculate the factors
+        factors = self._get_factors(inx, iny, radians=radians, errcheck=errcheck)
+
+        # if inputs were lists, tuples or floats, convert back.
+        return Factors(
+            meridional_scale=_convertback(
+                xisfloat, xislist, xistuple, factors.meridional_scale
+            ),
+            parallel_scale=_convertback(
+                xisfloat, xislist, xistuple, factors.parallel_scale
+            ),
+            areal_scale=_convertback(xisfloat, xislist, xistuple, factors.areal_scale),
+            angular_distortion=_convertback(
+                xisfloat, xislist, xistuple, factors.angular_distortion
+            ),
+            meridian_parallel_angle=_convertback(
+                xisfloat, xislist, xistuple, factors.meridian_parallel_angle
+            ),
+            meridian_convergence=_convertback(
+                xisfloat, xislist, xistuple, factors.meridian_convergence
+            ),
+            tissot_semimajor=_convertback(
+                xisfloat, xislist, xistuple, factors.tissot_semimajor
+            ),
+            tissot_semiminor=_convertback(
+                xisfloat, xislist, xistuple, factors.tissot_semiminor
+            ),
+            dx_dlam=_convertback(xisfloat, xislist, xistuple, factors.dx_dlam),
+            dx_dphi=_convertback(xisfloat, xislist, xistuple, factors.dx_dphi),
+            dy_dlam=_convertback(xisfloat, xislist, xistuple, factors.dy_dlam),
+            dy_dphi=_convertback(xisfloat, xislist, xistuple, factors.dy_dphi),
+        )
+
     def definition_string(self):
         """Returns formal definition string for projection