#66 Add "reverse" for general polynomial method

CHANGELOG.md

@@ -2,7 +2,10 @@
 
 ## [Unreleased]
 
+## [5.9.1] - 2024-08-15
+### Changed
 - Updates to data for USA and WGS84
+- Support for Irish polynomial transformation in the ETRS89 to TM75 direction (TM75 to ETRS89 was already supported)
 
 ## [5.9.0] - 2024-08-04
 ### Changed

docs/reflection/coordinateoperation/geographic2d.txt

@@ -18975,6 +18975,53 @@ to TM75 (Geographic2D)
  :header: "EPSG", "PHPCoord"
  :widths: 40, 60
 
+ "| Name: TM75 to ETRS89 (1)
+ | Code: ``urn:ogc:def:coordinateOperation:EPSG::1041``
+ | Extent: Europe - Ireland (Republic and Ulster) - onshore",".. code-block:: php
+
+ $point->generalPolynomial(
+ to: Geographic2D::fromSRID(Geographic2D::EPSG_TM75),
+ ordinate1OfEvaluationPointInSourceCRS: new Degree(53.5),
+ ordinate2OfEvaluationPointInSourceCRS: new Degree(-7.7),
+ ordinate1OfEvaluationPointInTargetCRS: new Degree(53.5),
+ ordinate2OfEvaluationPointInTargetCRS: new Degree(-7.7),
+ scalingFactorForSourceCRSCoordDifferences: new Unity(0.1),
+ scalingFactorForTargetCRSCoordDifferences: new Unity(3600),
+ A0: new Coefficient(0.763),
+ Au1v0: new Coefficient(-4.487),
+ Au0v1: new Coefficient(0.123),
+ Au2v0: new Coefficient(0.215),
+ Au1v1: new Coefficient(-0.515),
+ Au0v2: new Coefficient(0.183),
+ Au3v0: new Coefficient(-0.265),
+ Au2v1: new Coefficient(-0.57),
+ Au1v2: new Coefficient(0.414),
+ Au0v3: new Coefficient(-0.374),
+ Au3v1: new Coefficient(2.852),
+ Au2v2: new Coefficient(5.703),
+ Au1v3: new Coefficient(13.11),
+ Au3v2: new Coefficient(-61.678),
+ Au2v3: new Coefficient(113.743),
+ Au3v3: new Coefficient(-265.898),
+ B0: new Coefficient(-2.81),
+ Bu1v0: new Coefficient(-0.341),
+ Bu0v1: new Coefficient(-4.68),
+ Bu2v0: new Coefficient(1.196),
+ Bu1v1: new Coefficient(-0.119),
+ Bu0v2: new Coefficient(0.17),
+ Bu3v0: new Coefficient(-0.887),
+ Bu2v1: new Coefficient(4.877),
+ Bu1v2: new Coefficient(3.913),
+ Bu0v3: new Coefficient(2.163),
+ Bu3v1: new Coefficient(-46.666),
+ Bu2v2: new Coefficient(-27.795),
+ Bu1v3: new Coefficient(18.867),
+ Bu3v2: new Coefficient(-95.377),
+ Bu2v3: new Coefficient(-284.294),
+ Bu3v3: new Coefficient(-853.95)
+ )
+
+ "
  "| Name: TM75 to ETRS89 (2)
  | Code: ``urn:ogc:def:coordinateOperation:EPSG::1953``
  | Extent: Europe - Ireland (Republic and Ulster) - onshore",".. code-block:: php

src/CoordinateOperation/CoordinateOperationMethods.php

