Point::vincenty() - Code Metrics - Inspection of "Update EPSG data" - dvdoug/PHPCoord - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( bbeaaa...9ed6c0 )

by Doug

created 2022-09-25 19:55 UTC

Point::vincenty() B

↳ Parent: Point

Complexity

Conditions	9
Paths	20

Size

Total Lines	62
Code Lines	46

Duplication

Lines	0
Ratio	0 %

Code Coverage

Tests	43
CRAP Score	9.0009

Importance

Changes

Metric	Value
cc	9
eloc	46
nc	20
nop	3
dl	0
loc	62
ccs	43
cts	44
cp	0.9773
crap	9.0009
rs	7.6226
c	0
b	0
f	0

How to fix Long Method

<?php
/**
 * PHPCoord.
 *
 * @author Doug Wright
 */
declare(strict_types=1);

namespace PHPCoord;

use DateTimeImmutable;
use PHPCoord\CoordinateOperation\CoordinateOperationMethods;
use PHPCoord\CoordinateOperation\CoordinateOperations;
use PHPCoord\CoordinateOperation\GeographicValue;
use PHPCoord\CoordinateReferenceSystem\Compound;
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
use PHPCoord\CoordinateReferenceSystem\Geocentric;
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
use PHPCoord\CoordinateReferenceSystem\Geographic3D;
use PHPCoord\CoordinateReferenceSystem\Projected;
filter:
    dependency_paths: ["lib/*"]
use PHPCoord\CoordinateReferenceSystem\Vertical;
use PHPCoord\Datum\Ellipsoid;
use PHPCoord\UnitOfMeasure\Angle\Angle;
use PHPCoord\UnitOfMeasure\Length\Length;
use PHPCoord\UnitOfMeasure\Length\Metre;
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
use PHPCoord\UnitOfMeasure\Scale\Scale;
use PHPCoord\UnitOfMeasure\UnitOfMeasure;
use PHPCoord\UnitOfMeasure\UnitOfMeasureFactory;
use Stringable;

use function abs;
use function acos;
use function asin;
use function atan;
use function atan2;
use function cos;
use function sin;
use function sqrt;
use function sscanf;
use function str_starts_with;
use function tan;

use const M_PI;

abstract class Point implements Stringable
{
    protected const ITERATION_CONVERGENCE_FORMULA = 1e-10;
    protected const ITERATION_CONVERGENCE_GRID = 0.0001;
    protected const METHODS_REQUIRING_HORIZONTAL_POINT = [
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE,
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN,
    ];
    protected const METHODS_THAT_REQUIRE_DIRECTION = [
        CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION => CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
        CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION => CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
        CoordinateOperationMethods::EPSG_NADCON5_2D => CoordinateOperationMethods::EPSG_NADCON5_2D,
        CoordinateOperationMethods::EPSG_NADCON5_3D => CoordinateOperationMethods::EPSG_NADCON5_3D,
        CoordinateOperationMethods::EPSG_NTV2 => CoordinateOperationMethods::EPSG_NTV2,
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
        CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN => CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN,
    ];

    private static array $gridCache = [];

    /**
     * @internal
     */
    public function performOperation(string $srid, Compound|Geocentric|Geographic2D|Geographic3D|Projected|Vertical $to, bool $inReverse, array $additionalParams = []): self
    {
        $operation = CoordinateOperations::getOperationData($srid);

        if ($operation['method'] === CoordinateOperationMethods::EPSG_ALIAS) {
            $point = clone $this;
            $point->crs = $to;


            return $point;
        } else {
            $method = CoordinateOperationMethods::getFunctionName($operation['method']);
            $params = self::resolveParamsByOperation($srid, $operation['method'], $inReverse);

            if (isset(self::METHODS_REQUIRING_HORIZONTAL_POINT[$operation['method']])) {
                $params['horizontalPoint'] = $additionalParams['horizontalPoint'];
            }

            return $this->$method($to, ...$params);
        }
    }

    protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
    {
        $params = [];
        $powerCoefficients = [];
        foreach (CoordinateOperations::getParamData($operationSrid) as $paramName => $paramData) {
            if (isset($paramData['fileProvider'])) {
                $params[$paramName] = static::$gridCache[$paramData['fileProvider']] ??= (new $paramData['fileProvider']())->provideGrid();

            } else {
                if ($inReverse && $paramData['reverses']) {
                    $paramData['value'] *= -1;
                }
                if ($paramData['uom']) {
                    $param = UnitOfMeasureFactory::makeUnit($paramData['value'], $paramData['uom']);
                } else {
                    $param = $paramData['value'];
                }
                if (str_starts_with($paramName, 'Au') || str_starts_with($paramName, 'Bu')) {
                    $powerCoefficients[$paramName] = $param;
                } else {
                    $params[$paramName] = $param;
                }
            }
        }
        if ($powerCoefficients) {
            $params['powerCoefficients'] = $powerCoefficients;
        }
        if (isset(self::METHODS_THAT_REQUIRE_DIRECTION[$methodSrid])) {
            $params['inReverse'] = $inReverse;
        }

        return $params;
    }

    protected static function sign(float $number): int
    {
        if ($number < 0) {
            return -1;
        }

        return 1;
    }

    /**
     * Calculate surface distance between two points.
     */
    protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length
    {
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
        $b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
        $f = $ellipsoid->getFlattening();
        $U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue()));
        $U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue()));
        $L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue();

        $lambda = $L;
        do {
            $lambdaN = $lambda;

            $sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2);
            $cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda);
            $sigma = atan2($sinSigma, $cosSigma);

            $sinAlpha = $sinSigma ? (cos($U1) * cos($U2) * sin($lambda) / $sinSigma) : 0;
            $cosSqAlpha = (1 - $sinAlpha ** 2);
            $cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0;
            $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
            $lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
        } while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE_FORMULA && abs($lambda) < M_PI);

        // Antipodal case
        if (abs($lambda) >= M_PI) {
            if ($L >= 0) {
                $LPrime = M_PI - $L;
            } else {
                $LPrime = -M_PI - $L;
            }

            $lambdaPrime = 0;
            $sigma = M_PI - abs($U1 + $U2);
            $sinSigma = sin($sigma);
            $cosSqAlpha = 0.5;
            $sinAlpha = 0;

            do {
                $sinAlphaN = $sinAlpha;

                $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
                $cos2SigmaM = $cosSqAlpha ? cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha : 0;
                $D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2)));
                $sinAlpha = ($LPrime - $lambdaPrime) / $D;
                $cosSqAlpha = (1 - $sinAlpha ** 2);
                $sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2));
                $lambdaPrime = self::asin($sinLambdaPrime);
                $sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2;
                $sinSigma = sqrt($sinSqSigma);
            } while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE_FORMULA);
        }

        $E = sqrt(1 + (($a ** 2 - $b ** 2) / $b ** 2) * $cosSqAlpha);
        $F = ($E - 1) / ($E + 1);

        $A = (1 + $F ** 2 / 4) / (1 - $F);
        $B = $F * (1 - 3 / 8 * $F ** 2);

        $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2)));

        return new Metre($b * $A * ($sigma - $deltaSigma));
    }

    /**
     * General polynomial.
     * @param Coefficient[] $powerCoefficients
     */
    protected function generalPolynomialUnitless(
        float $xs,
        float $ys,
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
        Scale $scalingFactorForSourceCRSCoordDifferences,
        Scale $scalingFactorForTargetCRSCoordDifferences,
        Scale $A0,
        Scale $B0,
        array $powerCoefficients
    ): array {
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();

        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);

        $mTdX = $A0->getValue();
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
            if ($coefficientName[0] === 'A') {
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
            }
        }

        $mTdY = $B0->getValue();
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
            if ($coefficientName[0] === 'B') {
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
            }
        }

        $xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
        $yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();

        return ['xt' => $xt, 'yt' => $yt];
    }

    /**
     * Reversible polynomial.
     */
    protected function reversiblePolynomialUnitless(
        float $xs,
        float $ys,
        Angle $ordinate1OfEvaluationPoint,
        Angle $ordinate2OfEvaluationPoint,
        Scale $scalingFactorForCoordDifferences,
        Scale $A0,
        Scale $B0,
        array $powerCoefficients
    ): array {
        $xo = $ordinate1OfEvaluationPoint->getValue();
        $yo = $ordinate2OfEvaluationPoint->getValue();

        $U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
        $V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);

        $mTdX = $A0->getValue();
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
            if ($coefficientName[0] === 'A') {
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
            }
        }

        $mTdY = $B0->getValue();
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
            if ($coefficientName[0] === 'B') {
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
            }
        }

        $xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
        $yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();

        return ['xt' => $xt, 'yt' => $yt];
    }

    /**
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
     */
    protected static function acos(float $num): float
    {
        if ($num > 1.0) {
            $num = 1.0;
        } elseif ($num < -1) {
            $num = -1.0;
        }

        return acos($num);
    }

    /**
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
     */
    protected static function asin(float $num): float
    {
        if ($num > 1.0) {
            $num = 1.0;
        } elseif ($num < -1.0) {
            $num = -1.0;
        }

        return asin($num);
    }

    abstract public function getCRS(): CoordinateReferenceSystem;

    abstract public function getCoordinateEpoch(): ?DateTimeImmutable;

    abstract public function calculateDistance(self $to): Length;
}


