diff --git a/README.md b/README.md index dbe50f94..38de23d2 100644 --- a/README.md +++ b/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. diff --git a/src/Polygon/MultiPolygon.php b/src/Polygon/MultiPolygon.php index 8f5ecfe8..a5cbe620 100644 --- a/src/Polygon/MultiPolygon.php +++ b/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; + } } diff --git a/src/Polygon/Polygon.php b/src/Polygon/Polygon.php index 4f024afa..cf17915d 100644 --- a/src/Polygon/Polygon.php +++ b/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 */ diff --git a/src/Polygon/PolygonCenterTrait.php b/src/Polygon/PolygonCenterTrait.php new file mode 100644 index 00000000..b1f62f16 --- /dev/null +++ b/src/Polygon/PolygonCenterTrait.php @@ -0,0 +1,60 @@ + + * + * 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)); + } +} diff --git a/src/Polygon/PolygonInterface.php b/src/Polygon/PolygonInterface.php index 61eaa312..0ef76599 100644 --- a/src/Polygon/PolygonInterface.php +++ b/src/Polygon/PolygonInterface.php @@ -36,4 +36,20 @@ public function pointOnBoundary(CoordinateInterface $coordinate); * @return boolean */ public function pointOnVertex(CoordinateInterface $coordinate); + + /** + * Returns the geographic centroid of the polygon's coordinates, + * calculated by averaging x/y/z unit-sphere vectors. + * + * @return CoordinateInterface|null null when the polygon is empty + */ + public function getCenter(): ?CoordinateInterface; + + /** + * 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; } diff --git a/tests/Polygon/PolygonTest.php b/tests/Polygon/PolygonTest.php index 421c19d4..1e705ad0 100644 --- a/tests/Polygon/PolygonTest.php +++ b/tests/Polygon/PolygonTest.php @@ -225,4 +225,111 @@ public function testPointNotInPolygon( $this->assertFalse($this->polygon->pointInPolygon(new Coordinate($pointNotInPolygonCoordinate))); } } + + public function testGetCenterReturnsNullOnEmptyPolygon() + { + $this->assertNull($this->polygon->getCenter()); + } + + public function testGetRadiusReturnsZeroOnEmptyPolygon() + { + $this->assertSame(0.0, $this->polygon->getRadius()); + } + + public function testGetCenterSingleCoordinate() + { + $this->polygon->add(new Coordinate([48.8566, 2.3522])); + $center = $this->polygon->getCenter(); + + $this->assertNotNull($center); + $this->assertEqualsWithDelta(48.8566, (float) $center->getLatitude(), 0.0001); + $this->assertEqualsWithDelta(2.3522, (float) $center->getLongitude(), 0.0001); + } + + public function testGetCenterSymmetricSquare() + { + // A square centred at (0, 0) — centre must be (0, 0). + $this->polygon->set([ + [1.0, -1.0], + [1.0, 1.0], + [-1.0, 1.0], + [-1.0, -1.0], + ]); + + $center = $this->polygon->getCenter(); + $this->assertNotNull($center); + $this->assertEqualsWithDelta(0.0, (float) $center->getLatitude(), 0.0001); + $this->assertEqualsWithDelta(0.0, (float) $center->getLongitude(), 0.0001); + } + + public function testGetCenterWithKnownPolygon() + { + // The centroid of these four points should be approximately at their + // geographic mean (verified independently with spherical averaging). + $this->polygon->set([ + [48.9675969, 1.7440796], + [48.4711003, 2.5268555], + [48.9279131, 3.1448364], + [49.3895245, 2.6119995], + ]); + + $center = $this->polygon->getCenter(); + $this->assertNotNull($center); + // The centroid latitude/longitude must fall inside the bounding box. + $this->assertGreaterThan(48.4, (float) $center->getLatitude()); + $this->assertLessThan(49.4, (float) $center->getLatitude()); + $this->assertGreaterThan(1.7, (float) $center->getLongitude()); + $this->assertLessThan(3.2, (float) $center->getLongitude()); + } + + public function testGetRadiusGreaterThanZeroForNonTrivialPolygon() + { + $this->polygon->set([ + [48.9675969, 1.7440796], + [48.4711003, 2.5268555], + [48.9279131, 3.1448364], + [49.3895245, 2.6119995], + ]); + + $radius = $this->polygon->getRadius(); + $this->assertGreaterThan(0.0, $radius); + // Rough sanity check: radius should be in tens-of-km range for this polygon. + $this->assertGreaterThan(50000.0, $radius); // > 50 km + $this->assertLessThan(200000.0, $radius); // < 200 km + } + + public function testGetCenterUpdatesIncrementallyOnAdd() + { + $coord1 = new Coordinate([10.0, 20.0]); + $coord2 = new Coordinate([10.0, 22.0]); + + $this->polygon->add($coord1); + $center1 = $this->polygon->getCenter(); + $this->assertNotNull($center1); + $this->assertEqualsWithDelta(10.0, (float) $center1->getLatitude(), 0.0001); + $this->assertEqualsWithDelta(20.0, (float) $center1->getLongitude(), 0.0001); + + $this->polygon->add($coord2); + $center2 = $this->polygon->getCenter(); + $this->assertNotNull($center2); + // After adding the second point, the centroid longitude should be ~21. + $this->assertEqualsWithDelta(10.0, (float) $center2->getLatitude(), 0.01); + $this->assertEqualsWithDelta(21.0, (float) $center2->getLongitude(), 0.01); + } + + public function testGetCenterUpdatesAfterRemove() + { + $this->polygon->set([ + [10.0, 20.0], + [10.0, 22.0], + [10.0, 24.0], + ]); + + // Remove the last coordinate; centre should shift towards [10, 21]. + $this->polygon->remove(2); + $center = $this->polygon->getCenter(); + $this->assertNotNull($center); + $this->assertEqualsWithDelta(10.0, (float) $center->getLatitude(), 0.01); + $this->assertEqualsWithDelta(21.0, (float) $center->getLongitude(), 0.01); + } }