LineString::haversineLength() - Code Metrics - Inspection of "removes some small bugs" - swen100/geoPHP - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( 0fd618...5e6b25 )

by Swen

created 2021-12-06 07:42 UTC

LineString::haversineLength() A

↳ Parent: LineString

Complexity

Conditions	3
Paths	3

Size

Total Lines	20
Code Lines	14

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	14
c	0
b	0
f	0
dl	0
loc	20
rs	9.7998
cc	3
nc	3
nop	0

<?php
namespace geoPHP\Geometry;

use geoPHP\geoPHP;

/**
 * A LineString is defined by a sequence of points, (X,Y) pairs, which define the reference points of the line string.
 * Linear interpolation between the reference points defines the resulting linestring.
 *
 * @method   Point[] getComponents()
 * @property Point[] $components
 */
class LineString extends Curve
{

    /**
     * @return string "LineString"
     */
    public function geometryType(): string
    {
        return Geometry::LINESTRING;
    }

    /**
     * @param  array<array> $array
     * @return LineString
     */
    public static function fromArray(array $array): LineString
    {
        $points = [];
        foreach ($array as $point) {
            $points[] = Point::fromArray($point);
        }
        return new LineString($points);
    }

    /**
     * Returns the number of points of the LineString
     *
     * @return int
     */
    public function numPoints(): int
    {
        return count($this->components);
    }

    /**
     * Returns the 1-based Nth point of the LineString.
     * Negative values are counted backwards from the end of the LineString.
     *
     * @param  int $n Nth point of the LineString
     * @return Point|null
     */
    public function pointN(int $n)
    {
        $n = $n < $this->numPoints() ? $n : $this->numPoints();
        
        /** @var Point|null $point */
        $point = $n >= 0
                ? $this->geometryN($n)
                : $this->geometryN(count($this->components) - (int) abs($n + 1));
        
        return $point;
    }

    /**
     * @return Point
     */
    public function getCentroid(): Point
    {
        if ($this->isEmpty()) {
            return new Point();
        }

        $geosObj = $this->getGeos();
        if (is_object($geosObj)) {
            // @codeCoverageIgnoreStart
            /** @noinspection PhpUndefinedMethodInspection */
            /** @var Point|null $geometry */
            $geometry = geoPHP::geosToGeometry($geosObj->centroid());
            return $geometry !== null ? $geometry : new Point();
            // @codeCoverageIgnoreEnd
        }

        $x = 0;
        $y = 0;
        $length = 0.0;
        $points = $this->getPoints();
        $numPoints = count($points)-1;
        for ($i=0; $i<$numPoints; ++$i) {
            $currX = $points[$i]->getX();
            $currY = $points[$i]->getY();
            $nextX = $points[$i+1]->getX();
            $nextY = $points[$i+1]->getY();
            
            $dx = $nextX - $currX;
            $dy = $nextY - $currY;
            $segmentLength = sqrt($dx*$dx + $dy*$dy);
            $length += $segmentLength;
            $x += ($currX + $nextX) / 2 * $segmentLength;
            $y += ($currY + $nextY) / 2 * $segmentLength;
        }

        return $length === 0.0 ? $this->startPoint() : new Point($x / $length, $y / $length);

    }

    /**
     * Returns the length of this Curve in its associated spatial reference.
     * E.g. if Geometry is in geographical coordinate system it returns the length in degrees.
     *
     * @return float
     */
    public function getLength(): float
    {
        $geosObj = $this->getGeos();
        if (is_object($geosObj)) {
            // @codeCoverageIgnoreStart
            /** @noinspection PhpUndefinedMethodInspection */
            return $geosObj->length();
            // @codeCoverageIgnoreEnd
        }
        
        $length = 0.0;
        $points = $this->getPoints();
        $numPoints = count($points)-1;
        for ($i=0; $i<$numPoints; ++$i) {
            $length += sqrt(
                pow(($points[$i]->getX() - $points[$i+1]->getX()), 2) +
                pow(($points[$i]->getY() - $points[$i+1]->getY()), 2)
            );
        }
        
        return $length;
    }