1		<?php
2		/**
3		* PHPCoord.
4		*
5		* @author Doug Wright
6		*/
7		declare(strict_types=1);
8
9		namespace PHPCoord;
10
11		use DateTimeImmutable;
12		use PHPCoord\CoordinateOperation\CoordinateOperationMethods;
13		use PHPCoord\CoordinateOperation\CoordinateOperations;
14		use PHPCoord\CoordinateOperation\GeographicValue;
15		use PHPCoord\CoordinateReferenceSystem\Compound;
16		use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
17		use PHPCoord\CoordinateReferenceSystem\Geocentric;
18		use PHPCoord\CoordinateReferenceSystem\Geographic2D;
19		use PHPCoord\CoordinateReferenceSystem\Geographic3D;
20		use PHPCoord\CoordinateReferenceSystem\Projected;
		0 ignored issues – show Bug introduced 2021-10-10 00:52 UTC by Report Bug Copy Issue Report The type `PHPCoord\CoordinateReferenceSystem\Projected` was not found. Maybe you did not declare it correctly or list all dependencies? The issue could also be caused by a filter entry in the build configuration. If the path has been excluded in your configuration, e.g. `excluded_paths: ["lib/"]`, you can move it to the dependency path list as follows: filter: dependency_paths: ["lib/"] For further information see https://scrutinizer-ci.com/docs/tools/php/php-scrutinizer/#list-dependency-paths Loading history...
