Add getCenter() and getRadius() to Polygon

README.md

@@ -473,8 +473,19 @@ $polygon->pointOnBoundary(new \League\Geotools\Coordinate\Coordinate([47.1587188
 $polygon->pointOnVertex(new \League\Geotools\Coordinate\Coordinate([48.4711003, 2.5268555])); // true
 $polygon->pointOnVertex(new \League\Geotools\Coordinate\Coordinate([49.1785607, 2.4444580])); // false
 $polygon->getBoundingBox(); // return the BoundingBox object
+$polygon->getCenter(); // return a Coordinate representing the geographic centroid of the polygon
+$polygon->getRadius(); // return the radius in meters (max haversine distance from centroid to any vertex)
 ```
 
+`getCenter()` computes the centroid by averaging the unit-sphere (x, y, z) vectors of all
+vertices and projecting the result back to (latitude, longitude).  Running sums of x/y/z are
+maintained incrementally as coordinates are added or removed via `add()` and `remove()`, so
+repeated calls after single-step mutations are O(1).
+
+`getRadius()` returns the maximum haversine distance (in metres) between the centroid and any
+vertex, which can be used to detect when a sliding-window polygon has grown beyond a desired
+size.
+
 ## CLI
 
 It provides command lines to compute methods provided by **Distance**, **Point**, **Geohash** and **Convert** classes.

src/Polygon/MultiPolygon.php

@@ -2,11 +2,13 @@
 
 namespace League\Geotools\Polygon;
 
+use League\Geotools\Coordinate\Coordinate;
 use League\Geotools\Coordinate\CoordinateInterface;
 use League\Geotools\GeometryCollection;
 
 class MultiPolygon extends GeometryCollection implements PolygonInterface
 {
+    use PolygonCenterTrait;
     const TYPE = 'MULTIPOLYGON';
 
     /**
@@ -64,4 +66,77 @@ public function pointOnVertex(CoordinateInterface $coordinate)
 
         return false;
     }
+
+    /**
+     * Returns the geographic centroid of all coordinates across every polygon
+     * in this collection.
+     *
+     * @return CoordinateInterface|null  null when the collection is empty
+     */
+    public function getCenter(): ?CoordinateInterface
+    {
+        $sumX = $sumY = $sumZ = 0.0;
+        $count = 0;
+        $ellipsoid = null;
+
+        foreach ($this->getCoordinates() as $coordinate) {
+            [$x, $y, $z] = $this->coordinateToXYZ($coordinate);
+            $sumX += $x;
+            $sumY += $y;
+            $sumZ += $z;
+            $count++;
+            if ($ellipsoid === null) {
+                $ellipsoid = $coordinate->getEllipsoid();
+            }
+        }
+
+        if ($count === 0) {
+            return null;
+        }
+
+        $x = $sumX / $count;
+        $y = $sumY / $count;
+        $z = $sumZ / $count;
+
+        $lng = rad2deg(atan2($y, $x));
+        $hyp = sqrt($x * $x + $y * $y);
+        $lat = rad2deg(atan2($z, $hyp));
+
+        return new Coordinate([$lat, $lng], $ellipsoid);
+    }
+
+    /**
+     * Returns the radius of the multi-polygon in meters, defined as the
+     * maximum haversine distance from the centroid to any vertex.
+     *
+     * @return float
+     */
+    public function getRadius(): float
+    {
+        $center = $this->getCenter();
+        if ($center === null) {
+            return 0.0;
+        }
+
+        $radius = 0.0;
+        $R    = $center->getEllipsoid()->getArithmeticMeanRadius();
+        $latC = deg2rad($center->getLatitude());
+        $lngC = deg2rad($center->getLongitude());
+
+        foreach ($this->getCoordinates() as $coordinate) {
+            $d = $this->haversineDistance(
+                $R,
+                $latC,
+                $lngC,
+                deg2rad($coordinate->getLatitude()),
+                deg2rad($coordinate->getLongitude())
+            );
+
+            if ($d > $radius) {
+                $radius = $d;
+            }
+        }
+
+        return $radius;
+    }
 }

src/Polygon/Polygon.php

@@ -23,6 +23,7 @@
  */
 class Polygon implements PolygonInterface, \Countable, \IteratorAggregate, \ArrayAccess, \JsonSerializable
 {
+    use PolygonCenterTrait;
     const TYPE = 'POLYGON';
 
     /**
@@ -45,6 +46,27 @@ class Polygon implements PolygonInterface, \Countable, \IteratorAggregate, \Arra
      */
     private $precision = 8;
 
+    /**
+     * Running sum of the x-components of coordinates on the unit sphere.
+     *
+     * @var float
+     */
+    private $sumX = 0.0;
+
+    /**
+     * Running sum of the y-components of coordinates on the unit sphere.
+     *
+     * @var float
+     */
+    private $sumY = 0.0;
+
+    /**
+     * Running sum of the z-components of coordinates on the unit sphere.
+     *
+     * @var float
+     */
+    private $sumZ = 0.0;
+
     /**
      * @param null|array|CoordinateCollection $coordinates
      */
@@ -55,6 +77,7 @@ public function __construct($coordinates = null)
         } elseif ($coordinates instanceof CoordinateCollection) {
             $this->coordinates = $coordinates;
             $this->hasCoordinate = $coordinates->count() > 0;
+            $this->recalculateSums();
         } else {
             throw new \InvalidArgumentException;
         }
@@ -294,6 +317,7 @@ public function setCoordinates(CoordinateCollection $coordinates)
     {
         $this->coordinates = $coordinates;
         $this->boundingBox->setPolygon($this);
+        $this->recalculateSums();
 
         return $this;
     }
@@ -339,7 +363,12 @@ public function offsetGet($offset)
     #[\ReturnTypeWillChange]
     public function offsetSet($offset, $value)
     {
+        $existing = $this->coordinates->get($offset);
+        if ($existing !== null) {
+            $this->subtractFromSums($existing);
+        }
         $this->coordinates->offsetSet($offset, $value);
+        $this->addToSums($value);
         $this->boundingBox->setPolygon($this);
     }
 
@@ -349,6 +378,10 @@ public function offsetSet($offset, $value)
     #[\ReturnTypeWillChange]
     public function offsetUnset($offset)
     {
+        $existing = $this->coordinates->get($offset);
+        if ($existing !== null) {
+            $this->subtractFromSums($existing);
+        }
         $retval = $this->coordinates->offsetUnset($offset);
         $this->boundingBox->setPolygon($this);
         return $retval;
@@ -402,6 +435,7 @@ public function set($key, ?CoordinateInterface $coordinate = null)
 
         $this->hasCoordinate = true;
         $this->boundingBox->setPolygon($this);
+        $this->recalculateSums();
     }
 
     /**
@@ -413,6 +447,7 @@ public function add(CoordinateInterface $coordinate)
 
         $this->hasCoordinate = true;
         $this->boundingBox->setPolygon($this);
+        $this->addToSums($coordinate);
 
         return $retval;
     }
@@ -422,16 +457,120 @@ public function add(CoordinateInterface $coordinate)
      */
     public function remove($key)
     {
+        $coordinate = $this->coordinates->get($key);
         $retval = $this->coordinates->remove($key);
 
         if (!count($this->coordinates)) {
             $this->hasCoordinate = false;
+            $this->sumX = $this->sumY = $this->sumZ = 0.0;
+        } elseif ($coordinate !== null) {
+            $this->subtractFromSums($coordinate);
         }
         $this->boundingBox->setPolygon($this);
 
         return $retval;
     }
 
+    /**
+     * Returns the geographic centroid of the polygon's coordinates.
+     *
+     * The centroid is computed by averaging the unit-sphere (x, y, z) vectors
+     * of all vertices and projecting the result back to a (latitude, longitude)
+     * coordinate.  Running sums of x/y/z are maintained incrementally so that
+     * repeated calls are O(1) after the initial population.
+     *
+     * @return CoordinateInterface|null  null when the polygon is empty
+     */
+    public function getCenter(): ?CoordinateInterface
+    {
+        $count = $this->count();
+        if ($count === 0) {
+            return null;
+        }
+
+        $x = $this->sumX / $count;
+        $y = $this->sumY / $count;
+        $z = $this->sumZ / $count;
+
+        $lng = rad2deg(atan2($y, $x));
+        $hyp = sqrt($x * $x + $y * $y);
+        $lat = rad2deg(atan2($z, $hyp));
+
+        return new Coordinate([$lat, $lng], $this->getEllipsoid());
+    }
+
+    /**
+     * Returns the radius of the polygon in meters, defined as the maximum
+     * haversine distance from the centroid to any vertex.
+     *
+     * @return float
+     */
+    public function getRadius(): float
+    {
+        $center = $this->getCenter();
+        if ($center === null) {
+            return 0.0;
+        }
+
+        $radius = 0.0;
+        $R = $center->getEllipsoid()->getArithmeticMeanRadius();
+        $latC = deg2rad($center->getLatitude());
+        $lngC = deg2rad($center->getLongitude());
+
+        foreach ($this->coordinates as $coordinate) {
+            $d = $this->haversineDistance(
+                $R,
+                $latC,
+                $lngC,
+                deg2rad($coordinate->getLatitude()),
+                deg2rad($coordinate->getLongitude())
+            );
+
+            if ($d > $radius) {
+                $radius = $d;
+            }
+        }
+
+        return $radius;
+    }
+
+    /**
+     * Adds a coordinate's unit-sphere contribution to the running sums.
+     *
+     * @param CoordinateInterface $coordinate
+     */
+    private function addToSums(CoordinateInterface $coordinate): void
+    {
+        [$x, $y, $z] = $this->coordinateToXYZ($coordinate);
+        $this->sumX += $x;
+        $this->sumY += $y;
+        $this->sumZ += $z;
+    }
+
+    /**
+     * Subtracts a coordinate's unit-sphere contribution from the running sums.
+     *
+     * @param CoordinateInterface $coordinate
+     */
+    private function subtractFromSums(CoordinateInterface $coordinate): void
+    {
+        [$x, $y, $z] = $this->coordinateToXYZ($coordinate);
+        $this->sumX -= $x;
+        $this->sumY -= $y;
+        $this->sumZ -= $z;
+    }
+
+    /**
+     * Recalculates the running sums from scratch based on current coordinates.
+     */
+    private function recalculateSums(): void
+    {
+        $this->sumX = $this->sumY = $this->sumZ = 0.0;
+        foreach ($this->coordinates as $coordinate) {
+            $this->addToSums($coordinate);
+        }
+    }
+
     /**
      * @return integer
      */

src/Polygon/PolygonCenterTrait.php

@@ -0,0 +1,60 @@
+<?php
+
+/*
+ * This file is part of the Geotools library.
+ *
+ * (c) Antoine Corcy <contact@sbin.dk>
+ *
+ * For the full copyright and license information, please view the LICENSE
+ * file that was distributed with this source code.
+ */
+
+namespace League\Geotools\Polygon;
+
+use League\Geotools\Coordinate\CoordinateInterface;
+
+/**
+ * Shared helpers for computing the geographic centroid and radius of a
+ * collection of coordinates.
+ */
+trait PolygonCenterTrait
+{
+    /**
+     * Converts a lat/lng coordinate to its (x, y, z) unit-sphere representation.
+     *
+     * @param  CoordinateInterface $coordinate
+     * @return float[]  [x, y, z]
+     */
+    private function coordinateToXYZ(CoordinateInterface $coordinate): array
+    {
+        $lat = deg2rad($coordinate->getLatitude());
+        $lng = deg2rad($coordinate->getLongitude());
+
+        return [
+            cos($lat) * cos($lng),
+            cos($lat) * sin($lng),
+            sin($lat),
+        ];
+    }
+
+    /**
+     * Returns the haversine distance in metres between a centre coordinate and
+     * a vertex coordinate.
+     *
+     * @param  float $R     Mean radius of the earth in metres
+     * @param  float $latC  Centre latitude in radians
+     * @param  float $lngC  Centre longitude in radians
+     * @param  float $latV  Vertex latitude in radians
+     * @param  float $lngV  Vertex longitude in radians
+     * @return float        Distance in metres
+     */
+    private function haversineDistance(float $R, float $latC, float $lngC, float $latV, float $lngV): float
+    {
+        $dLat = $latV - $latC;
+        $dLng = $lngV - $lngC;
+
+        $a = (sin($dLat / 2) ** 2) + cos($latC) * cos($latV) * (sin($dLng / 2) ** 2);
+
+        return 2 * $R * asin(sqrt($a));
+    }
+}