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