1<?php 2 3namespace geoPHP\Geometry; 4 5use geoPHP\Exception\InvalidGeometryException; 6use geoPHP\geoPHP; 7 8/** 9 * A LineString is defined by a sequence of points, (X,Y) pairs, which define the reference points of the line string. 10 * Linear interpolation between the reference points defines the resulting linestring. 11 * 12 * @method Point[] getComponents() 13 * @property Point[] $components 14 * @method Point geometryN($n) 15 */ 16class LineString extends Curve 17{ 18 19 public function geometryType() 20 { 21 return Geometry::LINE_STRING; 22 } 23 24 /** 25 * Constructor 26 * 27 * @param Point[] $points An array of at least two points with 28 * which to build the LineString 29 * @throws \Exception 30 */ 31 public function __construct($points = []) 32 { 33 parent::__construct($points); 34 } 35 36 public static function fromArray($array) 37 { 38 $points = []; 39 foreach ($array as $point) { 40 $points[] = Point::fromArray($point); 41 } 42 return new static($points); 43 } 44 45 /** 46 * Returns the number of points of the LineString 47 * 48 * @return int 49 */ 50 public function numPoints() 51 { 52 return count($this->components); 53 } 54 55 /** 56 * Returns the 1-based Nth point of the LineString. 57 * Negative values are counted backwards from the end of the LineString. 58 * 59 * @param int $n Nth point of the LineString 60 * @return Point|null 61 */ 62 public function pointN($n) 63 { 64 return $n >= 0 65 ? $this->geometryN($n) 66 : $this->geometryN(count($this->components) - abs($n + 1)); 67 } 68 69 public function centroid() 70 { 71 return $this->getCentroidAndLength(); 72 } 73 74 public function getCentroidAndLength(&$length = 0.0) 75 { 76 if ($this->isEmpty()) { 77 return new Point(); 78 } 79 80 if ($this->getGeos()) { 81 // @codeCoverageIgnoreStart 82 /** @noinspection PhpUndefinedMethodInspection */ 83 return geoPHP::geosToGeometry($this->getGeos()->centroid()); 84 // @codeCoverageIgnoreEnd 85 } 86 87 $x = 0; 88 $y = 0; 89 $length = 0.0; 90 /** @var Point $previousPoint */ 91 $previousPoint = null; 92 foreach ($this->getPoints() as $point) { 93 if ($previousPoint) { 94 // Equivalent to $previousPoint->distance($point) but much faster 95 $segmentLength = sqrt( 96 pow(($previousPoint->x() - $point->x()), 2) + 97 pow(($previousPoint->y() - $point->y()), 2) 98 ); 99 $length += $segmentLength; 100 $x += ($previousPoint->x() + $point->x()) / 2 * $segmentLength; 101 $y += ($previousPoint->y() + $point->y()) / 2 * $segmentLength; 102 } 103 $previousPoint = $point; 104 } 105 if ($length === 0.0) { 106 return $this->startPoint(); 107 } 108 return new Point($x / $length, $y / $length); 109 } 110 111 /** 112 * Returns the length of this Curve in its associated spatial reference. 113 * Eg. if Geometry is in geographical coordinate system it returns the length in degrees 114 * @return float|int 115 */ 116 public function length() 117 { 118 if ($this->getGeos()) { 119 // @codeCoverageIgnoreStart 120 /** @noinspection PhpUndefinedMethodInspection */ 121 return $this->getGeos()->length(); 122 // @codeCoverageIgnoreEnd 123 } 124 $length = 0.0; 125 /** @var Point $previousPoint */ 126 $previousPoint = null; 127 foreach ($this->getPoints() as $point) { 128 if ($previousPoint) { 129 $length += sqrt( 130 pow(($previousPoint->x() - $point->x()), 2) + 131 pow(($previousPoint->y() - $point->y()), 2) 132 ); 133 } 134 $previousPoint = $point; 135 } 136 return $length; 137 } 138 139 public function length3D() 140 { 141 $length = 0.0; 142 /** @var Point $previousPoint */ 143 $previousPoint = null; 144 foreach ($this->getPoints() as $point) { 145 if ($previousPoint) { 146 $length += sqrt( 147 pow(($previousPoint->x() - $point->x()), 2) + 148 pow(($previousPoint->y() - $point->y()), 2) + 149 pow(($previousPoint->z() - $point->z()), 2) 150 ); 151 } 152 $previousPoint = $point; 153 } 154 return $length; 155 } 156 157 /** 158 * @param float|null $radius Earth radius 159 * @return float Length in meters 160 */ 161 public function greatCircleLength($radius = geoPHP::EARTH_WGS84_SEMI_MAJOR_AXIS) 162 { 163 $length = 0.0; 164 $rad = M_PI / 180; 165 $points = $this->getPoints(); 166 $numPoints = $this->numPoints() - 1; 167 for ($i = 0; $i < $numPoints; ++$i) { 168 // Simplified Vincenty formula with equal major and minor axes (a sphere) 169 $lat1 = $points[$i]->y() * $rad; 170 $lat2 = $points[$i + 1]->y() * $rad; 171 $lon1 = $points[$i]->x() * $rad; 172 $lon2 = $points[$i + 1]->x() * $rad; 173 $deltaLon = $lon2 - $lon1; 174 $d = 175 $radius * 176 atan2( 177 sqrt( 178 pow(cos($lat2) * sin($deltaLon), 2) + 179 pow(cos($lat1) * sin($lat2) - sin($lat1) * cos($lat2) * cos($deltaLon), 2) 180 ), 181 sin($lat1) * sin($lat2) + 182 cos($lat1) * cos($lat2) * cos($deltaLon) 183 ); 184 if ($points[$i]->is3D()) { 185 $d = sqrt( 186 pow($d, 2) + 187 pow($points[$i + 1]->z() - $points[$i]->z(), 2) 188 ); 189 } 190 191 $length += $d; 192 } 193 // Returns length in meters. 194 return $length; 195 } 196 197 /** 198 * @return float Haversine length of geometry in degrees 199 */ 200 public function haversineLength() 201 { 202 $distance = 0.0; 203 $points = $this->getPoints(); 204 $numPoints = $this->numPoints() - 1; 205 for ($i = 0; $i < $numPoints; ++$i) { 206 $point = $points[$i]; 207 $nextPoint = $points[$i + 1]; 208 $degree = (geoPHP::EARTH_WGS84_SEMI_MAJOR_AXIS * 209 acos( 210 sin(deg2rad($point->y())) * sin(deg2rad($nextPoint->y())) + 211 cos(deg2rad($point->y())) * cos(deg2rad($nextPoint->y())) * 212 cos(deg2rad(abs($point->x() - $nextPoint->x()))) 213 ) 214 ); 215 if (!is_nan($degree)) { 216 $distance += $degree; 217 } 218 } 219 return $distance; 220 } 221 222 /** 223 * @source https://github.com/mjaschen/phpgeo/blob/master/src/Location/Distance/Vincenty.php 224 * @author Marcus Jaschen <mjaschen@gmail.com> 225 * @license https://opensource.org/licenses/GPL-3.0 GPL 226 * (note: geoPHP uses "GPL version 2 (or later)" license which is compatible with GPLv3) 227 * 228 * @return float Length in meters 229 */ 230 public function vincentyLength() 231 { 232 $length = 0.0; 233 $rad = M_PI / 180; 234 $points = $this->getPoints(); 235 $numPoints = $this->numPoints() - 1; 236 for ($i = 0; $i < $numPoints; ++$i) { 237 // Inverse Vincenty formula 238 $lat1 = $points[$i]->y() * $rad; 239 $lat2 = $points[$i + 1]->y() * $rad; 240 $lng1 = $points[$i]->x() * $rad; 241 $lng2 = $points[$i + 1]->x() * $rad; 242 243 $a = geoPHP::EARTH_WGS84_SEMI_MAJOR_AXIS; 244 $b = geoPHP::EARTH_WGS84_SEMI_MINOR_AXIS; 245 $f = 1 / geoPHP::EARTH_WGS84_FLATTENING; 246 $L = $lng2 - $lng1; 247 $U1 = atan((1 - $f) * tan($lat1)); 248 $U2 = atan((1 - $f) * tan($lat2)); 249 $iterationLimit = 100; 250 $lambda = $L; 251 $sinU1 = sin($U1); 252 $sinU2 = sin($U2); 253 $cosU1 = cos($U1); 254 $cosU2 = cos($U2); 255 do { 256 $sinLambda = sin($lambda); 257 $cosLambda = cos($lambda); 258 $sinSigma = sqrt( 259 ($cosU2 * $sinLambda) * 260 ($cosU2 * $sinLambda) + 261 ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) * 262 ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) 263 ); 264 if ($sinSigma == 0) { 265 return 0.0; 266 } 267 $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda; 268 $sigma = atan2($sinSigma, $cosSigma); 269 $sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma; 270 $cosSqAlpha = 1 - $sinAlpha * $sinAlpha; 271 $cos2SigmaM = 0; 272 if ($cosSqAlpha <> 0) { 273 $cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha; 274 } 275 $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha)); 276 $lambdaP = $lambda; 277 $lambda = $L + (1 - $C) * $f * $sinAlpha * 278 ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (- 1 + 2 * $cos2SigmaM * $cos2SigmaM))); 279 } while (abs($lambda - $lambdaP) > 1e-12 && --$iterationLimit > 0); 280 if ($iterationLimit == 0) { 281 return null; // not converging 282 } 283 $uSq = $cosSqAlpha * ($a * $a - $b * $b) / ($b * $b); 284 $A = 1 + $uSq / 16384 * (4096 + $uSq * (- 768 + $uSq * (320 - 175 * $uSq))); 285 $B = $uSq / 1024 * (256 + $uSq * (- 128 + $uSq * (74 - 47 * $uSq))); 286 $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * 287 ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6 288 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) 289 * (-3 + 4 * $cos2SigmaM * $cos2SigmaM))); 290 291 $length += $b * $A * ($sigma - $deltaSigma); 292 } 293 // Returns length in meters. 294 return $length; 295 } 296 297 public function minimumZ() 298 { 299 $min = PHP_INT_MAX; 300 foreach ($this->getPoints() as $point) { 301 if ($point->hasZ() && $point->z() < $min) { 302 $min = $point->z(); 303 } 304 } 305 return $min < PHP_INT_MAX ? $min : null; 306 } 307 308 public function maximumZ() 309 { 310 $max = ~PHP_INT_MAX; 311 foreach ($this->getPoints() as $point) { 312 if ($point->hasZ() && $point->z() > $max) { 313 $max = $point->z(); 314 } 315 } 316 317 return $max > ~PHP_INT_MAX ? $max : null; 318 } 319 320 public function zDifference() 321 { 322 if ($this->startPoint()->hasZ() && $this->endPoint()->hasZ()) { 323 return abs($this->startPoint()->z() - $this->endPoint()->z()); 324 } else { 325 return null; 326 } 327 } 328 329 /** 330 * Returns the cumulative elevation gain of the LineString 331 * 332 * @param int|float|null $verticalTolerance Smoothing factor filtering noisy elevation data. 333 * Its unit equals to the z-coordinates unit (meters for geographical coordinates) 334 * If the elevation data comes from a DEM, a value around 3.5 can be acceptable. 335 * 336 * @return float 337 */ 338 public function elevationGain($verticalTolerance = 0) 339 { 340 $gain = 0.0; 341 $lastEle = $this->startPoint()->z(); 342 $pointCount = $this->numPoints(); 343 foreach ($this->getPoints() as $i => $point) { 344 if (abs($point->z() - $lastEle) > $verticalTolerance || $i === $pointCount - 1) { 345 if ($point->z() > $lastEle) { 346 $gain += $point->z() - $lastEle; 347 } 348 $lastEle = $point->z(); 349 } 350 } 351 return $gain; 352 } 353 354 /** 355 * Returns the cumulative elevation loss of the LineString 356 * 357 * @param int|float|null $verticalTolerance Smoothing factor filtering noisy elevation data. 358 * Its unit equals to the z-coordinates unit (meters for geographical coordinates) 359 * If the elevation data comes from a DEM, a value around 3.5 can be acceptable. 360 * 361 * @return float 362 */ 363 public function elevationLoss($verticalTolerance = 0) 364 { 365 $loss = 0.0; 366 $lastEle = $this->startPoint()->z(); 367 $pointCount = $this->numPoints(); 368 foreach ($this->getPoints() as $i => $point) { 369 if (abs($point->z() - $lastEle) > $verticalTolerance || $i === $pointCount - 1) { 370 if ($point->z() < $lastEle) { 371 $loss += $lastEle - $point->z(); 372 } 373 $lastEle = $point->z(); 374 } 375 } 376 return $loss; 377 } 378 379 public function minimumM() 380 { 381 $min = PHP_INT_MAX; 382 foreach ($this->getPoints() as $point) { 383 if ($point->isMeasured() && $point->m() < $min) { 384 $min = $point->m(); 385 } 386 } 387 return $min < PHP_INT_MAX ? $min : null; 388 } 389 390 public function maximumM() 391 { 392 $max = ~PHP_INT_MAX; 393 foreach ($this->getPoints() as $point) { 394 if ($point->isMeasured() && $point->m() > $max) { 395 $max = $point->m(); 396 } 397 } 398 399 return $max > ~PHP_INT_MAX ? $max : null; 400 } 401 402 /** 403 * Get all line segments 404 * @param bool $toArray return segments as LineString or array of start and end points 405 * 406 * @return LineString[]|array[Point] 407 */ 408 public function explode($toArray = false) 409 { 410 $points = $this->getPoints(); 411 $numPoints = count($points); 412 if ($numPoints < 2) { 413 return []; 414 } 415 $parts = []; 416 for ($i = 1; $i < $numPoints; ++$i) { 417 $segment = [$points[$i - 1], $points[$i]]; 418 $parts[] = $toArray ? $segment : new LineString($segment); 419 } 420 return $parts; 421 } 422 423 /** 424 * Checks that LineString is a Simple Geometry 425 * 426 * @return boolean 427 */ 428 public function isSimple() 429 { 430 if ($this->getGeos()) { 431 // @codeCoverageIgnoreStart 432 /** @noinspection PhpUndefinedMethodInspection */ 433 return $this->getGeos()->isSimple(); 434 // @codeCoverageIgnoreEnd 435 } 436 437 // As of OGR specification a ring is simple only if its start and end points equals in all coordinates 438 // Neither GEOS, nor PostGIS support it 439// if ($this->hasZ() 440// && $this->startPoint()->equals($this->endPoint()) 441// && $this->startPoint()->z() !== $this->endPoint()->z() 442// ) { 443// return false; 444// } 445 446 $segments = $this->explode(true); 447 foreach ($segments as $i => $segment) { 448 foreach ($segments as $j => $checkSegment) { 449 if ($i != $j) { 450 if (Geometry::segmentIntersects($segment[0], $segment[1], $checkSegment[0], $checkSegment[1])) { 451 return false; 452 } 453 } 454 } 455 } 456 return true; 457 } 458 459 /** 460 * @param LineString $segment 461 * @return bool 462 */ 463 public function lineSegmentIntersect($segment) 464 { 465 return Geometry::segmentIntersects( 466 $this->startPoint(), 467 $this->endPoint(), 468 $segment->startPoint(), 469 $segment->endPoint() 470 ); 471 } 472 473 /** 474 * @param Geometry|Collection $geometry 475 * @return float|null 476 */ 477 public function distance($geometry) 478 { 479 if ($this->getGeos()) { 480 // @codeCoverageIgnoreStart 481 /** @noinspection PhpUndefinedMethodInspection */ 482 return $this->getGeos()->distance($geometry->getGeos()); 483 // @codeCoverageIgnoreEnd 484 } 485 486 if ($geometry->geometryType() == Geometry::POINT) { 487 // This is defined in the Point class nicely 488 return $geometry->distance($this); 489 } 490 if ($geometry->geometryType() == Geometry::LINE_STRING) { 491 $distance = null; 492 $geometrySegments = $geometry->explode(); 493 foreach ($this->explode() as $seg1) { 494 /** @var LineString $seg2 */ 495 foreach ($geometrySegments as $seg2) { 496 if ($seg1->lineSegmentIntersect($seg2)) { 497 return 0.0; 498 } 499 // Because line-segments are straight, the shortest distance will occur at an endpoint. 500 // If they are parallel an endpoint calculation is still accurate. 501 $checkDistance1 = $seg1->startPoint()->distance($seg2); 502 $checkDistance2 = $seg1->endPoint()->distance($seg2); 503 $checkDistance3 = $seg2->startPoint()->distance($seg1); 504 $checkDistance4 = $seg2->endPoint()->distance($seg1); 505 506 $checkDistance = min($checkDistance1, $checkDistance2, $checkDistance3, $checkDistance4); 507 if ($checkDistance === 0.0) { 508 return 0.0; 509 } 510 if ($distance === null || $checkDistance < $distance) { 511 $distance = $checkDistance; 512 } 513 } 514 } 515 return $distance; 516 } else { 517 // It can be treated as collection 518 return parent::distance($geometry); 519 } 520 } 521} 522