Geodesic: Added polygon area computations & support for shapely geome…

…tries

Makefile

@@ -63,7 +63,7 @@ lint: ## check style with flake8
 
 check: lint ## flake8 black isort check
 	black --check setup.py pyproj/ test/ docs/
-	isort --recursive -m 3 -w 88 -tc -p test setup.py pyproj/ test/ docs/
+	isort --check --recursive -m 3 -w 88 -tc -p test setup.py pyproj/ test/ docs/
 
 isort: ## order imports
 	isort --recursive -m 3 -w 88 -tc -p test setup.py pyproj/ test/ docs/

docs/examples.rst

@@ -356,28 +356,89 @@ Step 2: Create Transformer to convert from geodetic CRS to CRS
 Geodesic calculations
 ---------------------
 This is useful if you need to calculate the distance between two
-points on Earth's surface.
+points or the area of a geometry on Earth's surface.
 
 For more examples of usage and documentation, see :class:`~pyproj.Geod`.
 
+
+Creating Geod class
+~~~~~~~~~~~~~~~~~~~
+
 This example demonstrates creating a :class:`~pyproj.Geod` using an
 ellipsoid name as well as deriving one using a :class:`~pyproj.crs.CRS`.
 
 .. code:: python
 
-    >>> from pyproj import CRS, Geod
     >>> from pyproj import CRS, Geod
     >>> geod_clrk = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid.
     >>> geod_clrk
     Geod(ellps='clrk66')
     >>> geod_wgs84 = CRS("epsg:4326").get_geod()
     >>> geod_wgs84
     Geod('+a=6378137 +f=0.0033528106647475126')
-    >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
-    >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
-    >>> az12, az21, dist = geod_clrk.inv(boston_lon, boston_lat, portland_lon, portland_lat)
-    >>> az12, az21, dist
-    (-66.5305947876623, 75.65363415556968, 4164192.7080994667)
-    >>> az12_wgs, az21_wgs, dist_wgs = geod_wgs84.inv(boston_lon, boston_lat, portland_lon, portland_lat)
-    >>> az12_wgs, az21_wgs, dist_wgs
-    (-66.53043696156747, 75.65384304856798, 4164074.2392955828)
+
+
+Geodesic line length
+~~~~~~~~~~~~~~~~~~~~
+
+Calculate the geodesic length of a line (See: :meth:`~pyproj.Geod.line_length`):
+
+.. code:: python
+
+    >>> from pyproj import Geod
+    >>> lats = [-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
+    ...         -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7]
+    >>> lons = [-74, -102, -102, -131, -163, 163, 172, 140, 113,
+    ...         88, 59, 25, -4, -14, -33, -46, -61]
+    >>> geod = Geod(ellps="WGS84")
+    >>> total_length = geod.line_length(lons, lats)
+    >>> "{:.3f}".format(total_length)
+    '14259605.611'
+
+Calculate the geodesic length of a shapely geometry (See: :meth:`~pyproj.Geod.geometry_length`):
+
+.. code:: python
+
+    >>> from pyproj import Geod
+    >>> from shapely.geometry import Point, LineString
+    >>> line_string = LineString([Point(1, 2), Point(3, 4)]))
+    >>> geod = Geod(ellps="WGS84")
+    >>> total_length = geod.geometry_length(line_string)
+    >>> "{:.3f}".format(total_length)
+    '313588.397'
+
+
+Geodesic area
+~~~~~~~~~~~~~
+
+Calculate the geodesic area and perimeter of a polygon (See: :meth:`~pyproj.Geod.polygon_area_perimeter`):
+
+.. code:: python
+
+    >>> from pyproj import Geod
+    >>> geod = Geod('+a=6378137 +f=0.0033528106647475126')
+    >>> lats = [-72.9, -71.9, -74.9, -74.3, -77.5, -77.4, -71.7, -65.9, -65.7,
+    ...         -66.6, -66.9, -69.8, -70.0, -71.0, -77.3, -77.9, -74.7]
+    >>> lons = [-74, -102, -102, -131, -163, 163, 172, 140, 113,
+    ...         88, 59, 25, -4, -14, -33, -46, -61]
+    >>> poly_area, poly_perimeter = geod.polygon_area_perimeter(lons, lats)
+    >>> "{:.3f} {:.3f}".format(poly_area, poly_perimeter)
+    '13376856682207.406 14710425.407'
+
+
+Calculate the geodesic area and perimeter of a shapely polygon (See: :meth:`~pyproj.Geod.geometry_area_perimeter`):
+
+
+.. code:: python
+
+    >>> from pyproj import Geod
+    >>> from shapely.geometry import LineString, Point, Polygon
+    >>> geod = Geod('+a=6378137 +f=0.0033528106647475126')
+    >>> poly_area, poly_perimeter = geod.geometry_area_perimeter(
+            Polygon(
+                LineString([Point(1, 1), Point(1, 10), Point(10, 10), Point(10, 1)]),
+                holes=[LineString([Point(1, 2), Point(3, 4), Point(5, 2)])],
+            )
+        )
+    >>> "{:.3f} {:.3f}".format(poly_area, poly_perimeter)
+    '-944373881400.339 3979008.036'    

