Passed
Push — master ( 34d579...0aeeb2 )
by Swen
03:06
created

LineString::getCentroid()   B

Complexity

Conditions 6
Paths 7

Size

Total Lines 36
Code Lines 23

Duplication

Lines 0
Ratio 0 %

Importance

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