    /**
     * Returns the length of a 3-dimensional geometry.
     *
     * @return float
     */
    public function length3D(): float
    {
        $length = 0.0;
        
        $previousPoint = null;
        foreach ($this->getPoints() as $point) {
            if ($previousPoint) {
                $length += sqrt(
                    pow(($previousPoint->getX() - $point->getX()), 2) +
                    pow(($previousPoint->getY() - $point->getY()), 2) +
                    pow(($previousPoint->getZ() - $point->getZ()), 2)
                );
            }
            /** @var Point $previousPoint */
            $previousPoint = $point;
        }
        return $length;
    }

    /**
     * @param  float|int $radius Default is the semi-major axis of WGS84.
     * @return float length in meters
     */
    public function greatCircleLength($radius = geoPHP::EARTH_WGS84_SEMI_MAJOR_AXIS): float
    {
        $length = 0.0;
        $rad = M_PI / 180;
        $points = $this->getPoints();
        $numPoints = $this->numPoints() - 1;
        for ($i = 0; $i < $numPoints; ++$i) {
            // Simplified Vincenty formula with equal major and minor axes (a sphere)
            $lat1 = $points[$i]->getY() * $rad;
            $lat2 = $points[$i + 1]->getY() * $rad;
            $lon1 = $points[$i]->getX() * $rad;
            $lon2 = $points[$i + 1]->getX() * $rad;
            $deltaLon = $lon2 - $lon1;
            $d = $radius *
                atan2(
                    sqrt(
                        pow(cos($lat2) * sin($deltaLon), 2) +
                        pow(cos($lat1) * sin($lat2) - sin($lat1) * cos($lat2) * cos($deltaLon), 2)
                    ),
                    sin($lat1) * sin($lat2) +
                    cos($lat1) * cos($lat2) * cos($deltaLon)
                );
            if ($points[$i]->is3D()) {
                $d = sqrt(
                    pow($d, 2) +
                    pow($points[$i + 1]->getZ() - $points[$i]->getZ(), 2)
                );
            }

            $length += $d;
        }
        
        return $length;
    }

    /**
     * @return float Haversine length of geometry in degrees
     */
    public function haversineLength(): float
    {
        $distance = 0.0;
        $points = $this->getPoints();
        $numPoints = $this->numPoints() - 1;
        for ($i = 0; $i < $numPoints; ++$i) {
            $point = $points[$i];
            $nextPoint = $points[$i + 1];
            $degree = (geoPHP::EARTH_WGS84_SEMI_MAJOR_AXIS *
                acos(
                    sin(deg2rad($point->getY() ?? 0)) * sin(deg2rad($nextPoint->getY() ?? 0)) +
                    cos(deg2rad($point->getY() ?? 0)) * cos(deg2rad($nextPoint->getY() ?? 0)) *
                    cos(deg2rad(abs($point->getX() - $nextPoint->getX())))
                )
            );
            if (!is_nan($degree)) {
                $distance += $degree;
            }
        }
        return $distance;
    }