@@ -3018,7 +3018,7 @@ class CoordinateOperationMethods
  ],
  'urn:ogc:def:method:EPSG::9648' => [
  'name' => 'General polynomial of degree 6',
- 'reversible' => false,
+ 'reversible' => true,
  'paramData' => [
  'ordinate1OfEvaluationPointInSourceCRS' => [
  'reverses' => false,

src/EPSG/Import/EPSGImporter.php

@@ -148,6 +148,11 @@ public function dataFromSQLFiles(): void
  */
  $sqlite->exec('UPDATE epsg_usage SET extent_code = 1262 WHERE extent_code IN (1263, 2346, 2830, 4393, 4520, 4523)');
 
+ /*
+ * Ireland TM75 Polynomial is not "reversible" but can be reversed via iteration
+ */
+ $sqlite->exec('UPDATE epsg_coordoperationmethod SET reverse_op = 1 WHERE coord_op_method_code = 9648');
+
  $sqlite->close();
  }
 

src/Point/GeographicPoint.php

@@ -2041,7 +2041,8 @@ public function generalPolynomial(
  Scale $scalingFactorForTargetCRSCoordDifferences,
  Scale $A0,
  Scale $B0,
- array $powerCoefficients
+ array $powerCoefficients,
+ bool $inReverse
  ): self {
  $xs = $this->latitude->getValue();
  $ys = $this->longitude->getValue();
@@ -2057,7 +2058,8 @@ public function generalPolynomial(
  $scalingFactorForTargetCRSCoordDifferences,
  $A0,
  $B0,
- $powerCoefficients
+ $powerCoefficients,
+ $inReverse
  );
 
  $xtUnit = $to->getCoordinateSystem()->getAxes()[0]->getUnitOfMeasureId();

src/Point/Point.php

@@ -63,6 +63,8 @@ abstract class Point implements Stringable
  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_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN => CoordinateOperationMethods::EPSG_VERTICAL_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN,
+ CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_2 => CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_2,
+ CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_6 => CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_6,
  ];
 
  /**
@@ -147,6 +149,51 @@ protected static function sign(float $number): int
  * @return array{xt: float, yt: float}
  */
  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,
+ bool $inReverse
+ ): array {
+ if (!$inReverse) {
+ return $this->generalPolynomialUnitlessForward(
+ $xs,
+ $ys,
+ $ordinate1OfEvaluationPointInSourceCRS,
+ $ordinate2OfEvaluationPointInSourceCRS,
+ $ordinate1OfEvaluationPointInTargetCRS,
+ $ordinate2OfEvaluationPointInTargetCRS,
+ $scalingFactorForSourceCRSCoordDifferences,
+ $scalingFactorForTargetCRSCoordDifferences,
+ $A0,
+ $B0,
+ $powerCoefficients,
+ );
+ } else {
+ return $this->generalPolynomialUnitlessReverse(
+ $xs,
+ $ys,
+ $ordinate1OfEvaluationPointInSourceCRS,
+ $ordinate2OfEvaluationPointInSourceCRS,
+ $ordinate1OfEvaluationPointInTargetCRS,
+ $ordinate2OfEvaluationPointInTargetCRS,
+ $scalingFactorForSourceCRSCoordDifferences,
+ $scalingFactorForTargetCRSCoordDifferences,
+ $A0,
+ $B0,
+ $powerCoefficients,
+ );
+ }
+ }
+
+ protected function generalPolynomialUnitlessForward(
  float $xs,
  float $ys,
  UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
@@ -189,6 +236,45 @@ protected function generalPolynomialUnitless(
  return ['xt' => $xt, 'yt' => $yt];
  }
 
+ protected function generalPolynomialUnitlessReverse(
+ float $xs,
+ float $ys,
+ UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
+ UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
+ UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
+ UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
+ Scale $scalingFactorForSourceCRSCoordDifferences,
+ Scale $scalingFactorForTargetCRSCoordDifferences,
+ Scale $A0,
+ Scale $B0,
+ array $powerCoefficients
+ ): array {
+ $result = ['xt' => $xs, 'yt' => $ys];
+ for ($i = 0; $i <= 15; ++$i) {
+ $forwardShiftedCoordinates = $this->generalPolynomialUnitlessForward(
+ $result['xt'],
+ $result['yt'],
+ $ordinate1OfEvaluationPointInSourceCRS,
+ $ordinate2OfEvaluationPointInSourceCRS,
+ $ordinate1OfEvaluationPointInTargetCRS,
+ $ordinate2OfEvaluationPointInTargetCRS,
+ $scalingFactorForSourceCRSCoordDifferences,
+ $scalingFactorForTargetCRSCoordDifferences,
+ $A0,
+ $B0,
+ $powerCoefficients
+ );
+ $deltaError = [
+ 'xt' => $forwardShiftedCoordinates['xt'] - $xs,
+ 'yt' => $forwardShiftedCoordinates['yt'] - $ys,
+ ];
+ $result['xt'] -= $deltaError['xt'];
+ $result['yt'] -= $deltaError['yt'];
+ }
+
+ return $result;
+ }
+
  /**
  * Reversible polynomial.
  * @param array<string, Coefficient> $powerCoefficients

src/Point/ProjectedPoint.php

@@ -2043,7 +2043,8 @@ public function generalPolynomial(
  Scale $scalingFactorForTargetCRSCoordDifferences,
  Scale $A0,
  Scale $B0,
- array $powerCoefficients
+ array $powerCoefficients,
+ bool $inReverse
  ): self {
  $xs = $this->easting->getValue();
  $ys = $this->northing->getValue();
@@ -2059,7 +2060,8 @@ public function generalPolynomial(
  $scalingFactorForTargetCRSCoordDifferences,
  $A0,
  $B0,
- $powerCoefficients
+ $powerCoefficients,
+ $inReverse
  );
 
  $xtUnit = $to->getCoordinateSystem()->getAxes()[0]->getUnitOfMeasureId();

tests/CoordinateOperation/AutoConversionTest.php

@@ -71,10 +71,10 @@ public function testWGS84ToIrishTM(): void
  $toCRS = Projected::fromSRID(Projected::EPSG_TM75_IRISH_GRID);
  $to = $from->convert($toCRS);
 
- self::assertEqualsWithDelta(318977.585, $to->getEasting()->asMetres()->getValue(), 0.01);
- self::assertEqualsWithDelta(311611.131, $to->getNorthing()->asMetres()->getValue(), 0.01);
+ self::assertEqualsWithDelta(318977.534, $to->getEasting()->asMetres()->getValue(), 0.01);
+ self::assertEqualsWithDelta(311611.219, $to->getNorthing()->asMetres()->getValue(), 0.01);
 
- self::assertEquals('J 18978 11611', $to->asGridReferenceWithSpaces(12));
+ self::assertEquals('J 18977 11611', $to->asGridReferenceWithSpaces(12));
  }
 
  public function testETRS89ToIrishTM(): void
@@ -85,10 +85,10 @@ public function testETRS89ToIrishTM(): void
  /** @var IrishGridPoint $to */
  $to = $from->convert($toCRS);
 
- self::assertEqualsWithDelta(318977.528, $to->getEasting()->asMetres()->getValue(), 0.01);
- self::assertEqualsWithDelta(311611.188, $to->getNorthing()->asMetres()->getValue(), 0.01);
+ self::assertEqualsWithDelta(318977.475, $to->getEasting()->asMetres()->getValue(), 0.01);
+ self::assertEqualsWithDelta(311611.277, $to->getNorthing()->asMetres()->getValue(), 0.01);
 
- self::assertEquals('J 18978 11611', $to->asGridReferenceWithSpaces(12));
+ self::assertEquals('J 18977 11611', $to->asGridReferenceWithSpaces(12));
  }
 
  public function testNoop(): void
@@ -540,6 +540,16 @@ public function testLambert93RGF2BCorsica(): void
  self::assertEqualsWithDelta(9.274852923, $to->getLongitude()->asDegrees()->getValue(), 0.0000001);
  }
 
+ public function testIRENET95TM75(): void
+ {
+ $from = GeographicPoint::create(Geographic2D::fromSRID(Geographic2D::EPSG_IRENET95), new Degree(53.349803), new Degree(-6.2628240));
+ $toCRS = Projected::fromSRID(Projected::EPSG_TM75_IRISH_GRID);
+ $to = $from->convert($toCRS);
+
+ self::assertEqualsWithDelta(315732.8042, $to->getEasting()->getValue(), 0.001);
+ self::assertEqualsWithDelta(234667.6783, $to->getNorthing()->getValue(), 0.001);
+ }
+
  public function testRDNAPToITRF2014(): void
  {
  if (!class_exists(GTXETRS89NAPProvider::class)) {

tests/Point/GeographicPointTest.php

@@ -958,7 +958,8 @@ public function testGeneralPolynomial(): void
  'Bu3v2' => new Coefficient(-95.377),
  'Bu2v3' => new Coefficient(-284.294),
  'Bu3v3' => new Coefficient(-853.95),
- ]
+ ],
+ false,
  );
 
  self::assertEqualsWithDelta(55.00002972, $to->getLatitude()->getValue(), 0.00000001);