dvdoug /
PHPCoord
| Total Complexity | 130 |
| Total Lines | 2025 |
| Duplicated Lines | 0 % |
| Changes | 8 | ||
| Bugs | 0 | Features | 0 |
Complex classes like ProjectedPoint often do a lot of different things. To break such a class down, we need to identify a cohesive component within that class. A common approach to find such a component is to look for fields/methods that share the same prefixes, or suffixes.
Once you have determined the fields that belong together, you can apply the Extract Class refactoring. If the component makes sense as a sub-class, Extract Subclass is also a candidate, and is often faster.
While breaking up the class, it is a good idea to analyze how other classes use ProjectedPoint, and based on these observations, apply Extract Interface, too.
| 1 | <?php |
||
| 60 | class ProjectedPoint extends Point implements ConvertiblePoint |
||
| 61 | { |
||
| 62 | use AutoConversion; |
||
| 63 | |||
| 64 | /** |
||
| 65 | * Easting. |
||
| 66 | */ |
||
| 67 | protected Length $easting; |
||
| 68 | |||
| 69 | /** |
||
| 70 | * Northing. |
||
| 71 | */ |
||
| 72 | protected Length $northing; |
||
| 73 | |||
| 74 | /** |
||
| 75 | * Westing. |
||
| 76 | */ |
||
| 77 | protected Length $westing; |
||
| 78 | |||
| 79 | /** |
||
| 80 | * Southing. |
||
| 81 | */ |
||
| 82 | protected Length $southing; |
||
| 83 | |||
| 84 | /** |
||
| 85 | * Coordinate reference system. |
||
| 86 | */ |
||
| 87 | protected Projected $crs; |
||
| 88 | |||
| 89 | /** |
||
| 90 | * Coordinate epoch (date for which the specified coordinates represented this point). |
||
| 91 | */ |
||
| 92 | protected ?DateTimeImmutable $epoch; |
||
| 93 | |||
| 94 | protected function __construct(?Length $easting, ?Length $northing, ?Length $westing, ?Length $southing, Projected $crs, ?DateTimeInterface $epoch = null) |
||
| 95 | { |
||
| 96 | $this->crs = $crs; |
||
| 97 | |||
| 98 | $eastingAxis = $this->getAxisByName(Axis::EASTING); |
||
| 99 | $westingAxis = $this->getAxisByName(Axis::WESTING); |
||
| 100 | $northingAxis = $this->getAxisByName(Axis::NORTHING); |
||
| 101 | $southingAxis = $this->getAxisByName(Axis::SOUTHING); |
||
| 102 | |||
| 103 | if ($easting && $eastingAxis) { |
||
| 104 | $this->easting = Length::convert($easting, $eastingAxis->getUnitOfMeasureId()); |
||
| 105 | $this->westing = $this->easting->multiply(-1); |
||
| 106 | } elseif ($westing && $westingAxis) { |
||
| 107 | $this->westing = Length::convert($westing, $westingAxis->getUnitOfMeasureId()); |
||
| 108 | $this->easting = $this->westing->multiply(-1); |
||
| 109 | } else { |
||
| 110 | throw new InvalidAxesException($crs->getCoordinateSystem()->getAxes()); |
||
| 111 | } |
||
| 112 | |||
| 113 | if ($northing && $northingAxis) { |
||
| 114 | $this->northing = Length::convert($northing, $northingAxis->getUnitOfMeasureId()); |
||
| 115 | $this->southing = $this->northing->multiply(-1); |
||
| 116 | } elseif ($southing && $southingAxis) { |
||
| 117 | $this->southing = Length::convert($southing, $southingAxis->getUnitOfMeasureId()); |
||
| 118 | $this->northing = $this->southing->multiply(-1); |
||
| 119 | } else { |
||
| 120 | throw new InvalidAxesException($crs->getCoordinateSystem()->getAxes()); |
||
| 121 | } |
||
| 122 | |||
| 123 | if ($epoch instanceof DateTime) { |
||
| 124 | $epoch = DateTimeImmutable::createFromMutable($epoch); |
||
| 125 | } |
||
| 126 | $this->epoch = $epoch; |
||
| 127 | } |
||
| 128 | |||
| 129 | public static function create(?Length $easting, ?Length $northing, ?Length $westing, ?Length $southing, Projected $crs, ?DateTimeInterface $epoch = null): self |
||
| 130 | { |
||
| 131 | if ($crs->getSRID() === Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID) { |
||
| 132 | return new BritishNationalGridPoint($easting, $northing, $epoch); |
||
| 133 | } |
||
| 134 | |||
| 135 | if ($crs->getSRID() === Projected::EPSG_TM75_IRISH_GRID) { |
||
| 136 | return new IrishGridPoint($easting, $northing, $epoch); |
||
| 137 | } |
||
| 138 | |||
| 139 | if ($crs->getSRID() === Projected::EPSG_IRENET95_IRISH_TRANSVERSE_MERCATOR) { |
||
| 140 | return new IrishTransverseMercatorPoint($easting, $northing, $epoch); |
||
| 141 | } |
||
| 142 | |||
| 143 | return new static($easting, $northing, $westing, $southing, $crs, $epoch); |
||
| 144 | } |
||
| 145 | |||
| 146 | public static function createFromEastingNorthing(Length $easting, Length $northing, Projected $crs, ?DateTimeInterface $epoch = null): self |
||
| 147 | { |
||
| 148 | return static::create($easting, $northing, null, null, $crs, $epoch); |
||
| 149 | } |
||
| 150 | |||
| 151 | public static function createFromWestingNorthing(Length $westing, Length $northing, Projected $crs, ?DateTimeInterface $epoch = null): self |
||
| 152 | { |
||
| 153 | return static::create(null, $northing, $westing, null, $crs, $epoch); |
||
| 154 | } |
||
| 155 | |||
| 156 | public static function createFromWestingSouthing(Length $westing, Length $southing, Projected $crs, ?DateTimeInterface $epoch = null): self |
||
| 157 | { |
||
| 158 | return static::create(null, null, $westing, $southing, $crs, $epoch); |
||
| 159 | } |
||
| 160 | |||
| 161 | public function getEasting(): Length |
||
| 162 | { |
||
| 163 | return $this->easting; |
||
| 164 | } |
||
| 165 | |||
| 166 | public function getNorthing(): Length |
||
| 167 | { |
||
| 168 | return $this->northing; |
||
| 169 | } |
||
| 170 | |||
| 171 | public function getWesting(): Length |
||
| 172 | { |
||
| 173 | return $this->westing; |
||
| 174 | } |
||
| 175 | |||
| 176 | public function getSouthing(): Length |
||
| 177 | { |
||
| 178 | return $this->southing; |
||
| 179 | } |
||
| 180 | |||
| 181 | public function getCRS(): Projected |
||
| 182 | { |
||
| 183 | return $this->crs; |
||
| 184 | } |
||
| 185 | |||
| 186 | public function getCoordinateEpoch(): ?DateTimeImmutable |
||
| 189 | } |
||
| 190 | |||
| 191 | /** |
||
| 192 | * Calculate distance between two points. |
||
| 193 | * Because this is a simple grid, we can use Pythagoras. |
||
| 194 | */ |
||
| 195 | public function calculateDistance(Point $to): Length |
||
| 196 | { |
||
| 197 | try { |
||
| 198 | if ($to instanceof ConvertiblePoint) { |
||
| 199 | $to = $to->convert($this->crs); |
||
| 200 | } |
||
| 201 | } finally { |
||
| 202 | if ($to->getCRS()->getSRID() !== $this->crs->getSRID()) { |
||
| 203 | throw new InvalidCoordinateReferenceSystemException('Can only calculate distances between two points in the same CRS'); |
||
| 204 | } |
||
| 205 | |||
| 206 | /* @var ProjectedPoint $to */ |
||
| 207 | return new Metre( |
||
| 208 | sqrt( |
||
| 209 | ($to->getEasting()->getValue() - $this->getEasting()->getValue()) ** 2 + |
||
| 210 | ($to->getNorthing()->getValue() - $this->getNorthing()->getValue()) ** 2 |
||
| 211 | ) |
||
| 212 | ); |
||
| 213 | } |
||
| 214 | } |
||
| 215 | |||
| 216 | public function __toString(): string |
||
| 217 | { |
||
| 218 | $values = []; |
||
| 219 | foreach ($this->getCRS()->getCoordinateSystem()->getAxes() as $axis) { |
||
| 220 | if ($axis->getName() === Axis::EASTING) { |
||
| 221 | $values[] = $this->easting; |
||
| 222 | } elseif ($axis->getName() === Axis::NORTHING) { |
||
| 223 | $values[] = $this->northing; |
||
| 224 | } elseif ($axis->getName() === Axis::WESTING) { |
||
| 225 | $values[] = $this->westing; |
||
| 226 | } elseif ($axis->getName() === Axis::SOUTHING) { |
||
| 227 | $values[] = $this->southing; |
||
| 228 | } else { |
||
| 229 | throw new UnknownAxisException(); // @codeCoverageIgnore |
||
| 230 | } |
||
| 231 | } |
||
| 232 | |||
| 233 | return '(' . implode(', ', $values) . ')'; |
||
| 234 | } |
||
| 235 | |||
| 236 | /** |
||
| 237 | * Affine parametric transformation. |
||
| 238 | */ |
||
| 239 | public function affineParametricTransform( |
||
| 240 | Projected $to, |
||
| 241 | Length $A0, |
||
| 242 | Coefficient $A1, |
||
| 243 | Coefficient $A2, |
||
| 244 | Length $B0, |
||
| 245 | Coefficient $B1, |
||
| 246 | Coefficient $B2, |
||
| 247 | bool $inReverse |
||
| 248 | ): self { |
||
| 249 | $xs = $this->easting->getValue(); // native unit to metre conversion already embedded in the scale factor |
||
| 250 | $ys = $this->northing->getValue(); // native unit to metre conversion already embedded in the scale factor |
||
| 251 | |||
| 252 | if ($inReverse) { |
||
| 253 | $D = ($A1->getValue() * $B2->getValue()) - ($A2->getValue() * $B1->getValue()); |
||
| 254 | $a0 = (($A2->getValue() * $B0->asMetres()->getValue()) - ($B2->getValue() * $A0->asMetres()->getValue())) / $D; |
||
| 255 | $b0 = (($B1->getValue() * $A0->asMetres()->getValue()) - ($A1->getValue() * $B0->asMetres()->getValue())) / $D; |
||
| 256 | $a1 = $B2->getValue() / $D; |
||
| 257 | $a2 = -$A2->getValue() / $D; |
||
| 258 | $b1 = -$B1->getValue() / $D; |
||
| 259 | $b2 = $A1->getValue() / $D; |
||
| 260 | } else { |
||
| 261 | $a0 = $A0->asMetres()->getValue(); |
||
| 262 | $a1 = $A1->getValue(); |
||
| 263 | $a2 = $A2->getValue(); |
||
| 264 | $b0 = $B0->asMetres()->getValue(); |
||
| 265 | $b1 = $B1->getValue(); |
||
| 266 | $b2 = $B2->getValue(); |
||
| 267 | } |
||
| 268 | |||
| 269 | $xt = $a0 + ($a1 * $xs) + ($a2 * $ys); |
||
| 270 | $yt = $b0 + ($b1 * $xs) + ($b2 * $ys); |
||
| 271 | |||
| 272 | return static::create(new Metre($xt), new Metre($yt), new Metre(-$xt), new Metre(-$yt), $to, $this->epoch); |
||
| 273 | } |
||
| 274 | |||
| 275 | /** |
||
| 276 | * Albers Equal Area. |
||
| 277 | */ |
||
| 278 | public function albersEqualArea( |
||
| 279 | Geographic $to, |
||
| 280 | Angle $latitudeOfFalseOrigin, |
||
| 281 | Angle $longitudeOfFalseOrigin, |
||
| 282 | Angle $latitudeOf1stStandardParallel, |
||
| 283 | Angle $latitudeOf2ndStandardParallel, |
||
| 284 | Length $eastingAtFalseOrigin, |
||
| 285 | Length $northingAtFalseOrigin |
||
| 286 | ): GeographicPoint { |
||
| 287 | $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue(); |
||
| 288 | $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue(); |
||
| 289 | $phiOrigin = $latitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 290 | $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 291 | $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue(); |
||
| 292 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 293 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 294 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 295 | $e4 = $e ** 4; |
||
| 296 | $e6 = $e ** 6; |
||
| 297 | |||
| 298 | $centralMeridianFirstParallel = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2); |
||
| 299 | $centralMeridianSecondParallel = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2); |
||
| 300 | |||
| 301 | $alphaOrigin = (1 - $e2) * (sin($phiOrigin) / (1 - $e2 * sin($phiOrigin) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phiOrigin)) / (1 + $e * sin($phiOrigin)))); |
||
| 302 | $alphaFirstParallel = (1 - $e2) * (sin($phi1) / (1 - $e2 * sin($phi1) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phi1)) / (1 + $e * sin($phi1)))); |
||
| 303 | $alphaSecondParallel = (1 - $e2) * (sin($phi2) / (1 - $e2 * sin($phi2) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phi2)) / (1 + $e * sin($phi2)))); |
||
| 304 | |||
| 305 | $n = ($centralMeridianFirstParallel ** 2 - $centralMeridianSecondParallel ** 2) / ($alphaSecondParallel - $alphaFirstParallel); |
||
| 306 | $C = $centralMeridianFirstParallel ** 2 + $n * $alphaFirstParallel; |
||
| 307 | $rhoOrigin = $a * sqrt($C - $n * $alphaOrigin) / $n; |
||
| 308 | $rhoPrime = hypot($easting, $rhoOrigin - $northing); |
||
| 309 | $alphaPrime = ($C - $rhoPrime ** 2 * $n ** 2 / $a ** 2) / $n; |
||
| 310 | $betaPrime = self::asin($alphaPrime / (1 - (1 - $e2) / 2 / $e * log((1 - $e) / (1 + $e)))); |
||
| 311 | if ($n > 0) { |
||
| 312 | $theta = atan2($easting, $rhoOrigin - $northing); |
||
| 313 | } else { |
||
| 314 | $theta = atan2(-$easting, $northing - $rhoOrigin); |
||
| 315 | } |
||
| 316 | |||
| 317 | $latitude = $betaPrime + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $betaPrime)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $betaPrime)) + ((761 * $e6 / 45360) * sin(6 * $betaPrime)); |
||
| 318 | $longitude = $longitudeOfFalseOrigin->asRadians()->getValue() + ($theta / $n); |
||
| 319 | |||
| 320 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 321 | } |
||
| 322 | |||
| 323 | /** |
||
| 324 | * American Polyconic. |
||
| 325 | */ |
||
| 326 | public function americanPolyconic( |
||
| 327 | Geographic $to, |
||
| 328 | Angle $latitudeOfNaturalOrigin, |
||
| 329 | Angle $longitudeOfNaturalOrigin, |
||
| 330 | Length $falseEasting, |
||
| 331 | Length $falseNorthing |
||
| 332 | ): GeographicPoint { |
||
| 333 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 334 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 335 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 336 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 337 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 338 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 339 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 340 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 341 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 342 | |||
| 343 | $i = (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256); |
||
| 344 | $ii = (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024); |
||
| 345 | $iii = (15 * $e4 / 256 + 45 * $e6 / 1024); |
||
| 346 | $iv = (35 * $e6 / 3072); |
||
| 347 | |||
| 348 | $MO = $a * ($i * $latitudeOrigin - $ii * sin(2 * $latitudeOrigin) + $iii * sin(4 * $latitudeOrigin) - $iv * sin(6 * $latitudeOrigin)); |
||
| 349 | |||
| 350 | if ($MO === $northing) { |
||
| 351 | $latitude = 0; |
||
| 352 | $longitude = $longitudeOrigin + $easting / $a; |
||
| 353 | } else { |
||
| 354 | $A = ($MO + $northing) / $a; |
||
| 355 | $B = $A ** 2 + $easting ** 2 / $a ** 2; |
||
| 356 | |||
| 357 | $latitude = $A; |
||
| 358 | $C = sqrt(1 - $e2 * sin($latitude) ** 2) * tan($latitude); |
||
| 359 | do { |
||
| 360 | $latitudeN = $latitude; |
||
| 361 | $Ma = $i * $latitude - $ii * sin(2 * $latitude) + $iii * sin(4 * $latitude) - $iv * sin(6 * $latitude); |
||
| 362 | $MnPrime = $i - 2 * $ii * cos(2 * $latitude) + 4 * $iii * cos(4 * $latitude) - 6 * $iv * cos(6 * $latitude); |
||
| 363 | $latitude = $latitude - ($A * ($C * $Ma + 1) - $Ma - $C * ($Ma ** 2 + $B) / 2) / ($e2 * sin(2 * $latitude) * ($Ma ** 2 + $B - 2 * $A * $Ma) / 4 * $C + ($A - $Ma) * ($C * $MnPrime - (2 / sin(2 * $latitude))) - $MnPrime); |
||
| 364 | $C = sqrt(1 - $e2 * sin($latitude) ** 2) * tan($latitude); |
||
| 365 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 366 | |||
| 367 | $longitude = $longitudeOrigin + (self::asin($easting * $C / $a)) / sin($latitude); |
||
| 368 | } |
||
| 369 | |||
| 370 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 371 | } |
||
| 372 | |||
| 373 | /** |
||
| 374 | * Bonne. |
||
| 375 | */ |
||
| 376 | public function bonne( |
||
| 377 | Geographic $to, |
||
| 378 | Angle $latitudeOfNaturalOrigin, |
||
| 379 | Angle $longitudeOfNaturalOrigin, |
||
| 380 | Length $falseEasting, |
||
| 381 | Length $falseNorthing |
||
| 382 | ): GeographicPoint { |
||
| 383 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 384 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 385 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 386 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 387 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 388 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 389 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 390 | $e4 = $e ** 4; |
||
| 391 | $e6 = $e ** 6; |
||
| 392 | |||
| 393 | $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2); |
||
| 394 | $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin)); |
||
| 395 | $rho = hypot($easting, $a * $mO / sin($latitudeOrigin) - $northing) * static::sign($latitudeOrigin); |
||
| 396 | |||
| 397 | $M = $a * $mO / sin($latitudeOrigin) + $MO - $rho; |
||
| 398 | $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256)); |
||
| 399 | $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2)); |
||
| 400 | |||
| 401 | $latitude = $mu + ((3 * $e1 / 2) - (27 * $e1 ** 3 / 32)) * sin(2 * $mu) + ((21 * $e1 ** 2 / 16) - (55 * $e1 ** 4 / 32)) * sin(4 * $mu) + ((151 * $e1 ** 3 / 96)) * sin(6 * $mu) + ((1097 * $e1 ** 4 / 512)) * sin(8 * $mu); |
||
| 402 | $m = cos($latitude) / sqrt(1 - $e2 * sin($latitude) ** 2); |
||
| 403 | |||
| 404 | if ($m === 0.0) { |
||
| 405 | $longitude = $longitudeOrigin; // pole |
||
| 406 | } elseif ($latitudeOrigin >= 0) { |
||
| 407 | $longitude = $longitudeOrigin + $rho * atan2($easting, $a * $mO / sin($latitudeOrigin) - $northing) / $a / $m; |
||
| 408 | } else { |
||
| 409 | $longitude = $longitudeOrigin + $rho * atan2(-$easting, -($a * $mO / sin($latitudeOrigin) - $northing)) / $a / $m; |
||
| 410 | } |
||
| 411 | |||
| 412 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 413 | } |
||
| 414 | |||
| 415 | /** |
||
| 416 | * Bonne South Orientated. |
||
| 417 | */ |
||
| 418 | public function bonneSouthOrientated( |
||
| 419 | Geographic $to, |
||
| 420 | Angle $latitudeOfNaturalOrigin, |
||
| 421 | Angle $longitudeOfNaturalOrigin, |
||
| 422 | Length $falseEasting, |
||
| 423 | Length $falseNorthing |
||
| 424 | ): GeographicPoint { |
||
| 425 | $westing = $falseEasting->asMetres()->getValue() - $this->westing->asMetres()->getValue(); |
||
| 426 | $southing = $falseNorthing->asMetres()->getValue() - $this->southing->asMetres()->getValue(); |
||
| 427 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 428 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 429 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 430 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 431 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 432 | $e4 = $e ** 4; |
||
| 433 | $e6 = $e ** 6; |
||
| 434 | |||
| 435 | $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2); |
||
| 436 | $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin)); |
||
| 437 | $rho = hypot($westing, $a * $mO / sin($latitudeOrigin) - $southing) * static::sign($latitudeOrigin); |
||
| 438 | |||
| 439 | $M = $a * $mO / sin($latitudeOrigin) + $MO - $rho; |
||
| 440 | $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256)); |
||
| 441 | $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2)); |
||
| 442 | |||
| 443 | $latitude = $mu + ((3 * $e1 / 2) - (27 * $e1 ** 3 / 32)) * sin(2 * $mu) + ((21 * $e1 ** 2 / 16) - (55 * $e1 ** 4 / 32)) * sin(4 * $mu) + ((151 * $e1 ** 3 / 96)) * sin(6 * $mu) + ((1097 * $e1 ** 4 / 512)) * sin(8 * $mu); |
||
| 444 | $m = cos($latitude) / sqrt(1 - $e2 * sin($latitude) ** 2); |
||
| 445 | |||
| 446 | if ($m === 0.0) { |
||
| 447 | $longitude = $longitudeOrigin; // pole |
||
| 448 | } elseif ($latitudeOrigin >= 0) { |
||
| 449 | $longitude = $longitudeOrigin + $rho * atan2($westing, $a * $mO / sin($latitudeOrigin) - $southing) / $a / $m; |
||
| 450 | } else { |
||
| 451 | $longitude = $longitudeOrigin + $rho * atan2(-$westing, -($a * $mO / sin($latitudeOrigin) - $southing)) / $a / $m; |
||
| 452 | } |
||
| 453 | |||
| 454 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 455 | } |
||
| 456 | |||
| 457 | /** |
||
| 458 | * Cartesian Grid Offsets |
||
| 459 | * This transformation allows calculation of coordinates in the target system by adding the parameter value to the |
||
| 460 | * coordinate values of the point in the source system. |
||
| 461 | */ |
||
| 462 | public function offsets( |
||
| 463 | Projected $to, |
||
| 464 | Length $eastingOffset, |
||
| 465 | Length $northingOffset |
||
| 466 | ): self { |
||
| 467 | $easting = $this->easting->asMetres()->getValue() + $eastingOffset->asMetres()->getValue(); |
||
| 468 | $northing = $this->northing->asMetres()->getValue() + $northingOffset->asMetres()->getValue(); |
||
| 469 | |||
| 470 | return static::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch); |
||
| 471 | } |
||
| 472 | |||
| 473 | /** |
||
| 474 | * Cassini-Soldner. |
||
| 475 | */ |
||
| 476 | public function cassiniSoldner( |
||
| 477 | Geographic $to, |
||
| 478 | Angle $latitudeOfNaturalOrigin, |
||
| 479 | Angle $longitudeOfNaturalOrigin, |
||
| 480 | Length $falseEasting, |
||
| 481 | Length $falseNorthing |
||
| 482 | ): GeographicPoint { |
||
| 483 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 484 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 485 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 486 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 487 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 488 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 489 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 490 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 491 | |||
| 492 | $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin)); |
||
| 493 | |||
| 494 | $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2)); |
||
| 495 | $M = $MO + $northing; |
||
| 496 | $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256)); |
||
| 497 | $latitudeCentralMeridian = $mu + (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) + (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) + (151 * $e1 ** 3 / 96) * sin(6 * $mu) + (1097 * $e1 ** 4 / 512) * sin(8 * $mu); |
||
| 498 | |||
| 499 | $nu = $a / sqrt((1 - $e2 * sin($latitudeCentralMeridian) ** 2)); |
||
| 500 | $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitudeCentralMeridian) ** 2) ** 1.5; |
||
| 501 | |||
| 502 | $T = tan($latitudeCentralMeridian) ** 2; |
||
| 503 | $D = $easting / $nu; |
||
| 504 | |||
| 505 | $latitude = $latitudeCentralMeridian - ($nu * tan($latitudeCentralMeridian) / $rho) * ($D ** 2 / 2 - (1 + 3 * $T) * $D ** 4 / 24); |
||
| 506 | $longitude = $longitudeOrigin + ($D - $T * $D ** 3 / 3 + (1 + 3 * $T) * $T * $D ** 5 / 15) / cos($latitudeCentralMeridian); |
||
| 507 | |||
| 508 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 509 | } |
||
| 510 | |||
| 511 | /** |
||
| 512 | * Hyperbolic Cassini-Soldner. |
||
| 513 | */ |
||
| 514 | public function hyperbolicCassiniSoldner( |
||
| 515 | Geographic $to, |
||
| 516 | Angle $latitudeOfNaturalOrigin, |
||
| 517 | Angle $longitudeOfNaturalOrigin, |
||
| 518 | Length $falseEasting, |
||
| 519 | Length $falseNorthing |
||
| 520 | ): GeographicPoint { |
||
| 521 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 522 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 523 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 524 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 525 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 526 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 527 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 528 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 529 | |||
| 530 | $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin)); |
||
| 531 | $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2)); |
||
| 532 | |||
| 533 | $latitude1 = $latitudeOrigin + $northing / 1567446; |
||
| 534 | |||
| 535 | $nu = $a / sqrt((1 - $e2 * sin($latitude1) ** 2)); |
||
| 536 | $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitude1) ** 2) ** 1.5; |
||
| 537 | |||
| 538 | $qPrime = $northing ** 3 / (6 * $rho * $nu); |
||
| 539 | $q = ($northing + $qPrime) ** 3 / (6 * $rho * $nu); |
||
| 540 | $M = $MO + $northing + $q; |
||
| 541 | |||
| 542 | $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256)); |
||
| 543 | $latitudeCentralMeridian = $mu + (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) + (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) + (151 * $e1 ** 3 / 96) * sin(6 * $mu) + (1097 * $e1 ** 4 / 512) * sin(8 * $mu); |
||
| 544 | |||
| 545 | $T = tan($latitudeCentralMeridian) ** 2; |
||
| 546 | $D = $easting / $nu; |
||
| 547 | |||
| 548 | $latitude = $latitudeCentralMeridian - ($nu * tan($latitudeCentralMeridian) / $rho) * ($D ** 2 / 2 - (1 + 3 * $T) * $D ** 4 / 24); |
||
| 549 | $longitude = $longitudeOrigin + ($D - $T * $D ** 3 / 3 + (1 + 3 * $T) * $T * $D ** 5 / 15) / cos($latitudeCentralMeridian); |
||
| 550 | |||
| 551 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 552 | } |
||
| 553 | |||
| 554 | /** |
||
| 555 | * Colombia Urban. |
||
| 556 | */ |
||
| 557 | public function columbiaUrban( |
||
| 558 | Geographic $to, |
||
| 559 | Angle $latitudeOfNaturalOrigin, |
||
| 560 | Angle $longitudeOfNaturalOrigin, |
||
| 561 | Length $falseEasting, |
||
| 562 | Length $falseNorthing, |
||
| 563 | Length $projectionPlaneOriginHeight |
||
| 564 | ): GeographicPoint { |
||
| 565 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 566 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 567 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 568 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 569 | $heightOrigin = $projectionPlaneOriginHeight->asMetres()->getValue(); |
||
| 570 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 571 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 572 | |||
| 573 | $rhoOrigin = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** 1.5; |
||
| 574 | |||
| 575 | $nuOrigin = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2)); |
||
| 576 | |||
| 577 | $B = tan($latitudeOrigin) / (2 * $rhoOrigin * $nuOrigin); |
||
| 578 | $C = 1 + $heightOrigin / $a; |
||
| 579 | $D = $rhoOrigin * (1 + $heightOrigin / ($a * (1 - $e2))); |
||
| 580 | |||
| 581 | $latitude = $latitudeOrigin + ($northing / $D) - $B * ($easting / $C) ** 2; |
||
| 582 | $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2)); |
||
| 583 | $longitude = $longitudeOrigin + $easting / ($C * $nu * cos($latitude)); |
||
| 584 | |||
| 585 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 586 | } |
||
| 587 | |||
| 588 | /** |
||
| 589 | * Equal Earth. |
||
| 590 | */ |
||
| 591 | public function equalEarth( |
||
| 592 | Geographic $to, |
||
| 593 | Angle $longitudeOfNaturalOrigin, |
||
| 594 | Length $falseEasting, |
||
| 595 | Length $falseNorthing |
||
| 596 | ): GeographicPoint { |
||
| 597 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 598 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 599 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 600 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 601 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 602 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 603 | $e4 = $e ** 4; |
||
| 604 | $e6 = $e ** 6; |
||
| 605 | |||
| 606 | $qP = (1 - $e2) * ((1 / (1 - $e2)) - (1 / (2 * $e) * log((1 - $e) / (1 + $e)))); |
||
| 607 | $Rq = $a * sqrt($qP / 2); |
||
| 608 | |||
| 609 | $theta = $northing / $Rq; |
||
| 610 | do { |
||
| 611 | $thetaN = $theta; |
||
| 612 | $correctionFactor = ($theta * (1.340264 - 0.081106 * $theta ** 2 + $theta ** 6 * (0.000893 + 0.003796 * $theta ** 2)) - $northing / $Rq) / (1.340264 - 0.243318 * $theta ** 2 + $theta ** 6 * (0.006251 + 0.034164 * $theta ** 2)); |
||
| 613 | $theta = $theta - $correctionFactor; |
||
| 614 | } while (abs($theta - $thetaN) >= static::ITERATION_CONVERGENCE); |
||
| 615 | |||
| 616 | $beta = self::asin(2 * sin($theta) / sqrt(3)); |
||
| 617 | |||
| 618 | $latitude = $beta + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $beta)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $beta)) + ((761 * $e6 / 45360) * sin(6 * $beta)); |
||
| 619 | $longitude = $longitudeOrigin + sqrt(3) * $easting * (1.340264 - 0.243318 * $theta ** 2 + $theta ** 6 * (0.006251 + 0.034164 * $theta ** 2)) / (2 * $Rq * cos($theta)); |
||
| 620 | |||
| 621 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 622 | } |
||
| 623 | |||
| 624 | /** |
||
| 625 | * Equidistant Cylindrical |
||
| 626 | * See method code 1029 for spherical development. See also Pseudo Plate Carree, method code 9825. |
||
| 627 | */ |
||
| 628 | public function equidistantCylindrical( |
||
| 629 | Geographic $to, |
||
| 630 | Angle $latitudeOf1stStandardParallel, |
||
| 631 | Angle $longitudeOfNaturalOrigin, |
||
| 632 | Length $falseEasting, |
||
| 633 | Length $falseNorthing |
||
| 634 | ): GeographicPoint { |
||
| 635 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 636 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 637 | $latitudeFirstParallel = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 638 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 639 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 640 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 641 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 642 | $e4 = $e ** 4; |
||
| 643 | $e6 = $e ** 6; |
||
| 644 | $e8 = $e ** 8; |
||
| 645 | $e10 = $e ** 10; |
||
| 646 | $e12 = $e ** 12; |
||
| 647 | $e14 = $e ** 14; |
||
| 648 | |||
| 649 | $n = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2)); |
||
| 650 | $n2 = $n ** 2; |
||
| 651 | $n3 = $n ** 3; |
||
| 652 | $n4 = $n ** 4; |
||
| 653 | $n5 = $n ** 5; |
||
| 654 | $n6 = $n ** 6; |
||
| 655 | $n7 = $n ** 7; |
||
| 656 | $mu = $northing / ($a * (1 - 1 / 4 * $e2 - 3 / 64 * $e4 - 5 / 256 * $e6 - 175 / 16384 * $e8 - 441 / 65536 * $e10 - 4851 / 1048576 * $e12 - 14157 / 4194304 * $e14)); |
||
| 657 | |||
| 658 | $latitude = $mu + (3 / 2 * $n - 27 / 32 * $n3 + 269 / 512 * $n5 - 6607 / 24576 * $n7) * sin(2 * $mu) |
||
| 659 | + (21 / 16 * $n2 - 55 / 32 * $n4 + 6759 / 4096 * $n6) * sin(4 * $mu) |
||
| 660 | + (151 / 96 * $n3 - 417 / 128 * $n5 + 87963 / 20480 * $n7) * sin(6 * $mu) |
||
| 661 | + (1097 / 512 * $n4 - 15543 / 2560 * $n6) * sin(8 * $mu) |
||
| 662 | + (8011 / 2560 * $n5 - 69119 / 6144 * $n7) * sin(10 * $mu) |
||
| 663 | + (293393 / 61440 * $n6) * sin(12 * $mu) |
||
| 664 | + (6845701 / 860160 * $n7) * sin(14 * $mu); |
||
| 665 | |||
| 666 | $longitude = $longitudeOrigin + $easting * sqrt(1 - $e2 * sin($latitudeFirstParallel) ** 2) / ($a * cos($latitudeFirstParallel)); |
||
| 667 | |||
| 668 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 669 | } |
||
| 670 | |||
| 671 | /** |
||
| 672 | * Guam Projection |
||
| 673 | * Simplified form of Oblique Azimuthal Equidistant projection method. |
||
| 674 | */ |
||
| 675 | public function guamProjection( |
||
| 676 | Geographic $to, |
||
| 677 | Angle $latitudeOfNaturalOrigin, |
||
| 678 | Angle $longitudeOfNaturalOrigin, |
||
| 679 | Length $falseEasting, |
||
| 680 | Length $falseNorthing |
||
| 681 | ): GeographicPoint { |
||
| 682 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 683 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 684 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 685 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 686 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 687 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 688 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 689 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 690 | |||
| 691 | $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin)); |
||
| 692 | $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2)); |
||
| 693 | $i = (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256); |
||
| 694 | |||
| 695 | $latitude = $latitudeOrigin; |
||
| 696 | do { |
||
| 697 | $latitudeN = $latitude; |
||
| 698 | $M = $MO + $northing - ($easting ** 2 * tan($latitude) * sqrt(1 - $e2 * sin($latitude) ** 2) / (2 * $a)); |
||
| 699 | $mu = $M / ($a * $i); |
||
| 700 | $latitude = $mu + (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) + (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) + (151 * $e1 ** 3 / 96) * sin(6 * $mu) + (1097 * $e1 ** 4 / 512) * sin(8 * $mu); |
||
| 701 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 702 | |||
| 703 | $longitude = $longitudeOrigin + $easting * sqrt(1 - $e2 * sin($latitude) ** 2) / ($a * cos($latitude)); |
||
| 704 | |||
| 705 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 706 | } |
||
| 707 | |||
| 708 | /** |
||
| 709 | * Krovak. |
||
| 710 | */ |
||
| 711 | public function krovak( |
||
| 712 | Geographic $to, |
||
| 713 | Angle $latitudeOfProjectionCentre, |
||
| 714 | Angle $longitudeOfOrigin, |
||
| 715 | Angle $coLatitudeOfConeAxis, |
||
| 716 | Angle $latitudeOfPseudoStandardParallel, |
||
| 717 | Scale $scaleFactorOnPseudoStandardParallel, |
||
| 718 | Length $falseEasting, |
||
| 719 | Length $falseNorthing |
||
| 720 | ): GeographicPoint { |
||
| 721 | $longitudeOffset = $this->crs->getDatum()->getPrimeMeridian()->getGreenwichLongitude()->asRadians()->getValue() - $to->getDatum()->getPrimeMeridian()->getGreenwichLongitude()->asRadians()->getValue(); |
||
| 722 | $westing = $this->westing->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 723 | $southing = $this->southing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 724 | $latitudeC = $latitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 725 | $longitudeO = $longitudeOfOrigin->asRadians()->getValue(); |
||
| 726 | $alphaC = $coLatitudeOfConeAxis->asRadians()->getValue(); |
||
| 727 | $latitudeP = $latitudeOfPseudoStandardParallel->asRadians()->getValue(); |
||
| 728 | $kP = $scaleFactorOnPseudoStandardParallel->asUnity()->getValue(); |
||
| 729 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 730 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 731 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 732 | |||
| 733 | $A = $a * sqrt(1 - $e2) / (1 - $e2 * sin($latitudeC) ** 2); |
||
| 734 | $B = sqrt(1 + $e2 * cos($latitudeC) ** 4 / (1 - $e2)); |
||
| 735 | $upsilonO = self::asin(sin($latitudeC) / $B); |
||
| 736 | $tO = tan(M_PI / 4 + $upsilonO / 2) * ((1 + $e * sin($latitudeC)) / (1 - $e * sin($latitudeC))) ** ($e * $B / 2) / (tan(M_PI / 4 + $latitudeC / 2) ** $B); |
||
| 737 | $n = sin($latitudeP); |
||
| 738 | $rO = $kP * $A / tan($latitudeP); |
||
| 739 | |||
| 740 | $r = hypot($southing, $westing) ?: 1; |
||
| 741 | $theta = atan2($westing, $southing); |
||
| 742 | $D = $theta / sin($latitudeP); |
||
| 743 | $T = 2 * (atan(($rO / $r) ** (1 / $n) * tan(M_PI / 4 + $latitudeP / 2)) - M_PI / 4); |
||
| 744 | $U = self::asin(cos($alphaC) * sin($T) - sin($alphaC) * cos($T) * cos($D)); |
||
| 745 | $V = self::asin(cos($T) * sin($D) / cos($U)); |
||
| 746 | |||
| 747 | $latitude = $U; |
||
| 748 | do { |
||
| 749 | $latitudeN = $latitude; |
||
| 750 | $latitude = 2 * (atan($tO ** (-1 / $B) * tan($U / 2 + M_PI / 4) ** (1 / $B) * ((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2)) - M_PI / 4); |
||
| 751 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 752 | |||
| 753 | $longitude = $longitudeO + $longitudeOffset - $V / $B; |
||
| 754 | |||
| 755 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 756 | } |
||
| 757 | |||
| 758 | /** |
||
| 759 | * Krovak Modified |
||
| 760 | * Incorporates a polynomial transformation which is defined to be exact and for practical purposes is considered |
||
| 761 | * to be a map projection. |
||
| 762 | */ |
||
| 763 | public function krovakModified( |
||
| 764 | Geographic $to, |
||
| 765 | Angle $latitudeOfProjectionCentre, |
||
| 766 | Angle $longitudeOfOrigin, |
||
| 767 | Angle $coLatitudeOfConeAxis, |
||
| 768 | Angle $latitudeOfPseudoStandardParallel, |
||
| 769 | Scale $scaleFactorOnPseudoStandardParallel, |
||
| 770 | Length $falseEasting, |
||
| 771 | Length $falseNorthing, |
||
| 772 | Length $ordinate1OfEvaluationPoint, |
||
| 773 | Length $ordinate2OfEvaluationPoint, |
||
| 774 | Coefficient $C1, |
||
| 775 | Coefficient $C2, |
||
| 776 | Coefficient $C3, |
||
| 777 | Coefficient $C4, |
||
| 778 | Coefficient $C5, |
||
| 779 | Coefficient $C6, |
||
| 780 | Coefficient $C7, |
||
| 781 | Coefficient $C8, |
||
| 782 | Coefficient $C9, |
||
| 783 | Coefficient $C10 |
||
| 784 | ): GeographicPoint { |
||
| 785 | $Xr = $this->getSouthing()->asMetres()->getValue() - $falseNorthing->asMetres()->getValue() - $ordinate1OfEvaluationPoint->asMetres()->getValue(); |
||
| 786 | $Yr = $this->getWesting()->asMetres()->getValue() - $falseEasting->asMetres()->getValue() - $ordinate2OfEvaluationPoint->asMetres()->getValue(); |
||
| 787 | $c1 = $C1->asUnity()->getValue(); |
||
| 788 | $c2 = $C2->asUnity()->getValue(); |
||
| 789 | $c3 = $C3->asUnity()->getValue(); |
||
| 790 | $c4 = $C4->asUnity()->getValue(); |
||
| 791 | $c5 = $C5->asUnity()->getValue(); |
||
| 792 | $c6 = $C6->asUnity()->getValue(); |
||
| 793 | $c7 = $C7->asUnity()->getValue(); |
||
| 794 | $c8 = $C8->asUnity()->getValue(); |
||
| 795 | $c9 = $C9->asUnity()->getValue(); |
||
| 796 | $c10 = $C10->asUnity()->getValue(); |
||
| 797 | |||
| 798 | $dX = $c1 + $c3 * $Xr - $c4 * $Yr - 2 * $c6 * $Xr * $Yr + $c5 * ($Xr ** 2 - $Yr ** 2) + $c7 * $Xr * ($Xr ** 2 - 3 * $Yr ** 2) - $c8 * $Yr * (3 * $Xr ** 2 - $Yr ** 2) + 4 * $c9 * $Xr * $Yr * ($Xr ** 2 - $Yr ** 2) + $c10 * ($Xr ** 4 + $Yr ** 4 - 6 * $Xr ** 2 * $Yr ** 2); |
||
| 799 | $dY = $c2 + $c3 * $Yr + $c4 * $Xr + 2 * $c5 * $Xr * $Yr + $c6 * ($Xr ** 2 - $Yr ** 2) + $c8 * $Xr * ($Xr ** 2 - 3 * $Yr ** 2) + $c7 * $Yr * (3 * $Xr ** 2 - $Yr ** 2) - 4 * $c10 * $Xr * $Yr * ($Xr ** 2 - $Yr ** 2) + $c9 * ($Xr ** 4 + $Yr ** 4 - 6 * $Xr ** 2 * $Yr ** 2); |
||
| 800 | |||
| 801 | $Xp = $this->getSouthing()->asMetres()->getValue() - $falseNorthing->asMetres()->getValue() + $dX; |
||
| 802 | $Yp = $this->getWesting()->asMetres()->getValue() - $falseEasting->asMetres()->getValue() + $dY; |
||
| 803 | |||
| 804 | $asKrovak = self::create(new Metre(-$Yp), new Metre(-$Xp), new Metre($Yp), new Metre($Xp), $this->crs, $this->epoch); |
||
| 805 | |||
| 806 | return $asKrovak->krovak($to, $latitudeOfProjectionCentre, $longitudeOfOrigin, $coLatitudeOfConeAxis, $latitudeOfPseudoStandardParallel, $scaleFactorOnPseudoStandardParallel, new Metre(0), new Metre(0)); |
||
| 807 | } |
||
| 808 | |||
| 809 | /** |
||
| 810 | * Lambert Azimuthal Equal Area |
||
| 811 | * This is the ellipsoidal form of the projection. |
||
| 812 | */ |
||
| 813 | public function lambertAzimuthalEqualArea( |
||
| 814 | Geographic $to, |
||
| 815 | Angle $latitudeOfNaturalOrigin, |
||
| 816 | Angle $longitudeOfNaturalOrigin, |
||
| 817 | Length $falseEasting, |
||
| 818 | Length $falseNorthing |
||
| 819 | ): GeographicPoint { |
||
| 820 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 821 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 822 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 823 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 824 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 825 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 826 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 827 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 828 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 829 | |||
| 830 | $qO = (1 - $e2) * ((sin($latitudeOrigin) / (1 - $e2 * sin($latitudeOrigin) ** 2)) - ((1 / (2 * $e)) * log((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))))); |
||
| 831 | $qP = (1 - $e2) * ((1 / (1 - $e2)) - ((1 / (2 * $e)) * log((1 - $e) / (1 + $e)))); |
||
| 832 | $betaO = self::asin($qO / $qP); |
||
| 833 | $Rq = $a * sqrt($qP / 2); |
||
| 834 | $D = $a * (cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2)) / ($Rq * cos($betaO)); |
||
| 835 | $rho = hypot($easting / $D, $D * $northing) ?: 1; |
||
| 836 | $C = 2 * self::asin($rho / (2 * $Rq)); |
||
| 837 | $beta = self::asin(cos($C) * sin($betaO) + ($D * $northing * sin($C) * cos($betaO)) / $rho); |
||
| 838 | |||
| 839 | $latitude = $beta + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $beta)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $beta)) + ((761 * $e6 / 45360) * sin(6 * $beta)); |
||
| 840 | $longitude = $longitudeOrigin + atan2($easting * sin($C), $D * $rho * cos($betaO) * cos($C) - $D ** 2 * $northing * sin($betaO) * sin($C)); |
||
| 841 | |||
| 842 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 843 | } |
||
| 844 | |||
| 845 | /** |
||
| 846 | * Lambert Azimuthal Equal Area (Spherical) |
||
| 847 | * This is the spherical form of the projection. See coordinate operation method Lambert Azimuthal Equal Area |
||
| 848 | * (code 9820) for ellipsoidal form. Differences of several tens of metres result from comparison of the two |
||
| 849 | * methods. |
||
| 850 | */ |
||
| 851 | public function lambertAzimuthalEqualAreaSpherical( |
||
| 852 | Geographic $to, |
||
| 853 | Angle $latitudeOfNaturalOrigin, |
||
| 854 | Angle $longitudeOfNaturalOrigin, |
||
| 855 | Length $falseEasting, |
||
| 856 | Length $falseNorthing |
||
| 857 | ): GeographicPoint { |
||
| 858 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 859 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 860 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 861 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 862 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 863 | |||
| 864 | $rho = hypot($easting, $northing) ?: 1; |
||
| 865 | $c = 2 * self::asin($rho / (2 * $a)); |
||
| 866 | |||
| 867 | $latitude = self::asin(cos($c) * sin($latitudeOrigin) + ($northing * sin($c) * cos($latitudeOrigin) / $rho)); |
||
| 868 | |||
| 869 | if ($latitudeOrigin === 90) { |
||
| 870 | $longitude = $longitudeOrigin + atan($easting / -$northing); |
||
| 871 | } elseif ($latitudeOrigin === -90) { |
||
| 872 | $longitude = $longitudeOrigin + atan($easting / $northing); |
||
| 873 | } else { |
||
| 874 | $longitudeDenominator = ($rho * cos($latitudeOrigin) * cos($c) - $northing * sin($latitudeOrigin) * sin($c)); |
||
| 875 | $longitude = $longitudeOrigin + atan($easting * sin($c) / $longitudeDenominator); |
||
| 876 | if ($longitudeDenominator < 0) { |
||
| 877 | $longitude += M_PI; |
||
| 878 | } |
||
| 879 | } |
||
| 880 | |||
| 881 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 882 | } |
||
| 883 | |||
| 884 | /** |
||
| 885 | * Lambert Conic Conformal (1SP). |
||
| 886 | */ |
||
| 887 | public function lambertConicConformal1SP( |
||
| 888 | Geographic $to, |
||
| 889 | Angle $latitudeOfNaturalOrigin, |
||
| 890 | Angle $longitudeOfNaturalOrigin, |
||
| 891 | Scale $scaleFactorAtNaturalOrigin, |
||
| 892 | Length $falseEasting, |
||
| 893 | Length $falseNorthing |
||
| 894 | ): GeographicPoint { |
||
| 895 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 896 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 897 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 898 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 899 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 900 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 901 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 902 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 903 | |||
| 904 | $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2); |
||
| 905 | $tO = tan(M_PI / 4 - $latitudeOrigin / 2) / ((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))) ** ($e / 2); |
||
| 906 | $n = sin($latitudeOrigin); |
||
| 907 | $F = $mO / ($n * $tO ** $n); |
||
| 908 | $rO = $a * $F * $tO ** $n * $scaleFactorOrigin; |
||
| 909 | $r = hypot($easting, $rO - $northing); |
||
| 910 | if ($n >= 0) { |
||
| 911 | $theta = atan2($easting, $rO - $northing); |
||
| 912 | } else { |
||
| 913 | $r = -$r; |
||
| 914 | $theta = atan2(-$easting, -($rO - $northing)); |
||
| 915 | } |
||
| 916 | |||
| 917 | $t = ($r / ($a * $scaleFactorOrigin * $F)) ** (1 / $n); |
||
| 918 | |||
| 919 | $latitude = M_PI / (2 - 2 * atan($t)); |
||
| 920 | do { |
||
| 921 | $latitudeN = $latitude; |
||
| 922 | $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2)); |
||
| 923 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 924 | |||
| 925 | $longitude = $theta / $n + $longitudeOrigin; |
||
| 926 | |||
| 927 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 928 | } |
||
| 929 | |||
| 930 | /** |
||
| 931 | * Lambert Conic Conformal (west orientated). |
||
| 932 | */ |
||
| 933 | public function lambertConicConformalWestOrientated( |
||
| 934 | Geographic $to, |
||
| 935 | Angle $latitudeOfNaturalOrigin, |
||
| 936 | Angle $longitudeOfNaturalOrigin, |
||
| 937 | Scale $scaleFactorAtNaturalOrigin, |
||
| 938 | Length $falseEasting, |
||
| 939 | Length $falseNorthing |
||
| 940 | ): GeographicPoint { |
||
| 941 | $westing = $falseEasting->asMetres()->getValue() - $this->westing->asMetres()->getValue(); |
||
| 942 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 943 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 944 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 945 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 946 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 947 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 948 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 949 | |||
| 950 | $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2); |
||
| 951 | $tO = tan(M_PI / 4 - $latitudeOrigin / 2) / ((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))) ** ($e / 2); |
||
| 952 | $n = sin($latitudeOrigin); |
||
| 953 | $F = $mO / ($n * $tO ** $n); |
||
| 954 | $rO = $a * $F * $tO ** $n ** $scaleFactorOrigin; |
||
| 955 | $r = hypot($westing, $rO - $northing); |
||
| 956 | if ($n >= 0) { |
||
| 957 | $theta = atan2($westing, $rO - $northing); |
||
| 958 | } else { |
||
| 959 | $r = -$r; |
||
| 960 | $theta = atan2(-$westing, -($rO - $northing)); |
||
| 961 | } |
||
| 962 | |||
| 963 | $t = ($r / ($a * $scaleFactorOrigin * $F)) ** (1 / $n); |
||
| 964 | |||
| 965 | $latitude = M_PI / (2 - 2 * atan($t)); |
||
| 966 | do { |
||
| 967 | $latitudeN = $latitude; |
||
| 968 | $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2)); |
||
| 969 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 970 | |||
| 971 | $longitude = $theta / $n + $longitudeOrigin; |
||
| 972 | |||
| 973 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 974 | } |
||
| 975 | |||
| 976 | /** |
||
| 977 | * Lambert Conic Conformal (1SP) Variant B. |
||
| 978 | */ |
||
| 979 | public function lambertConicConformal1SPVariantB( |
||
| 980 | Geographic $to, |
||
| 981 | Angle $latitudeOfNaturalOrigin, |
||
| 982 | Scale $scaleFactorAtNaturalOrigin, |
||
| 983 | Angle $latitudeOfFalseOrigin, |
||
| 984 | Angle $longitudeOfFalseOrigin, |
||
| 985 | Length $eastingAtFalseOrigin, |
||
| 986 | Length $northingAtFalseOrigin |
||
| 987 | ): GeographicPoint { |
||
| 988 | $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue(); |
||
| 989 | $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue(); |
||
| 990 | $latitudeNaturalOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 991 | $latitudeFalseOrigin = $latitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 992 | $longitudeFalseOrigin = $longitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 993 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 994 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 995 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 996 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 997 | |||
| 998 | $mO = cos($latitudeNaturalOrigin) / sqrt(1 - $e2 * sin($latitudeNaturalOrigin) ** 2); |
||
| 999 | $tO = tan(M_PI / 4 - $latitudeNaturalOrigin / 2) / ((1 - $e * sin($latitudeNaturalOrigin)) / (1 + $e * sin($latitudeNaturalOrigin))) ** ($e / 2); |
||
| 1000 | $tF = tan(M_PI / 4 - $latitudeFalseOrigin / 2) / ((1 - $e * sin($latitudeFalseOrigin)) / (1 + $e * sin($latitudeFalseOrigin))) ** ($e / 2); |
||
| 1001 | $n = sin($latitudeNaturalOrigin); |
||
| 1002 | $F = $mO / ($n * $tO ** $n); |
||
| 1003 | $rF = $a * $F * $tF ** $n * $scaleFactorOrigin; |
||
| 1004 | $r = hypot($easting, $rF - $northing); |
||
| 1005 | if ($n >= 0) { |
||
| 1006 | $theta = atan2($easting, $rF - $northing); |
||
| 1007 | } else { |
||
| 1008 | $r = -$r; |
||
| 1009 | $theta = atan2(-$easting, -($rF - $northing)); |
||
| 1010 | } |
||
| 1011 | |||
| 1012 | $t = ($r / ($a * $scaleFactorOrigin * $F)) ** (1 / $n); |
||
| 1013 | |||
| 1014 | $latitude = M_PI / (2 - 2 * atan($t)); |
||
| 1015 | do { |
||
| 1016 | $latitudeN = $latitude; |
||
| 1017 | $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2)); |
||
| 1018 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1019 | |||
| 1020 | $longitude = $theta / $n + $longitudeFalseOrigin; |
||
| 1021 | |||
| 1022 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1023 | } |
||
| 1024 | |||
| 1025 | /** |
||
| 1026 | * Lambert Conic Conformal (2SP). |
||
| 1027 | */ |
||
| 1028 | public function lambertConicConformal2SP( |
||
| 1029 | Geographic $to, |
||
| 1030 | Angle $latitudeOfFalseOrigin, |
||
| 1031 | Angle $longitudeOfFalseOrigin, |
||
| 1032 | Angle $latitudeOf1stStandardParallel, |
||
| 1033 | Angle $latitudeOf2ndStandardParallel, |
||
| 1034 | Length $eastingAtFalseOrigin, |
||
| 1035 | Length $northingAtFalseOrigin |
||
| 1036 | ): GeographicPoint { |
||
| 1037 | $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue(); |
||
| 1038 | $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue(); |
||
| 1039 | $lambdaOrigin = $longitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 1040 | $phiF = $latitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 1041 | $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 1042 | $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue(); |
||
| 1043 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1044 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1045 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1046 | |||
| 1047 | $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2); |
||
| 1048 | $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2); |
||
| 1049 | $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2); |
||
| 1050 | $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2); |
||
| 1051 | $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2); |
||
| 1052 | $n = (log($m1) - log($m2)) / (log($t1) - log($t2)); |
||
| 1053 | $F = $m1 / ($n * $t1 ** $n); |
||
| 1054 | $rF = $a * $F * $tF ** $n; |
||
| 1055 | $r = hypot($easting, $rF - $northing) * static::sign($n); |
||
| 1056 | $t = ($r / ($a * $F)) ** (1 / $n); |
||
| 1057 | if ($n >= 0) { |
||
| 1058 | $theta = atan2($easting, $rF - $northing); |
||
| 1059 | } else { |
||
| 1060 | $theta = atan2(-$easting, -($rF - $northing)); |
||
| 1061 | } |
||
| 1062 | |||
| 1063 | $latitude = M_PI / 2 - 2 * atan($t); |
||
| 1064 | do { |
||
| 1065 | $latitudeN = $latitude; |
||
| 1066 | $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2)); |
||
| 1067 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1068 | |||
| 1069 | $longitude = $theta / $n + $lambdaOrigin; |
||
| 1070 | |||
| 1071 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1072 | } |
||
| 1073 | |||
| 1074 | /** |
||
| 1075 | * Lambert Conic Conformal (2SP Michigan). |
||
| 1076 | */ |
||
| 1077 | public function lambertConicConformal2SPMichigan( |
||
| 1078 | Geographic $to, |
||
| 1079 | Angle $latitudeOfFalseOrigin, |
||
| 1080 | Angle $longitudeOfFalseOrigin, |
||
| 1081 | Angle $latitudeOf1stStandardParallel, |
||
| 1082 | Angle $latitudeOf2ndStandardParallel, |
||
| 1083 | Length $eastingAtFalseOrigin, |
||
| 1084 | Length $northingAtFalseOrigin, |
||
| 1085 | Scale $ellipsoidScalingFactor |
||
| 1086 | ): GeographicPoint { |
||
| 1087 | $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue(); |
||
| 1088 | $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue(); |
||
| 1089 | $lambdaOrigin = $longitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 1090 | $phiF = $latitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 1091 | $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 1092 | $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue(); |
||
| 1093 | $K = $ellipsoidScalingFactor->asUnity()->getValue(); |
||
| 1094 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1095 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1096 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1097 | |||
| 1098 | $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2); |
||
| 1099 | $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2); |
||
| 1100 | $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2); |
||
| 1101 | $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2); |
||
| 1102 | $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2); |
||
| 1103 | $n = (log($m1) - log($m2)) / (log($t1) - log($t2)); |
||
| 1104 | $F = $m1 / ($n * $t1 ** $n); |
||
| 1105 | $rF = $a * $K * $F * $tF ** $n; |
||
| 1106 | $r = sqrt($easting ** 2 + ($rF - $northing) ** 2) * static::sign($n); |
||
| 1107 | $t = ($r / ($a * $K * $F)) ** (1 / $n); |
||
| 1108 | if ($n >= 0) { |
||
| 1109 | $theta = atan2($easting, $rF - $northing); |
||
| 1110 | } else { |
||
| 1111 | $theta = atan2(-$easting, -($rF - $northing)); |
||
| 1112 | } |
||
| 1113 | |||
| 1114 | $latitude = M_PI / 2 - 2 * atan($t); |
||
| 1115 | do { |
||
| 1116 | $latitudeN = $latitude; |
||
| 1117 | $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2)); |
||
| 1118 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1119 | |||
| 1120 | $longitude = $theta / $n + $lambdaOrigin; |
||
| 1121 | |||
| 1122 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1123 | } |
||
| 1124 | |||
| 1125 | /** |
||
| 1126 | * Lambert Conic Conformal (2SP Belgium) |
||
| 1127 | * In 2000 this modification was replaced through use of the regular Lambert Conic Conformal (2SP) method [9802] |
||
| 1128 | * with appropriately modified parameter values. |
||
| 1129 | */ |
||
| 1130 | public function lambertConicConformal2SPBelgium( |
||
| 1131 | Geographic $to, |
||
| 1132 | Angle $latitudeOfFalseOrigin, |
||
| 1133 | Angle $longitudeOfFalseOrigin, |
||
| 1134 | Angle $latitudeOf1stStandardParallel, |
||
| 1135 | Angle $latitudeOf2ndStandardParallel, |
||
| 1136 | Length $eastingAtFalseOrigin, |
||
| 1137 | Length $northingAtFalseOrigin |
||
| 1138 | ): GeographicPoint { |
||
| 1139 | $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue(); |
||
| 1140 | $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue(); |
||
| 1141 | $lambdaOrigin = $longitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 1142 | $phiF = $latitudeOfFalseOrigin->asRadians()->getValue(); |
||
| 1143 | $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 1144 | $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue(); |
||
| 1145 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1146 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1147 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1148 | |||
| 1149 | $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2); |
||
| 1150 | $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2); |
||
| 1151 | $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2); |
||
| 1152 | $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2); |
||
| 1153 | $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2); |
||
| 1154 | $n = (log($m1) - log($m2)) / (log($t1) - log($t2)); |
||
| 1155 | $F = $m1 / ($n * $t1 ** $n); |
||
| 1156 | $rF = $a * $F * $tF ** $n; |
||
| 1157 | if (is_nan($rF)) { |
||
| 1158 | $rF = 0; |
||
| 1159 | } |
||
| 1160 | $r = hypot($easting, $rF - $northing) * static::sign($n); |
||
| 1161 | $t = ($r / ($a * $F)) ** (1 / $n); |
||
| 1162 | if ($n >= 0) { |
||
| 1163 | $theta = atan2($easting, $rF - $northing); |
||
| 1164 | } else { |
||
| 1165 | $theta = atan2(-$easting, -($rF - $northing)); |
||
| 1166 | } |
||
| 1167 | |||
| 1168 | $latitude = M_PI / 2 - 2 * atan($t); |
||
| 1169 | do { |
||
| 1170 | $latitudeN = $latitude; |
||
| 1171 | $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2)); |
||
| 1172 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1173 | |||
| 1174 | $longitude = ($theta + (new ArcSecond(29.2985))->asRadians()->getValue()) / $n + $lambdaOrigin; |
||
| 1175 | |||
| 1176 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1177 | } |
||
| 1178 | |||
| 1179 | /** |
||
| 1180 | * Lambert Conic Near-Conformal |
||
| 1181 | * The Lambert Near-Conformal projection is derived from the Lambert Conformal Conic projection by truncating the |
||
| 1182 | * series expansion of the projection formulae. |
||
| 1183 | */ |
||
| 1184 | public function lambertConicNearConformal( |
||
| 1185 | Geographic $to, |
||
| 1186 | Angle $latitudeOfNaturalOrigin, |
||
| 1187 | Angle $longitudeOfNaturalOrigin, |
||
| 1188 | Scale $scaleFactorAtNaturalOrigin, |
||
| 1189 | Length $falseEasting, |
||
| 1190 | Length $falseNorthing |
||
| 1191 | ): GeographicPoint { |
||
| 1192 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1193 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1194 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1195 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1196 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 1197 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1198 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1199 | $f = $this->crs->getDatum()->getEllipsoid()->getInverseFlattening(); |
||
| 1200 | |||
| 1201 | $n = $f / (2 - $f); |
||
| 1202 | $rhoO = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2); |
||
| 1203 | $nuO = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2)); |
||
| 1204 | $A = 1 / (6 * $rhoO * $nuO); |
||
| 1205 | $APrime = $a * (1 - $n + 5 * ($n ** 2 - $n ** 3) / 4 + 81 * ($n ** 4 - $n ** 5) / 64); |
||
| 1206 | $BPrime = 3 * $a * ($n - $n ** 2 + 7 * ($n ** 3 - $n ** 4) / 8 + 55 * $n ** 5 / 64) / 2; |
||
| 1207 | $CPrime = 15 * $a * ($n ** 2 - $n ** 3 + 3 * ($n ** 4 - $n ** 5) / 4) / 16; |
||
| 1208 | $DPrime = 35 * $a * ($n ** 3 - $n ** 4 + 11 * $n ** 5 / 16) / 48; |
||
| 1209 | $EPrime = 315 * $a * ($n ** 4 - $n ** 5) / 512; |
||
| 1210 | $rO = $scaleFactorOrigin * $nuO / tan($latitudeOrigin); |
||
| 1211 | $sO = $APrime * $latitudeOrigin - $BPrime * sin(2 * $latitudeOrigin) + $CPrime * sin(4 * $latitudeOrigin) - $DPrime * sin(6 * $latitudeOrigin) + $EPrime * sin(8 * $latitudeOrigin); |
||
| 1212 | |||
| 1213 | $theta = atan2($easting, $rO - $northing); |
||
| 1214 | $r = hypot($easting, $rO - $northing) * static::sign($latitudeOrigin); |
||
| 1215 | $M = $rO - $r; |
||
| 1216 | |||
| 1217 | $m = $M; |
||
| 1218 | do { |
||
| 1219 | $mN = $m; |
||
| 1220 | $m = $m - ($M - $scaleFactorOrigin * $m - $scaleFactorOrigin * $A * $m ** 3) / (-$scaleFactorOrigin - 3 * $scaleFactorOrigin * $A * $m ** 2); |
||
| 1221 | } while (abs($m - $mN) >= static::ITERATION_CONVERGENCE); |
||
| 1222 | |||
| 1223 | $latitude = $latitudeOrigin + $m / $A; |
||
| 1224 | do { |
||
| 1225 | $latitudeN = $latitude; |
||
| 1226 | $latitude = $latitude + ($m + $sO - ($APrime * $latitude - $BPrime * sin(2 * $latitude) + $CPrime * sin(4 * $latitude) - $DPrime * sin(6 * $latitude) + $EPrime * sin(8 * $latitude))) / $APrime; |
||
| 1227 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1228 | |||
| 1229 | $longitude = $longitudeOrigin + $theta / sin($latitudeOrigin); |
||
| 1230 | |||
| 1231 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1232 | } |
||
| 1233 | |||
| 1234 | /** |
||
| 1235 | * Lambert Cylindrical Equal Area |
||
| 1236 | * This is the ellipsoidal form of the projection. |
||
| 1237 | */ |
||
| 1238 | public function lambertCylindricalEqualArea( |
||
| 1239 | Geographic $to, |
||
| 1240 | Angle $latitudeOf1stStandardParallel, |
||
| 1241 | Angle $longitudeOfNaturalOrigin, |
||
| 1242 | Length $falseEasting, |
||
| 1243 | Length $falseNorthing |
||
| 1244 | ): GeographicPoint { |
||
| 1245 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1246 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1247 | $latitudeFirstParallel = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 1248 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1249 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1250 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1251 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1252 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1253 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1254 | |||
| 1255 | $k = cos($latitudeFirstParallel) / sqrt(1 - $e2 * sin($latitudeFirstParallel) ** 2); |
||
| 1256 | $qP = (1 - $e2) * ((sin(M_PI_2) / (1 - $e2 * sin(M_PI_2) ** 2)) - (1 / (2 * $e)) * log((1 - $e * sin(M_PI_2)) / (1 + $e * sin(M_PI_2)))); |
||
| 1257 | $beta = self::asin(2 * $northing * $k / ($a * $qP)); |
||
| 1258 | |||
| 1259 | $latitude = $beta + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $beta)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $beta)) + ((761 * $e6 / 45360) * sin(6 * $beta)); |
||
| 1260 | $longitude = $longitudeOrigin + $easting / ($a * $k); |
||
| 1261 | |||
| 1262 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1263 | } |
||
| 1264 | |||
| 1265 | /** |
||
| 1266 | * Modified Azimuthal Equidistant |
||
| 1267 | * Modified form of Oblique Azimuthal Equidistant projection method developed for Polynesian islands. For the |
||
| 1268 | * distances over which these projections are used (under 800km) this modification introduces no significant error. |
||
| 1269 | */ |
||
| 1270 | public function modifiedAzimuthalEquidistant( |
||
| 1271 | Geographic $to, |
||
| 1272 | Angle $latitudeOfNaturalOrigin, |
||
| 1273 | Angle $longitudeOfNaturalOrigin, |
||
| 1274 | Length $falseEasting, |
||
| 1275 | Length $falseNorthing |
||
| 1276 | ): GeographicPoint { |
||
| 1277 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1278 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1279 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1280 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1281 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1282 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1283 | |||
| 1284 | $nuO = $a / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2); |
||
| 1285 | $c = hypot($easting, $northing); |
||
| 1286 | $alpha = atan2($easting, $northing); |
||
| 1287 | $A = -$e2 * cos($latitudeOrigin) ** 2 * cos($alpha) ** 2 / (1 - $e2); |
||
| 1288 | $B = 3 * $e2 * (1 - $A) * sin($latitudeOrigin) * cos($latitudeOrigin) * cos($alpha) / (1 - $e2); |
||
| 1289 | $D = $c / $nuO; |
||
| 1290 | $J = $D - ($A * (1 + $A) * $D ** 3 / 6) - ($B * (1 + 3 * $A) * $D ** 4 / 24); |
||
| 1291 | $K = 1 - ($A * $J ** 2 / 2) - ($B * $J ** 3 / 6); |
||
| 1292 | $psi = self::asin(sin($latitudeOrigin) * cos($J) + cos($latitudeOrigin) * sin($J) * cos($alpha)); |
||
| 1293 | |||
| 1294 | $latitude = atan((1 - $e2 * $K * sin($latitudeOrigin) / sin($psi)) * tan($psi) / (1 - $e2)); |
||
| 1295 | $longitude = $longitudeOrigin + self::asin(sin($alpha) * sin($J) / cos($psi)); |
||
| 1296 | |||
| 1297 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1298 | } |
||
| 1299 | |||
| 1300 | /** |
||
| 1301 | * Oblique Stereographic |
||
| 1302 | * This is not the same as the projection method of the same name in USGS Professional Paper no. 1395, "Map |
||
| 1303 | * Projections - A Working Manual" by John P. Snyder. |
||
| 1304 | */ |
||
| 1305 | public function obliqueStereographic( |
||
| 1306 | Geographic $to, |
||
| 1307 | Angle $latitudeOfNaturalOrigin, |
||
| 1308 | Angle $longitudeOfNaturalOrigin, |
||
| 1309 | Scale $scaleFactorAtNaturalOrigin, |
||
| 1310 | Length $falseEasting, |
||
| 1311 | Length $falseNorthing |
||
| 1312 | ): GeographicPoint { |
||
| 1313 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1314 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1315 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1316 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1317 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 1318 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1319 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1320 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1321 | |||
| 1322 | $rhoOrigin = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2); |
||
| 1323 | $nuOrigin = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2)); |
||
| 1324 | $R = sqrt($rhoOrigin * $nuOrigin); |
||
| 1325 | |||
| 1326 | $n = sqrt(1 + ($e2 * cos($latitudeOrigin) ** 4 / (1 - $e2))); |
||
| 1327 | $S1 = (1 + sin($latitudeOrigin)) / (1 - sin($latitudeOrigin)); |
||
| 1328 | $S2 = (1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin)); |
||
| 1329 | $w1 = ($S1 * ($S2 ** $e)) ** $n; |
||
| 1330 | $c = (($n + sin($latitudeOrigin)) * (1 - ($w1 - 1) / ($w1 + 1))) / (($n - sin($latitudeOrigin)) * (1 + ($w1 - 1) / ($w1 + 1))); |
||
| 1331 | $w2 = $c * $w1; |
||
| 1332 | $chiOrigin = self::asin(($w2 - 1) / ($w2 + 1)); |
||
| 1333 | |||
| 1334 | $g = 2 * $R * $scaleFactorOrigin * tan(M_PI / 4 - $chiOrigin / 2); |
||
| 1335 | $h = 4 * $R * $scaleFactorOrigin * tan($chiOrigin) + $g; |
||
| 1336 | $i = atan2($easting, ($h + $northing)); |
||
| 1337 | $j = atan2($easting, ($g - $northing)) - $i; |
||
| 1338 | $chi = $chiOrigin + 2 * atan(($northing - $easting * tan($j / 2)) / (2 * $R * $scaleFactorOrigin)); |
||
| 1339 | $lambda = $j + 2 * $i + $longitudeOrigin; |
||
| 1340 | |||
| 1341 | $longitude = ($lambda - $longitudeOrigin) / $n + $longitudeOrigin; |
||
| 1342 | |||
| 1343 | $psi = 0.5 * log((1 + sin($chi)) / ($c * (1 - sin($chi)))) / $n; |
||
| 1344 | |||
| 1345 | $latitude = 2 * atan(M_E ** $psi) - M_PI / 2; |
||
| 1346 | do { |
||
| 1347 | $latitudeN = $latitude; |
||
| 1348 | $psiN = log((tan($latitudeN / 2 + M_PI / 4)) * ((1 - $e * sin($latitudeN)) / (1 + $e * sin($latitudeN))) ** ($e / 2)); |
||
| 1349 | $latitude = $latitudeN - ($psiN - $psi) * cos($latitudeN) * (1 - $e2 * sin($latitudeN) ** 2) / (1 - $e2); |
||
| 1350 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1351 | |||
| 1352 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1353 | } |
||
| 1354 | |||
| 1355 | /** |
||
| 1356 | * Polar Stereographic (variant A) |
||
| 1357 | * Latitude of natural origin must be either 90 degrees or -90 degrees (or equivalent in alternative angle unit). |
||
| 1358 | */ |
||
| 1359 | public function polarStereographicVariantA( |
||
| 1360 | Geographic $to, |
||
| 1361 | Angle $latitudeOfNaturalOrigin, |
||
| 1362 | Angle $longitudeOfNaturalOrigin, |
||
| 1363 | Scale $scaleFactorAtNaturalOrigin, |
||
| 1364 | Length $falseEasting, |
||
| 1365 | Length $falseNorthing |
||
| 1366 | ): GeographicPoint { |
||
| 1367 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1368 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1369 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1370 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1371 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 1372 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1373 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1374 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1375 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1376 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1377 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1378 | |||
| 1379 | $rho = hypot($easting, $northing); |
||
| 1380 | $t = $rho * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $a * $scaleFactorOrigin); |
||
| 1381 | |||
| 1382 | if ($latitudeOrigin < 0) { |
||
| 1383 | $chi = 2 * atan($t) - M_PI / 2; |
||
| 1384 | } else { |
||
| 1385 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1386 | } |
||
| 1387 | |||
| 1388 | $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi); |
||
| 1389 | |||
| 1390 | if ($easting === 0.0) { |
||
| 1391 | $longitude = $longitudeOrigin; |
||
| 1392 | } elseif ($latitudeOrigin < 0) { |
||
| 1393 | $longitude = $longitudeOrigin + atan2($easting, $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue()); |
||
| 1394 | } else { |
||
| 1395 | $longitude = $longitudeOrigin + atan2($easting, $falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue()); |
||
| 1396 | } |
||
| 1397 | |||
| 1398 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1399 | } |
||
| 1400 | |||
| 1401 | /** |
||
| 1402 | * Polar Stereographic (variant B). |
||
| 1403 | */ |
||
| 1404 | public function polarStereographicVariantB( |
||
| 1405 | Geographic $to, |
||
| 1406 | Angle $latitudeOfStandardParallel, |
||
| 1407 | Angle $longitudeOfOrigin, |
||
| 1408 | Length $falseEasting, |
||
| 1409 | Length $falseNorthing |
||
| 1410 | ): GeographicPoint { |
||
| 1411 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1412 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1413 | $standardParallel = $latitudeOfStandardParallel->asRadians()->getValue(); |
||
| 1414 | $longitudeOrigin = $longitudeOfOrigin->asRadians()->getValue(); |
||
| 1415 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1416 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1417 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1418 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1419 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1420 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1421 | |||
| 1422 | $rho = hypot($easting, $northing); |
||
| 1423 | if ($standardParallel < 0) { |
||
| 1424 | $tF = tan(M_PI / 4 + $standardParallel / 2) / (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2)); |
||
| 1425 | } else { |
||
| 1426 | $tF = tan(M_PI / 4 - $standardParallel / 2) * (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2)); |
||
| 1427 | } |
||
| 1428 | $mF = cos($standardParallel) / sqrt(1 - $e2 * sin($standardParallel) ** 2); |
||
| 1429 | $kO = $mF * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $tF); |
||
| 1430 | $t = $rho * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $a * $kO); |
||
| 1431 | |||
| 1432 | if ($standardParallel < 0) { |
||
| 1433 | $chi = 2 * atan($t) - M_PI / 2; |
||
| 1434 | } else { |
||
| 1435 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1436 | } |
||
| 1437 | |||
| 1438 | $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi); |
||
| 1439 | |||
| 1440 | if ($easting === 0.0) { |
||
| 1441 | $longitude = $longitudeOrigin; |
||
| 1442 | } elseif ($standardParallel < 0) { |
||
| 1443 | $longitude = $longitudeOrigin + atan2($easting, $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue()); |
||
| 1444 | } else { |
||
| 1445 | $longitude = $longitudeOrigin + atan2($easting, $falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue()); |
||
| 1446 | } |
||
| 1447 | |||
| 1448 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1449 | } |
||
| 1450 | |||
| 1451 | /** |
||
| 1452 | * Polar Stereographic (variant C). |
||
| 1453 | */ |
||
| 1454 | public function polarStereographicVariantC( |
||
| 1455 | Geographic $to, |
||
| 1456 | Angle $latitudeOfStandardParallel, |
||
| 1457 | Angle $longitudeOfOrigin, |
||
| 1458 | Length $eastingAtFalseOrigin, |
||
| 1459 | Length $northingAtFalseOrigin |
||
| 1460 | ): GeographicPoint { |
||
| 1461 | $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue(); |
||
| 1462 | $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue(); |
||
| 1463 | $standardParallel = $latitudeOfStandardParallel->asRadians()->getValue(); |
||
| 1464 | $longitudeOrigin = $longitudeOfOrigin->asRadians()->getValue(); |
||
| 1465 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1466 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1467 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1468 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1469 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1470 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1471 | |||
| 1472 | if ($standardParallel < 0) { |
||
| 1473 | $tF = tan(M_PI / 4 + $standardParallel / 2) / (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2)); |
||
| 1474 | } else { |
||
| 1475 | $tF = tan(M_PI / 4 - $standardParallel / 2) * (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2)); |
||
| 1476 | } |
||
| 1477 | $mF = cos($standardParallel) / sqrt(1 - $e2 * sin($standardParallel) ** 2); |
||
| 1478 | $rhoF = $a * $mF; |
||
| 1479 | if ($standardParallel < 0) { |
||
| 1480 | $rho = hypot($easting, $northing + $rhoF); |
||
| 1481 | $t = $rho * $tF / $rhoF; |
||
| 1482 | $chi = 2 * atan($t) - M_PI / 2; |
||
| 1483 | } else { |
||
| 1484 | $rho = hypot($easting, $northing - $rhoF); |
||
| 1485 | $t = $rho * $tF / $rhoF; |
||
| 1486 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1487 | } |
||
| 1488 | |||
| 1489 | $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi); |
||
| 1490 | |||
| 1491 | if ($easting === 0.0) { |
||
| 1492 | $longitude = $longitudeOrigin; |
||
| 1493 | } elseif ($standardParallel < 0) { |
||
| 1494 | $longitude = $longitudeOrigin + atan2($easting, $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue() + $rhoF); |
||
| 1495 | } else { |
||
| 1496 | $longitude = $longitudeOrigin + atan2($easting, $northingAtFalseOrigin->asMetres()->getValue() - $this->northing->asMetres()->getValue() + $rhoF); |
||
| 1497 | } |
||
| 1498 | |||
| 1499 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1500 | } |
||
| 1501 | |||
| 1502 | /** |
||
| 1503 | * Popular Visualisation Pseudo Mercator |
||
| 1504 | * Applies spherical formulas to the ellipsoid. As such does not have the properties of a true Mercator projection. |
||
| 1505 | */ |
||
| 1506 | public function popularVisualisationPseudoMercator( |
||
| 1524 | } |
||
| 1525 | |||
| 1526 | /** |
||
| 1527 | * Similarity transformation |
||
| 1528 | * Defined for two-dimensional coordinate systems. |
||
| 1529 | */ |
||
| 1530 | public function similarityTransformation( |
||
| 1531 | Projected $to, |
||
| 1532 | Length $ordinate1OfEvaluationPointInTargetCRS, |
||
| 1533 | Length $ordinate2OfEvaluationPointInTargetCRS, |
||
| 1534 | Scale $scaleFactorForSourceCRSAxes, |
||
| 1535 | Angle $rotationAngleOfSourceCRSAxes, |
||
| 1536 | bool $inReverse |
||
| 1537 | ): self { |
||
| 1538 | $xs = $this->easting->asMetres()->getValue(); |
||
| 1539 | $ys = $this->northing->asMetres()->getValue(); |
||
| 1540 | $xo = $ordinate1OfEvaluationPointInTargetCRS->asMetres()->getValue(); |
||
| 1541 | $yo = $ordinate2OfEvaluationPointInTargetCRS->asMetres()->getValue(); |
||
| 1542 | $M = $scaleFactorForSourceCRSAxes->asUnity()->getValue(); |
||
| 1543 | $theta = $rotationAngleOfSourceCRSAxes->asRadians()->getValue(); |
||
| 1544 | |||
| 1545 | if ($inReverse) { |
||
| 1546 | $easting = (($xs - $xo) * cos($theta) - ($ys - $yo) * sin($theta)) / $M; |
||
| 1547 | $northing = (($xs - $xo) * sin($theta) + ($ys - $yo) * cos($theta)) / $M; |
||
| 1548 | } else { |
||
| 1549 | $easting = $xo + $xs * $M * cos($theta) + $ys * $M * sin($theta); |
||
| 1550 | $northing = $yo - $xs * $M * sin($theta) + $ys * $M * cos($theta); |
||
| 1551 | } |
||
| 1552 | |||
| 1553 | return self::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch); |
||
| 1554 | } |
||
| 1555 | |||
| 1556 | /** |
||
| 1557 | * Mercator (variant A) |
||
| 1558 | * Note that in these formulas the parameter latitude of natural origin (latO) is not used. However for this |
||
| 1559 | * Mercator (variant A) method the EPSG dataset includes this parameter, which must have a value of zero, for |
||
| 1560 | * completeness in CRS labelling. |
||
| 1561 | */ |
||
| 1562 | public function mercatorVariantA( |
||
| 1563 | Geographic $to, |
||
| 1564 | Angle $latitudeOfNaturalOrigin, |
||
| 1565 | Angle $longitudeOfNaturalOrigin, |
||
| 1566 | Scale $scaleFactorAtNaturalOrigin, |
||
| 1567 | Length $falseEasting, |
||
| 1568 | Length $falseNorthing |
||
| 1569 | ): GeographicPoint { |
||
| 1570 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1571 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1572 | $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1573 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1574 | $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue(); |
||
| 1575 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1576 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1577 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1578 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1579 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1580 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1581 | |||
| 1582 | $t = M_E ** (($falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue()) / ($a * $scaleFactorOrigin)); |
||
| 1583 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1584 | |||
| 1585 | $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi); |
||
| 1586 | $longitude = $easting / ($a * $scaleFactorOrigin) + $longitudeOrigin; |
||
| 1587 | |||
| 1588 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1589 | } |
||
| 1590 | |||
| 1591 | /** |
||
| 1592 | * Mercator (variant B) |
||
| 1593 | * Used for most nautical charts. |
||
| 1594 | */ |
||
| 1595 | public function mercatorVariantB( |
||
| 1596 | Geographic $to, |
||
| 1597 | Angle $latitudeOf1stStandardParallel, |
||
| 1598 | Angle $longitudeOfNaturalOrigin, |
||
| 1599 | Length $falseEasting, |
||
| 1600 | Length $falseNorthing |
||
| 1601 | ): GeographicPoint { |
||
| 1602 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1603 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1604 | $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue(); |
||
| 1605 | $firstStandardParallel = $latitudeOf1stStandardParallel->asRadians()->getValue(); |
||
| 1606 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1607 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1608 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1609 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1610 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1611 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1612 | |||
| 1613 | $scaleFactorOrigin = cos($firstStandardParallel) / sqrt(1 - $e2 * sin($firstStandardParallel) ** 2); |
||
| 1614 | |||
| 1615 | $t = M_E ** (($falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue()) / ($a * $scaleFactorOrigin)); |
||
| 1616 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1617 | |||
| 1618 | $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi); |
||
| 1619 | $longitude = $easting / ($a * $scaleFactorOrigin) + $longitudeOrigin; |
||
| 1620 | |||
| 1621 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1622 | } |
||
| 1623 | |||
| 1624 | /** |
||
| 1625 | * Hotine Oblique Mercator (variant A). |
||
| 1626 | */ |
||
| 1627 | public function obliqueMercatorHotineVariantA( |
||
| 1628 | Geographic $to, |
||
| 1629 | Angle $latitudeOfProjectionCentre, |
||
| 1630 | Angle $longitudeOfProjectionCentre, |
||
| 1631 | Angle $azimuthOfInitialLine, |
||
| 1632 | Angle $angleFromRectifiedToSkewGrid, |
||
| 1633 | Scale $scaleFactorOnInitialLine, |
||
| 1634 | Length $falseEasting, |
||
| 1635 | Length $falseNorthing |
||
| 1636 | ): GeographicPoint { |
||
| 1637 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1638 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1639 | $latC = $latitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 1640 | $lonC = $longitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 1641 | $alphaC = $azimuthOfInitialLine->asRadians()->getValue(); |
||
| 1642 | $kC = $scaleFactorOnInitialLine->asUnity()->getValue(); |
||
| 1643 | $gammaC = $angleFromRectifiedToSkewGrid->asRadians()->getValue(); |
||
| 1644 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1645 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1646 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1647 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1648 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1649 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1650 | |||
| 1651 | $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2))); |
||
| 1652 | $A = $a * $B * $kC * sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2); |
||
| 1653 | $tO = tan(M_PI / 4 - $latC / 2) / ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2); |
||
| 1654 | $D = $B * sqrt((1 - $e2)) / (cos($latC) * sqrt(1 - $e2 * sin($latC) ** 2)); |
||
| 1655 | $DD = max(1, $D ** 2); |
||
| 1656 | $F = $D + sqrt($DD - 1) * static::sign($latC); |
||
| 1657 | $H = $F * ($tO) ** $B; |
||
| 1658 | $G = ($F - 1 / $F) / 2; |
||
| 1659 | $gammaO = self::asin(sin($alphaC) / $D); |
||
| 1660 | $lonO = $lonC - (self::asin($G * tan($gammaO))) / $B; |
||
| 1661 | |||
| 1662 | $v = $easting * cos($gammaC) - $northing * sin($gammaC); |
||
| 1663 | $u = $northing * cos($gammaC) + $easting * sin($gammaC); |
||
| 1664 | |||
| 1665 | $Q = M_E ** -($B * $v / $A); |
||
| 1666 | $S = ($Q - 1 / $Q) / 2; |
||
| 1667 | $T = ($Q + 1 / $Q) / 2; |
||
| 1668 | $V = sin($B * $u / $A); |
||
| 1669 | $U = ($V * cos($gammaO) + $S * sin($gammaO)) / $T; |
||
| 1670 | $t = ($H / sqrt((1 + $U) / (1 - $U))) ** (1 / $B); |
||
| 1671 | |||
| 1672 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1673 | |||
| 1674 | $latitude = $chi + sin(2 * $chi) * ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) + sin(4 * $chi) * (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) + sin(6 * $chi) * (7 * $e6 / 120 + 81 * $e8 / 1120) + sin(8 * $chi) * (4279 * $e8 / 161280); |
||
| 1675 | $longitude = $lonO - atan2(($S * cos($gammaO) - $V * sin($gammaO)), cos($B * $u / $A)) / $B; |
||
| 1676 | |||
| 1677 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1678 | } |
||
| 1679 | |||
| 1680 | /** |
||
| 1681 | * Hotine Oblique Mercator (variant B). |
||
| 1682 | */ |
||
| 1683 | public function obliqueMercatorHotineVariantB( |
||
| 1684 | Geographic $to, |
||
| 1685 | Angle $latitudeOfProjectionCentre, |
||
| 1686 | Angle $longitudeOfProjectionCentre, |
||
| 1687 | Angle $azimuthOfInitialLine, |
||
| 1688 | Angle $angleFromRectifiedToSkewGrid, |
||
| 1689 | Scale $scaleFactorOnInitialLine, |
||
| 1690 | Length $eastingAtProjectionCentre, |
||
| 1691 | Length $northingAtProjectionCentre |
||
| 1692 | ): GeographicPoint { |
||
| 1693 | $easting = $this->easting->asMetres()->getValue() - $eastingAtProjectionCentre->asMetres()->getValue(); |
||
| 1694 | $northing = $this->northing->asMetres()->getValue() - $northingAtProjectionCentre->asMetres()->getValue(); |
||
| 1695 | $latC = $latitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 1696 | $lonC = $longitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 1697 | $alphaC = $azimuthOfInitialLine->asRadians()->getValue(); |
||
| 1698 | $kC = $scaleFactorOnInitialLine->asUnity()->getValue(); |
||
| 1699 | $gammaC = $angleFromRectifiedToSkewGrid->asRadians()->getValue(); |
||
| 1700 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1701 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1702 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1703 | $e4 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 4; |
||
| 1704 | $e6 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 6; |
||
| 1705 | $e8 = $this->crs->getDatum()->getEllipsoid()->getEccentricity() ** 8; |
||
| 1706 | |||
| 1707 | $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2))); |
||
| 1708 | $A = $a * $B * $kC * sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2); |
||
| 1709 | $tO = tan(M_PI / 4 - $latC / 2) / ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2); |
||
| 1710 | $D = $B * sqrt((1 - $e2)) / (cos($latC) * sqrt(1 - $e2 * sin($latC) ** 2)); |
||
| 1711 | $DD = max(1, $D ** 2); |
||
| 1712 | $F = $D + sqrt($DD - 1) * static::sign($latC); |
||
| 1713 | $H = $F * ($tO) ** $B; |
||
| 1714 | $G = ($F - 1 / $F) / 2; |
||
| 1715 | $gammaO = self::asin(sin($alphaC) / $D); |
||
| 1716 | $lonO = $lonC - (self::asin($G * tan($gammaO))) / $B; |
||
| 1717 | $vC = 0; |
||
| 1718 | if ($alphaC === M_PI / 2) { |
||
| 1719 | $uC = $A * ($lonC - $lonO); |
||
| 1720 | } else { |
||
| 1721 | $uC = ($A / $B) * atan2(sqrt($DD - 1), cos($alphaC)) * static::sign($latC); |
||
| 1722 | } |
||
| 1723 | |||
| 1724 | $v = $easting * cos($gammaC) - $northing * sin($gammaC); |
||
| 1725 | $u = $northing * cos($gammaC) + $easting * sin($gammaC) + (abs($uC) * static::sign($latC)); |
||
| 1726 | |||
| 1727 | $Q = M_E ** -($B * $v / $A); |
||
| 1728 | $S = ($Q - 1 / $Q) / 2; |
||
| 1729 | $T = ($Q + 1 / $Q) / 2; |
||
| 1730 | $V = sin($B * $u / $A); |
||
| 1731 | $U = ($V * cos($gammaO) + $S * sin($gammaO)) / $T; |
||
| 1732 | $t = ($H / sqrt((1 + $U) / (1 - $U))) ** (1 / $B); |
||
| 1733 | |||
| 1734 | $chi = M_PI / 2 - 2 * atan($t); |
||
| 1735 | |||
| 1736 | $latitude = $chi + sin(2 * $chi) * ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) + sin(4 * $chi) * (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) + sin(6 * $chi) * (7 * $e6 / 120 + 81 * $e8 / 1120) + sin(8 * $chi) * (4279 * $e8 / 161280); |
||
| 1737 | $longitude = $lonO - atan2(($S * cos($gammaO) - $V * sin($gammaO)), cos($B * $u / $A)) / $B; |
||
| 1738 | |||
| 1739 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1740 | } |
||
| 1741 | |||
| 1742 | /** |
||
| 1743 | * Laborde Oblique Mercator. |
||
| 1744 | */ |
||
| 1745 | public function obliqueMercatorLaborde( |
||
| 1746 | Geographic $to, |
||
| 1747 | Angle $latitudeOfProjectionCentre, |
||
| 1748 | Angle $longitudeOfProjectionCentre, |
||
| 1749 | Angle $azimuthOfInitialLine, |
||
| 1750 | Scale $scaleFactorOnInitialLine, |
||
| 1751 | Length $falseEasting, |
||
| 1752 | Length $falseNorthing |
||
| 1753 | ): GeographicPoint { |
||
| 1754 | $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue(); |
||
| 1755 | $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue(); |
||
| 1756 | $latC = $latitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 1757 | $lonC = $longitudeOfProjectionCentre->asRadians()->getValue(); |
||
| 1758 | $alphaC = $azimuthOfInitialLine->asRadians()->getValue(); |
||
| 1759 | $kC = $scaleFactorOnInitialLine->asUnity()->getValue(); |
||
| 1760 | $a = $this->crs->getDatum()->getEllipsoid()->getSemiMajorAxis()->asMetres()->getValue(); |
||
| 1761 | $e = $this->crs->getDatum()->getEllipsoid()->getEccentricity(); |
||
| 1762 | $e2 = $this->crs->getDatum()->getEllipsoid()->getEccentricitySquared(); |
||
| 1763 | |||
| 1764 | $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2))); |
||
| 1765 | $latS = self::asin(sin($latC) / $B); |
||
| 1766 | $R = $a * $kC * (sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2)); |
||
| 1767 | $C = log(tan(M_PI / 4 + $latS / 2)) - $B * log(tan(M_PI / 4 + $latC / 2) * ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2)); |
||
| 1768 | |||
| 1769 | $G = (new ComplexNumber(1 - cos(2 * $alphaC), sin(2 * $alphaC)))->divide(new ComplexNumber(12, 0)); |
||
| 1770 | |||
| 1771 | $H0 = new ComplexNumber($northing / $R, $easting / $R); |
||
| 1772 | $H = $H0->divide($H0->pow(3)->multiply($G)->add($H0)); |
||
| 1773 | do { |
||
| 1774 | $HN = $H; |
||
| 1775 | $H = ($HN->pow(3)->multiply($G)->multiply(new ComplexNumber(2, 0))->add($H0))->divide($HN->pow(2)->multiply($G)->multiply(new ComplexNumber(3, 0))->add(new ComplexNumber(1, 0))); |
||
| 1776 | } while (abs($H0->subtract($H)->subtract($H->pow(3)->multiply($G))->getReal()) >= static::ITERATION_CONVERGENCE); |
||
| 1777 | |||
| 1778 | $LPrime = -1 * $H->getReal(); |
||
| 1779 | $PPrime = 2 * atan(M_E ** $H->getImaginary()) - M_PI / 2; |
||
| 1780 | $U = cos($PPrime) * cos($LPrime) * cos($latS) + cos($PPrime) * sin($LPrime) * sin($latS); |
||
| 1781 | $V = sin($PPrime); |
||
| 1782 | $W = cos($PPrime) * cos($LPrime) * sin($latS) - cos($PPrime) * sin($LPrime) * cos($latS); |
||
| 1783 | |||
| 1784 | $d = hypot($U, $V); |
||
| 1785 | if ($d === 0) { |
||
| 1786 | $L = 0; |
||
| 1787 | $P = static::sign($W) * M_PI / 2; |
||
| 1788 | } else { |
||
| 1789 | $L = 2 * atan($V / ($U + $d)); |
||
| 1790 | $P = atan($W / $d); |
||
| 1791 | } |
||
| 1792 | |||
| 1793 | $longitude = $lonC + ($L / $B); |
||
| 1794 | |||
| 1795 | $q = (log(tan(M_PI / 4 + $P / 2)) - $C) / $B; |
||
| 1796 | |||
| 1797 | $latitude = 2 * atan(M_E ** $q) - M_PI / 2; |
||
| 1798 | do { |
||
| 1799 | $latitudeN = $latitude; |
||
| 1800 | $latitude = 2 * atan(((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2) * M_E ** $q) - M_PI / 2; |
||
| 1801 | } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE); |
||
| 1802 | |||
| 1803 | return GeographicPoint::create(new Radian($latitude), new Radian($longitude), null, $to, $this->epoch); |
||
| 1804 | } |
||
| 1805 | |||
| 1806 | /** |
||
| 1807 | * Transverse Mercator. |
||
| 1808 | */ |
||
| 1809 | public function transverseMercator( |
||
| 1878 | } |
||
| 1879 | |||
| 1880 | /** |
||
| 1881 | * Transverse Mercator Zoned Grid System |
||
| 1882 | * If locations fall outwith the fixed zones the general Transverse Mercator method (code 9807) must be used for |
||
| 1883 | * each zone. |
||
| 1884 | */ |
||
| 1885 | public function transverseMercatorZonedGrid( |
||
| 1886 | Geographic $to, |
||
| 1887 | Angle $latitudeOfNaturalOrigin, |
||
| 1888 | Angle $initialLongitude, |
||
| 1889 | Angle $zoneWidth, |
||
| 1890 | Scale $scaleFactorAtNaturalOrigin, |
||
| 1891 | Length $falseEasting, |
||
| 1892 | Length $falseNorthing |
||
| 1893 | ): GeographicPoint { |
||
| 1894 | $Z = substr((string) $this->easting->asMetres()->getValue(), 0, 2); |
||
| 1895 | $falseEasting = $falseEasting->add(new Metre($Z * 1000000)); |
||
| 1896 | |||
| 1897 | $W = $zoneWidth->asDegrees()->getValue(); |
||
| 1898 | $longitudeOrigin = $initialLongitude->add(new Degree($Z * $W - $W / 2)); |
||
| 1899 | |||
| 1900 | return $this->transverseMercator($to, $latitudeOfNaturalOrigin, $longitudeOrigin, $scaleFactorAtNaturalOrigin, $falseEasting, $falseNorthing); |
||
| 1901 | } |
||
| 1902 | |||
| 1903 | /** |
||
| 1904 | * General polynomial. |
||
| 1905 | * @param Coefficient[] $powerCoefficients |
||
| 1906 | */ |
||
| 1907 | public function generalPolynomial( |
||
| 1944 | ); |
||
| 1945 | } |
||
| 1946 | |||
| 1947 | /** |
||
| 1948 | * New Zealand Map Grid. |
||
| 1949 | */ |
||
| 1950 | public function newZealandMapGrid( |
||
| 2019 | } |
||
| 2020 | |||
| 2021 | /** |
||
| 2022 | * Complex polynomial. |
||
| 2023 | * Coordinate pairs treated as complex numbers. This exploits the correlation between the polynomial coefficients |
||
| 2024 | * and leads to a smaller number of coefficients than the general polynomials. |
||
| 2025 | */ |
||
| 2026 | public function complexPolynomial( |
||
| 2027 | Projected $to, |
||
| 2028 | Length $ordinate1OfEvaluationPointInSourceCRS, |
||
| 2029 | Length $ordinate2OfEvaluationPointInSourceCRS, |
||
| 2030 | Length $ordinate1OfEvaluationPointInTargetCRS, |
||
| 2031 | Length $ordinate2OfEvaluationPointInTargetCRS, |
||
| 2032 | Scale $scalingFactorForSourceCRSCoordDifferences, |
||
| 2033 | Scale $scalingFactorForTargetCRSCoordDifferences, |
||
| 2034 | Scale $A1, |
||
| 2035 | Scale $A2, |
||
| 2036 | Scale $A3, |
||
| 2037 | Scale $A4, |
||
| 2038 | Scale $A5, |
||
| 2039 | Scale $A6, |
||
| 2070 | ); |
||
| 2071 | } |
||
| 2072 | |||
| 2073 | /** |
||
| 2074 | * Ordnance Survey National Transformation |
||
| 2075 | * Geodetic transformation between ETRS89 (or WGS 84) and OSGB36 / National Grid. Uses ETRS89 / National Grid as |
||
| 2076 | * an intermediate coordinate system for bi-linear interpolation of gridded grid coordinate differences. |
||
| 2077 | */ |
||
| 2078 | public function OSTN15( |
||
| 2079 | Geographic2D $to, |
||
| 2080 | OSTNOSGM15Grid $eastingAndNorthingDifferenceFile |
||
| 2081 | ): GeographicPoint { |
||
| 2082 | $asETRS89 = $eastingAndNorthingDifferenceFile->applyReverseAdjustment($this); |
||
| 2083 | |||
| 2084 | return $asETRS89->transverseMercator($to, new Degree(49), new Degree(-2), new Unity(0.9996012717), new Metre(400000), new Metre(-100000)); |
||
| 2087 |
Loading history...
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:For further information see https://scrutinizer-ci.com/docs/tools/php/php-scrutinizer/#list-dependency-paths