    /**
     * @source  https://github.com/mjaschen/phpgeo/blob/master/src/Location/Distance/Vincenty.php
     * @author  Marcus Jaschen <[email protected]>
     * @license https://opensource.org/licenses/GPL-3.0 GPL
     * (note: geoPHP uses "GPL version 2 (or later)" license which is compatible with GPLv3)
     *
     * @return float Length in meters
     */
    public function vincentyLength(): float
    {
        $length = 0.0;
        $rad = M_PI / 180;
        $points = $this->getPoints();
        $numPoints = count($points) - 1;
        for ($i = 0; $i < $numPoints; ++$i) {
            // Inverse Vincenty formula
            $lat1 = $points[$i]->getY() * $rad;
            $lat2 = $points[$i + 1]->getY() * $rad;
            $lng1 = $points[$i]->getX() * $rad;
            $lng2 = $points[$i + 1]->getX() * $rad;

            $a = geoPHP::EARTH_WGS84_SEMI_MAJOR_AXIS;
            $b = geoPHP::EARTH_WGS84_SEMI_MINOR_AXIS;
            $f = 1 / geoPHP::EARTH_WGS84_FLATTENING;
            $L = $lng2 - $lng1;
            $U1 = atan((1 - $f) * tan($lat1));
            $U2 = atan((1 - $f) * tan($lat2));
            $iterationLimit = 100;
            $lambda = $L;
            $sinU1 = sin($U1);
            $sinU2 = sin($U2);
            $cosU1 = cos($U1);
            $cosU2 = cos($U2);
            do {
                $sinLambda = sin($lambda);
                $cosLambda = cos($lambda);
                $sinSigma = sqrt(
                    ($cosU2 * $sinLambda) *
                    ($cosU2 * $sinLambda) +
                    ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) *
                    ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda)
                );
                if ($sinSigma == 0) {
                    return 0.0;
                }
                $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
                $sigma = atan2($sinSigma, $cosSigma);
                $sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
                $cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
                $cos2SigmaM = 0;
                if ($cosSqAlpha <> 0) {
                    $cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;
                }
                $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
                $lambdaP = $lambda;
                $lambda = $L + (1 - $C) * $f * $sinAlpha *
                    ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (- 1 + 2 * $cos2SigmaM * $cos2SigmaM)));
            } while (abs($lambda - $lambdaP) > 1e-12 && --$iterationLimit > 0);
            if ($iterationLimit == 0) {
                return 0.0; // not converging
            }
            $uSq = $cosSqAlpha * ($a * $a - $b * $b) / ($b * $b);
            $A = 1 + $uSq / 16384 * (4096 + $uSq * (- 768 + $uSq * (320 - 175 * $uSq)));
            $B = $uSq / 1024 * (256 + $uSq * (- 128 + $uSq * (74 - 47 * $uSq)));
            $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 *
                    ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6
                        * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma)
                        * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)));

            $length += $b * $A * ($sigma - $deltaSigma);
        }
        // Returns length in meters.
        return $length;
    }

    /**
     * @return int|float|null
     */
    public function minimumZ()
    {
        $min = PHP_INT_MAX;
        foreach ($this->getPoints() as $point) {
            if (null !== ($z = $point->getZ()) && $z < $min) {
                $min = $z;
            }
        }
        return $min !== PHP_INT_MAX ? $min : null;
    }

    /**
     * @return int|float|null
     */
    public function maximumZ()
    {
        $max = PHP_INT_MIN;
        foreach ($this->getPoints() as $point) {
            if (null !== ($z = $point->getZ()) && $z > $max) {
                $max = $z;
            }
        }

        return $max !== PHP_INT_MIN ? $max : null;
    }

    /**
     * @return int|float|null
     */
    public function zDifference()
    {
        $startPt = $this->startPoint();
        $endPt = $this->endPoint();
        
        if ($startPt->hasZ() && $endPt->hasZ()) {
            return abs($startPt->getZ() - $endPt->getZ());
        }
        
        return null;
    }

    /**
     * Returns the cumulative elevation gain of the LineString
     *
     * @param int|float|null $verticalTolerance Smoothing factor filtering noisy elevation data.
     *                                          Its unit equals to the z-coordinates unit (meters for geographical coordinates)
     *                                          If the elevation data comes from a DEM, a value around 3.5 can be acceptable.
     *
     * @return float
     */
    public function elevationGain($verticalTolerance = 0)
    {
        $gain = 0.0;
        $lastEle = $this->startPoint()->getZ();
        $numPoints = $this->numPoints();
        
        foreach ($this->getPoints() as $i => $point) {
            if (abs($point->getZ() - $lastEle) > $verticalTolerance || $i === $numPoints - 1) {
                if ($point->getZ() > $lastEle) {
                    $gain += $point->getZ() - $lastEle;
                }
                $lastEle = $point->getZ();
            }
        }
        
        return $gain;
    }

    /**
     * Returns the cumulative elevation loss of the LineString
     *
     * @param int|float|null $verticalTolerance Smoothing factor filtering noisy elevation data.
     *                                          Its unit equals to the z-coordinates unit (meters for geographical coordinates)
     *                                          If the elevation data comes from a DEM, a value around 3.5 can be acceptable.
     *
     * @return float
     */
    public function elevationLoss($verticalTolerance = 0)
    {
        $loss = 0.0;
        $lastEle = $this->startPoint()->getZ();
        $numPoints = $this->numPoints();
        
        foreach ($this->getPoints() as $i => $point) {
            if (abs($point->getZ() - $lastEle) > $verticalTolerance || $i === $numPoints - 1) {
                if ($point->getZ() < $lastEle) {
                    $loss += $lastEle - $point->getZ();
                }
                $lastEle = $point->getZ();
            }
        }
        
        return $loss;
    }

    public function minimumM()
    {
        $min = PHP_INT_MAX;
        foreach ($this->getPoints() as $point) {
            if ($point->isMeasured() && $point->m() < $min) {
                $min = $point->m();
            }
        }
        return $min !== PHP_INT_MAX ? $min : null;
    }

    public function maximumM()
    {
        $max = PHP_INT_MIN;
        foreach ($this->getPoints() as $point) {
            if ($point->isMeasured() && $point->m() > $max) {
                $max = $point->m();
            }
        }

        return $max !== PHP_INT_MIN ? $max : null;
    }

    /**
     * Get all line segments
     *
     * @param  bool $toArray return segments as LineString or array of start and end points.
     * @return LineString[]|Point[][]
     */
    public function explode(bool $toArray = false): array
    {
        $points = $this->getPoints();
        $numPoints = count($points);
        if ($numPoints < 2) {
            return [];
        }
        $parts = [];
        for ($i = 1; $i < $numPoints; ++$i) {
            $segment = [$points[$i - 1], $points[$i]];
            $parts[] = $toArray ? $segment : new LineString($segment);
        }
        return $parts;
    }

    /**
     * Checks that LineString is a Simple Geometry
     *
     * @return bool
     */
    public function isSimple(): bool
    {
        $geosObj = $this->getGeos();
        if (is_object($geosObj)) {
            // @codeCoverageIgnoreStart
            /** @noinspection PhpUndefinedMethodInspection */
            return $geosObj->isSimple();
            // @codeCoverageIgnoreEnd
        }
        
        $segments = $this->explode(true);
        /** @var Point[][] $segments */
        foreach ($segments as $i => $segment) {
            foreach ($segments as $j => $checkSegment) {
                if ($i != $j) {
                    if (Geometry::segmentIntersects($segment[0], $segment[1], $checkSegment[0], $checkSegment[1])) {
                        return false;
                    }
                }
            }
        }
        
        return true;
    }
    
    /**
     * @return bool
     */
    public function isValid(): bool
    {
        $geosObj = $this->getGeos();
        if (is_object($geosObj)) {
            /** @noinspection PhpUndefinedMethodInspection */
            return $geosObj->checkValidity()['valid'];
        }
        
        // there should only be unique points but there have to be at least 2 of them
        $pts = [];
        foreach ($this->components as $pt) {
            if ($pt->isEmpty()) {
                continue;
            }
            $pts[] = $pt->asText();
        }
        if (count(array_unique($pts)) < 2) {
            return false;
        }
        
        return $this->isSimple();
    }

    /**
     * @param  LineString $segment
     * @return bool
     */
    public function lineSegmentIntersect($segment): bool
    {
        return Geometry::segmentIntersects(
            $this->startPoint() ?? new Point(),
            $this->endPoint() ?? new Point(),
            $segment->startPoint() ?? new Point(),
            $segment->endPoint() ?? new Point()
        );
    }

    /**
     * @param  Geometry|Collection $geometry
     * @return float|null
     */
    public function distance(Geometry $geometry)
    {
        $geosObj = $this->getGeos();
        if (is_object($geosObj)) {
            // @codeCoverageIgnoreStart
            /** @noinspection PhpUndefinedMethodInspection */
            $geosObj2 = $geometry->getGeos();
            return $geosObj2 !== false ? $geosObj->distance($geosObj2) : null;
            // @codeCoverageIgnoreEnd
        }

        if ($geometry->geometryType() === Geometry::POINT) {
            // This is defined in the Point class nicely
            return $geometry->distance($this);
        }
        
        if ($geometry->geometryType() === Geometry::LINESTRING) {
            /** @var LineString $geometry */
            return $this->distanceToLinestring($geometry);
        }
        
        // It can be treated as a collection
        return parent::distance($geometry);
    }
    
    /**
     * @param  LineString $geometry
     * @return float
     */
    private function distanceToLinestring(LineString $geometry): float
    {
        $distance = PHP_INT_MAX;
        $geometrySegments = $geometry->explode();
        
        /** @var LineString $seg1 */
        foreach ($this->explode() as $seg1) {
            /** @var LineString $seg2 */
            foreach ($geometrySegments as $seg2) {
                if ($seg1->lineSegmentIntersect($seg2)) {
                    return 0.0;
                }
                // Because line-segments are straight, the shortest distance will occur at an endpoint.
                // If they are parallel, an endpoint calculation is still accurate.
                $distance = min(
                    $distance,
                    $seg1->startPoint()->distance($seg2),
                    $seg1->endPoint()->distance($seg2),
                    $seg2->startPoint()->distance($seg1),
                    $seg2->endPoint()->distance($seg1)
                );

                if ($distance === 0.0) {
                    return 0.0;
                }
            }
        }
        return (float) $distance;
    }
}


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 2021-11-10 08:46 UTC by Report Bug Copy Issue Report The condition `$length === 0.0` is always `true`. 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->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			$lastEle = $this->startPoint()->getZ();
356			$numPoints = $this->numPoints();
357
358			foreach ($this->getPoints() as $i => $point) {
359			if (abs($point->getZ() - $lastEle) > $verticalTolerance \|\| $i === $numPoints - 1) {
360			if ($point->getZ() > $lastEle) {
361			$gain += $point->getZ() - $lastEle;
362			}
363			$lastEle = $point->getZ();
364			}
365			}
366
367			return $gain;
368			}
369
370			/**
371			* Returns the cumulative elevation loss of the LineString
372			*
373			* @param int\|float\|null $verticalTolerance Smoothing factor filtering noisy elevation data.
374			* Its unit equals to the z-coordinates unit (meters for geographical coordinates)
375			* If the elevation data comes from a DEM, a value around 3.5 can be acceptable.
376			*
377			* @return float
378			*/
379			public function elevationLoss($verticalTolerance = 0)
380			{
381			$loss = 0.0;
382			$lastEle = $this->startPoint()->getZ();
383			$numPoints = $this->numPoints();
384
385			foreach ($this->getPoints() as $i => $point) {
386			if (abs($point->getZ() - $lastEle) > $verticalTolerance \|\| $i === $numPoints - 1) {
387			if ($point->getZ() < $lastEle) {
388			$loss += $lastEle - $point->getZ();
389			}
390			$lastEle = $point->getZ();
391			}
392			}
393
394			return $loss;
395			}
396
397			public function minimumM()
398			{
399			$min = PHP_INT_MAX;
400			foreach ($this->getPoints() as $point) {
401			if ($point->isMeasured() && $point->m() < $min) {
402			$min = $point->m();
403			}
404			}
405			return $min !== PHP_INT_MAX ? $min : null;
406			}
407
408			public function maximumM()
409			{
410			$max = PHP_INT_MIN;
411			foreach ($this->getPoints() as $point) {
412			if ($point->isMeasured() && $point->m() > $max) {
413			$max = $point->m();
414			}
415			}
416
417			return $max !== PHP_INT_MIN ? $max : null;
418			}
419
420			/**
421			* Get all line segments
422			*
423			* @param bool $toArray return segments as LineString or array of start and end points.
424			* @return LineString[]\|Point[][]
425			*/
426			public function explode(bool $toArray = false): array
427			{
428			$points = $this->getPoints();
429			$numPoints = count($points);
430			if ($numPoints < 2) {
431			return [];
432			}
433			$parts = [];
434			for ($i = 1; $i < $numPoints; ++$i) {
435			$segment = [$points[$i - 1], $points[$i]];
436			$parts[] = $toArray ? $segment : new LineString($segment);
437			}
438			return $parts;
439			}
440
441			/**
442			* Checks that LineString is a Simple Geometry
443			*
444			* @return bool
445			*/
446			public function isSimple(): bool
447			{
448			$geosObj = $this->getGeos();
449			if (is_object($geosObj)) {
450			// @codeCoverageIgnoreStart
451			/** @noinspection PhpUndefinedMethodInspection */
452			return $geosObj->isSimple();
453			// @codeCoverageIgnoreEnd
454			}
455
456			$segments = $this->explode(true);
457			/** @var Point[][] $segments */
458			foreach ($segments as $i => $segment) {
459			foreach ($segments as $j => $checkSegment) {
460			if ($i != $j) {
461			if (Geometry::segmentIntersects($segment[0], $segment[1], $checkSegment[0], $checkSegment[1])) {
462			return false;
463			}
464			}
465			}
466			}
467
468			return true;
469			}
470
471			/**
472			* @return bool
473			*/
474			public function isValid(): bool
475			{
476			$geosObj = $this->getGeos();
477			if (is_object($geosObj)) {
478			/** @noinspection PhpUndefinedMethodInspection */
479			return $geosObj->checkValidity()['valid'];
480			}
481
482			// there should only be unique points but there have to be at least 2 of them
483			$pts = [];
484			foreach ($this->components as $pt) {
485			if ($pt->isEmpty()) {
486			continue;
487			}
488			$pts[] = $pt->asText();
489			}
490			if (count(array_unique($pts)) < 2) {
491			return false;
492			}
493
494			return $this->isSimple();
495			}
496
497			/**
498			* @param LineString $segment
499			* @return bool
500			*/
501			public function lineSegmentIntersect($segment): bool
502			{
503			return Geometry::segmentIntersects(
504			$this->startPoint() ?? new Point(),
505			$this->endPoint() ?? new Point(),
506			$segment->startPoint() ?? new Point(),
507			$segment->endPoint() ?? new Point()
508			);
509			}
510
511			/**
512			* @param Geometry\|Collection $geometry
513			* @return float\|null
514			*/
515			public function distance(Geometry $geometry)
516			{
517			$geosObj = $this->getGeos();
518			if (is_object($geosObj)) {
519			// @codeCoverageIgnoreStart
520			/** @noinspection PhpUndefinedMethodInspection */
521			$geosObj2 = $geometry->getGeos();
522			return $geosObj2 !== false ? $geosObj->distance($geosObj2) : null;
523			// @codeCoverageIgnoreEnd
524			}
525
526			if ($geometry->geometryType() === Geometry::POINT) {
527			// This is defined in the Point class nicely
528			return $geometry->distance($this);
529			}
530
531			if ($geometry->geometryType() === Geometry::LINESTRING) {
532			/** @var LineString $geometry */
533			return $this->distanceToLinestring($geometry);
534			}
535
536			// It can be treated as a collection
537			return parent::distance($geometry);
538			}
539
540			/**
541			* @param LineString $geometry
542			* @return float
543			*/
544			private function distanceToLinestring(LineString $geometry): float
545			{
546			$distance = PHP_INT_MAX;
547			$geometrySegments = $geometry->explode();
548
549			/** @var LineString $seg1 */
550			foreach ($this->explode() as $seg1) {
551			/** @var LineString $seg2 */
552			foreach ($geometrySegments as $seg2) {
553			if ($seg1->lineSegmentIntersect($seg2)) {
554			return 0.0;
555			}
556			// Because line-segments are straight, the shortest distance will occur at an endpoint.
557			// If they are parallel, an endpoint calculation is still accurate.
558			$distance = min(
559			$distance,
560			$seg1->startPoint()->distance($seg2),
561			$seg1->endPoint()->distance($seg2),
562			$seg2->startPoint()->distance($seg1),
563			$seg2->endPoint()->distance($seg1)
564			);
565
566			if ($distance === 0.0) {
567			return 0.0;
568			}
569			}
570			}
571			return (float) $distance;
572			}
573			}
574

swen100 / geoPHP

Push — master ( 0fd618...5e6b25 )

LineString::haversineLength() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like