21		use PHPCoord\CoordinateReferenceSystem\Vertical;
22		use PHPCoord\Datum\Ellipsoid;
23		use PHPCoord\UnitOfMeasure\Angle\Angle;
24		use PHPCoord\UnitOfMeasure\Length\Length;
25		use PHPCoord\UnitOfMeasure\Length\Metre;
26		use PHPCoord\UnitOfMeasure\Scale\Coefficient;
27		use PHPCoord\UnitOfMeasure\Scale\Scale;
28		use PHPCoord\UnitOfMeasure\UnitOfMeasure;
29		use PHPCoord\UnitOfMeasure\UnitOfMeasureFactory;
30		use Stringable;
31
32		use function abs;
33		use function acos;
34		use function asin;
35		use function atan;
36		use function atan2;
37		use function cos;
38		use function sin;
39		use function sqrt;
40		use function sscanf;
41		use function str_starts_with;
42		use function tan;
43
44		use const M_PI;
45
46		abstract class Point implements Stringable
47		{
48		protected const ITERATION_CONVERGENCE_FORMULA = 1e-10;
49		protected const ITERATION_CONVERGENCE_GRID = 0.0001;
50		protected const METHODS_REQUIRING_HORIZONTAL_POINT = [
51		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE,
52		CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
53		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
54		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
55		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
56		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN,
57		];
58		protected const METHODS_THAT_REQUIRE_DIRECTION = [
59		CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION => CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
60		CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION => CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
61		CoordinateOperationMethods::EPSG_NADCON5_2D => CoordinateOperationMethods::EPSG_NADCON5_2D,
62		CoordinateOperationMethods::EPSG_NADCON5_3D => CoordinateOperationMethods::EPSG_NADCON5_3D,
63		CoordinateOperationMethods::EPSG_NTV2 => CoordinateOperationMethods::EPSG_NTV2,
64		CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
65		CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN => CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN,
66		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
67		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
68		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
69		CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN,
70		];
71
72		private static array $gridCache = [];
73
74		/**
75		* @internal
76		*/
77	5224	public function performOperation(string $srid, Compound\|Geocentric\|Geographic2D\|Geographic3D\|Projected\|Vertical $to, bool $inReverse, array $additionalParams = []): self
78		{
79	5224	$operation = CoordinateOperations::getOperationData($srid);
80
81	5224	if ($operation['method'] === CoordinateOperationMethods::EPSG_ALIAS) {
82		$point = clone $this;
83		$point->crs = $to;
		0 ignored issues – show Bug Best Practice introduced 2021-11-11 02:03 UTC by Report Bug Copy Issue Report The property `crs` does not exist. Although not strictly required by PHP, it is generally a best practice to declare properties explicitly. Loading history...
84
85		return $point;
86		} else {
87	5224	$method = CoordinateOperationMethods::getFunctionName($operation['method']);
88	5224	$params = self::resolveParamsByOperation($srid, $operation['method'], $inReverse);
89
90	5224	if (isset(self::METHODS_REQUIRING_HORIZONTAL_POINT[$operation['method']])) {
91	9	$params['horizontalPoint'] = $additionalParams['horizontalPoint'];
92		}
93
94	5224	return $this->$method($to, ...$params);
95		}
96		}
97
98	5224	protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
99		{
100	5224	$params = [];
101	5224	$powerCoefficients = [];
102	5224	foreach (CoordinateOperations::getParamData($operationSrid) as $paramName => $paramData) {
103	5224	if (isset($paramData['fileProvider'])) {
104	12	$params[$paramName] = static::$gridCache[$paramData['fileProvider']] ??= (new $paramData['fileProvider']())->provideGrid();
		0 ignored issues – show Bug introduced 2021-11-11 02:03 UTC by Report Bug Copy Issue Report Since `$gridCache` is declared private, accessing it with `static` will lead to errors in possible sub-classes; you can either use `self`, or increase the visibility of `$gridCache` to at least protected. Loading history...
105		} else {
106	5224	if ($inReverse && $paramData['reverses']) {
107	528	$paramData['value'] *= -1;
108		}
109	5224	if ($paramData['uom']) {
110	5224	$param = UnitOfMeasureFactory::makeUnit($paramData['value'], $paramData['uom']);
111		} else {
112	11	$param = $paramData['value'];
113		}
114	5224	if (str_starts_with($paramName, 'Au') \|\| str_starts_with($paramName, 'Bu')) {
115	27	$powerCoefficients[$paramName] = $param;
116		} else {
117	5224	$params[$paramName] = $param;
118		}
119		}
120		}
121	5224	if ($powerCoefficients) {
122	27	$params['powerCoefficients'] = $powerCoefficients;
123		}
124	5224	if (isset(self::METHODS_THAT_REQUIRE_DIRECTION[$methodSrid])) {
125	6	$params['inReverse'] = $inReverse;
126		}
127
128	5224	return $params;
129		}
130
131	802	protected static function sign(float $number): int
132		{
133	802	if ($number < 0) {
134		return -1;
135		}
136
137	802	return 1;
138		}
139
140		/**
141		* Calculate surface distance between two points.
142		*/
143	162	protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length
144		{
145	162	$a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
146	162	$b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
147	162	$f = $ellipsoid->getFlattening();
148	162	$U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue()));
149	162	$U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue()));
150	162	$L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue();
151
152	162	$lambda = $L;
153		do {
154	162	$lambdaN = $lambda;
155
156	162	$sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2);
157	162	$cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda);
158	162	$sigma = atan2($sinSigma, $cosSigma);
159
160	162	$sinAlpha = $sinSigma ? (cos($U1) * cos($U2) * sin($lambda) / $sinSigma) : 0;
161	162	$cosSqAlpha = (1 - $sinAlpha ** 2);
162	162	$cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0;
163	162	$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
164	162	$lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
165	162	} while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE_FORMULA && abs($lambda) < M_PI);
166
167		// Antipodal case
168	162	if (abs($lambda) >= M_PI) {
169	27	if ($L >= 0) {
170	27	$LPrime = M_PI - $L;
171		} else {
172		$LPrime = -M_PI - $L;
173		}
174
175	27	$lambdaPrime = 0;
176	27	$sigma = M_PI - abs($U1 + $U2);
177	27	$sinSigma = sin($sigma);
178	27	$cosSqAlpha = 0.5;
179	27	$sinAlpha = 0;
180
181		do {
182	27	$sinAlphaN = $sinAlpha;
183
184	27	$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
185	27	$cos2SigmaM = $cosSqAlpha ? cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha : 0;
186	27	$D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2)));
187	27	$sinAlpha = ($LPrime - $lambdaPrime) / $D;
188	27	$cosSqAlpha = (1 - $sinAlpha ** 2);
189	27	$sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2));
190	27	$lambdaPrime = self::asin($sinLambdaPrime);
191	27	$sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2;
192	27	$sinSigma = sqrt($sinSqSigma);
193	27	} while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE_FORMULA);
194		}
195
196	162	$E = sqrt(1 + (($a 2 - $b 2) / $b ** 2) * $cosSqAlpha);
197	162	$F = ($E - 1) / ($E + 1);
198
199	162	$A = (1 + $F ** 2 / 4) / (1 - $F);
200	162	$B = $F * (1 - 3 / 8 * $F ** 2);
201
202	162	$deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2)));
203
204	162	return new Metre($b * $A * ($sigma - $deltaSigma));
205		}
206
207		/**
208		* General polynomial.
209		* @param Coefficient[] $powerCoefficients
210		*/
211	18	protected function generalPolynomialUnitless(
212		float $xs,
213		float $ys,
214		UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
215		UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
216		UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
217		UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
218		Scale $scalingFactorForSourceCRSCoordDifferences,
219		Scale $scalingFactorForTargetCRSCoordDifferences,
220		Scale $A0,
221		Scale $B0,
222		array $powerCoefficients
223		): array {
224	18	$xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
225	18	$yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
226	18	$xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
227	18	$yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
228
229	18	$U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
230	18	$V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
231
232	18	$mTdX = $A0->getValue();
233	18	foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
234	18	if ($coefficientName[0] === 'A') {
235	18	sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
236	18	$mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
237		}
238		}
239
240	18	$mTdY = $B0->getValue();
241	18	foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
242	18	if ($coefficientName[0] === 'B') {
243	18	sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
244	18	$mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
245		}
246		}
247
248	18	$xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
249	18	$yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
250
251	18	return ['xt' => $xt, 'yt' => $yt];
252		}
253
254		/**
255		* Reversible polynomial.
256		*/
257	36	protected function reversiblePolynomialUnitless(
258		float $xs,
259		float $ys,
260		Angle $ordinate1OfEvaluationPoint,
261		Angle $ordinate2OfEvaluationPoint,
262		Scale $scalingFactorForCoordDifferences,
263		Scale $A0,
264		Scale $B0,
265		array $powerCoefficients
266		): array {
267	36	$xo = $ordinate1OfEvaluationPoint->getValue();
268	36	$yo = $ordinate2OfEvaluationPoint->getValue();
269
270	36	$U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
271	36	$V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);
272
273	36	$mTdX = $A0->getValue();
274	36	foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
275	36	if ($coefficientName[0] === 'A') {
276	36	sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
277	36	$mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
278		}
279		}
280
281	36	$mTdY = $B0->getValue();
282	36	foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
283	36	if ($coefficientName[0] === 'B') {
284	36	sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
285	36	$mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
286		}
287		}
288
289	36	$xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
290	36	$yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();
291
292	36	return ['xt' => $xt, 'yt' => $yt];
293		}
294
295		/**
296		* Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
297		* silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
298		*/
299		protected static function acos(float $num): float
300		{
301		if ($num > 1.0) {
302		$num = 1.0;
303		} elseif ($num < -1) {
304		$num = -1.0;
305		}
306
307		return acos($num);
308		}
309
310		/**
311		* Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
312		* silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
313		*/
314	2652	protected static function asin(float $num): float
315		{
316	2652	if ($num > 1.0) {
317		$num = 1.0;
318	2652	} elseif ($num < -1.0) {
319		$num = -1.0;
320		}
321
322	2652	return asin($num);
323		}
324
325		abstract public function getCRS(): CoordinateReferenceSystem;
326
327		abstract public function getCoordinateEpoch(): ?DateTimeImmutable;
328
329		abstract public function calculateDistance(self $to): Length;
330		}
331

dvdoug / PHPCoord

Push — master ( bbeaaa...9ed6c0 )

Point::vincenty() B

Complexity

Size

Duplication

Code Coverage

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like