docs/history.rst

@@ -1,6 +1,16 @@
 Change Log
 ==========
 
+2.3.0
+~~~~~
+* Added support for calculating geodesic area (:meth:`~pyproj.geod.Geod.polygon_area_perimeter`)
+  and added interface to calculate total length of a line 
+  (:meth:`~pyproj.geod.Geod.line_length` & :meth:`~pyproj.geod.Geod.line_lengths`) (issue #210).
+* Added support for calculating geodesic area and line lemgths with shapely geometries
+  (:meth:`~pyproj.geod.Geod.geometry_area_perimeter` & :meth:`~pyproj.geod.Geod.geometry_length`)
+  (pull #366)
+
+
 2.2.2
 ~~~~~
 * Add deprecation warning when using +init= syntax (pull #358)

pyproj/_geod.pxd

@@ -15,6 +15,10 @@ cdef extern from "geodesic.h":
                double lat1, double lon1, double azi1, unsigned caps)
   void geod_position(geod_geodesicline* l, double s12,
                double* plat2, double* plon2, double* pazi2);
+  void geod_polygonarea(geod_geodesic* g,
+                        double lats[], double lons[], int n,
+                        double* pA, double* pP);
+
   cdef enum:
       GEODESIC_VERSION_MAJOR
       GEODESIC_VERSION_MINOR

pyproj/_geod.pyx

@@ -28,12 +28,12 @@ cdef class Geod:
         """special method that allows pyproj.Geod instance to be pickled"""
         return self.__class__,(self.initstring,)
 
-    def _fwd(self, object lons, object lats, object az, object dist, radians=False):
+    def _fwd(self, object lons, object lats, object az, object dist, bint radians=False):
         """
- forward transformation - determine longitude, latitude and back azimuth
- of a terminus point given an initial point longitude and latitude, plus
- forward azimuth and distance.
- if radians=True, lons/lats are radians instead of degrees.
+        forward transformation - determine longitude, latitude and back azimuth
+        of a terminus point given an initial point longitude and latitude, plus
+        forward azimuth and distance.
+        if radians=True, lons/lats are radians instead of degrees.
         """
         cdef Py_ssize_t buflenlons, buflenlats, buflenaz, buflend, ndim, i
         cdef double lat1,lon1,az1,s12,plon2,plat2,pazi2
@@ -90,11 +90,11 @@ cdef class Geod:
                 latsdata[i] = _DG2RAD*plat2
                 azdata[i] = _DG2RAD*pazi2
 
-    def _inv(self, object lons1, object lats1, object lons2, object lats2, radians=False):
+    def _inv(self, object lons1, object lats1, object lons2, object lats2, bint radians=False):
         """
- inverse transformation - return forward and back azimuths, plus distance
- between an initial and terminus lat/lon pair.
- if radians=True, lons/lats are radians instead of degrees.
+        inverse transformation - return forward and back azimuths, plus distance
+        between an initial and terminus lat/lon pair.
+        if radians=True, lons/lats are radians instead of degree
         """
         cdef double lat1,lon1,lat2,lon2,pazi1,pazi2,ps12
         cdef Py_ssize_t buflenlons, buflenlats, buflenaz, buflend, ndim, i
@@ -150,9 +150,10 @@ cdef class Geod:
                 latsdata[i] = pazi2
             azdata[i] = ps12
 
-    def _npts(self, double lon1, double lat1, double lon2, double lat2, int npts, radians=False):
+    def _npts(self, double lon1, double lat1, double lon2, double lat2, int npts, bint radians=False):
+        """
+        given initial and terminus lat/lon, find npts intermediate points.
         """
- given initial and terminus lat/lon, find npts intermediate points."""
         cdef int i
         cdef double del_s,ps12,pazi1,pazi2,s12,plon2,plat2
         cdef geod_geodesicline line
@@ -184,6 +185,122 @@ cdef class Geod:
                 lons = lons + (plon2,)
         return lons, lats
 
+    def _line_length(self, object lons, object lats, bint radians=False):
+        """
+        Calculate the distance between points along a line.
+
+
+        Parameters
+        ----------
+        lons: array
+            The longitude points along a line.
+        lats: array
+            The latitude points along a line.
+        radians: bool, optional
+            If True, the input data is assumed to be in radians.
+
+        Returns
+        -------
+        float: The total distance.
+    
+        """
+        cdef double lat1,lon1,lat2,lon2,pazi1,pazi2,ps12
+        cdef double total_distance = 0.0
+        cdef Py_ssize_t buflenlons, buflenlats, ndim, iii
+        cdef double *lonsdata
+        cdef double *latsdata
+        cdef void *londata
+        cdef void *latdata
+        # if buffer api is supported, get pointer to data buffers.
+        if PyObject_AsWriteBuffer(lons, &londata, &buflenlons) <> 0:
+            raise GeodError
+        if PyObject_AsWriteBuffer(lats, &latdata, &buflenlats) <> 0:
+            raise GeodError
+        # process data in buffer
+        if buflenlons != buflenlats:
+            raise GeodError("Buffer lengths not the same")
+        ndim = buflenlons//_DOUBLESIZE
+        lonsdata = <double *>londata
+        latsdata = <double *>latdata
+
+        if ndim == 1:
+            lonsdata[0] = 0
+            return 0.0
+
+        for iii from 0 <= iii < ndim - 1:
+            if radians:
+                lon1 = _RAD2DG*lonsdata[iii]
+                lat1 = _RAD2DG*latsdata[iii]
+                lon2 = _RAD2DG*lonsdata[iii + 1]
+                lat2 = _RAD2DG*latsdata[iii + 1]
+            else:
+                lon1 = lonsdata[iii]
+                lat1 = latsdata[iii]
+                lon2 = lonsdata[iii + 1]
+                lat2 = latsdata[iii + 1]
+            geod_inverse(&self._geod_geodesic, lat1, lon1, lat2, lon2,
+                         &ps12, &pazi1, &pazi2)
+            lonsdata[iii] = ps12
+            total_distance += ps12
+        return total_distance
+
+    def _polygon_area_perimeter(self, object lons, object lats, bint radians=False):
+        """
+        A simple interface for computing the area of a geodesic polygon.
+        
+        lats should be in the range [-90 deg, 90 deg].
+        
+        Only simple polygons (which are not self-intersecting) are allowed.
+        There's no need to "close" the polygon by repeating the first vertex.
+        The area returned is signed with counter-clockwise traversal being treated as
+        positive.
+
+        Parameters
+        ----------
+        lons: array
+            An array of longitude values.
+        lats: array
+            An array of latitude values.
+        radians: bool, optional
+            If True, the input data is assumed to be in radians.
+
+        Returns
+        -------
+        (float, float): The area (meter^2) and permimeter (meters) of the polygon.
+
+        """
+        cdef Py_ssize_t buflenlons, buflenlats, ndim, iii
+        cdef void *londata
+        cdef void *latdata
+        cdef double *lonsdata
+        cdef double *latsdata
+        # if buffer api is supported, get pointer to data buffers.
+        if PyObject_AsWriteBuffer(lons, &londata, &buflenlons) <> 0:
+            raise GeodError
+        if PyObject_AsWriteBuffer(lats, &latdata, &buflenlats) <> 0:
+            raise GeodError
+        # process data in buffer
+        if not buflenlons == buflenlats:
+            raise GeodError("Buffer lengths not the same")
+
+        cdef double polygon_area
+        cdef double polygon_perimeter
+        ndim = buflenlons//_DOUBLESIZE
+
+        lonsdata = <double *>londata
+        latsdata = <double *>latdata
+        if radians:
+            for iii from 0 <= iii < ndim:
+                lonsdata[iii] *= _RAD2DG
+                latsdata[iii] *= _RAD2DG
+
+        geod_polygonarea(
+            &self._geod_geodesic, latsdata, lonsdata, ndim, 
+            &polygon_area, &polygon_perimeter
+        )
+        return (polygon_area, polygon_perimeter)
+
+
     def __repr__(self):
         return "{classname}({init!r})".format(classname=self.__class__.__name__,
                                               init=self.initstring)

pyproj/crs.py

@@ -44,9 +44,9 @@
     INVERSE_PROJ_PARAM_MAP,
     K_0_MAP,
     LON_0_MAP,
-    PROJ_PARAM_MAP,
-    PARAM_TO_CF_MAP,
     METHOD_NAME_TO_CF_MAP,
+    PARAM_TO_CF_MAP,
+    PROJ_PARAM_MAP,
 )
 from pyproj.compat import string_types
 from pyproj.exceptions import CRSError