Ja, man kan ju tycka att det borde räcka med en 3×3-matris för att konvertera mellan olika koordinatsystem, men matten blir ofta mer avancerad än så. Då kart-länkarna till hitta.se byggde på RT90 var jag tvungen att konvertera WGS84-koordinaterna i databasen till RT90 för att stödja kartlänkar till hitta.se. Den delen av phpoi blev såhär:
<?php
/******************************************************************
Copyright 2008 Henrik Carlqvist
This file is part of phpoi.
Phpoi is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Phpoi is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with phpoi. If not, see <http://www.gnu.org/licenses/>.
******************************************************************/
function RT90($lat, $long)
{
$sin_lat = sin(deg2rad($lat));
$cos_lat = cos(deg2rad($lat));
$sin_lon = sin(deg2rad($long));
$cos_lon = cos(deg2rad($long));
$a_wgs = 6378137.0;
$e2p_wgs = 0.00669437999013;
$Bessel_Ellipsoid_a = 6377397.1542;
$Bessel_Ellipsoid_e2p = 0.006674372232;
$Gd_dR0 = 424.3;
$Gd_dR1 = -80.5;
$Gd_dR2 = 613.1;
$Gd_dC0 = 0.000021315;
$Gd_dC1 = -0.0000096313;
$Gd_dC2 = 0.000025136;
/* Transform WGS84 Lat/Long/H to geo-centric coordinates */
$N_wgs = $a_wgs/sqrt( 1.0 - $e2p_wgs * $sin_lat*$sin_lat );
$X_wgs = $N_wgs*$cos_lat*$cos_lon;
$Y_wgs = $N_wgs*$cos_lat*$sin_lon;
$Z_wgs = $N_wgs*(1.0 - $e2p_wgs)*$sin_lat;
/* Bursa-Wolf transformation */
$X_gd = ( $X_wgs-$Gd_dC2*$Y_wgs+$Gd_dC1*$Z_wgs ) - $Gd_dR0;
$Y_gd = ( $Y_wgs+$Gd_dC2*$X_wgs-$Gd_dC0*$Z_wgs ) - $Gd_dR1;
$Z_gd = ( $Z_wgs-$Gd_dC1*$X_wgs+$Gd_dC0*$Y_wgs ) - $Gd_dR2;
/* Transform geo-centric coordinates to Lat/Long/H in Generic Datum */
$P = hypot($X_gd, $Y_gd);
$Q = sqrt(1.0 - $Bessel_Ellipsoid_e2p);
$R = $Bessel_Ellipsoid_a*$Bessel_Ellipsoid_e2p;
$Eta = atan( $Z_gd/($P*$Q) );
$lat_gd = atan( ($Z_gd + $R*pow(sin($Eta),3.0)/$Q) /
($P - $R*pow(cos($Eta),3.0) ) );
$long_gd = atan2( $Y_gd, $X_gd );
$N_gd = $Bessel_Ellipsoid_a /
sqrt( 1.0 - $Bessel_Ellipsoid_e2p*pow(sin($lat_gd), 2));
$h_gd = $P/cos($lat_gd) - $N_gd;
$a1 = 6.674372206E-3; /* Series expansion coefficients */
$a2 = 3.7073149E-5;
$a3 = 2.569374E-7;
$a4 = 1.9482E-9;
$b1 = 8.3522527E-4;
$b2 = 7.56302E-7;
$b3 = 1.193E-9;
$Y0 = 1500000.0; /* Mean Meridian Y Coordinate */
$RTH_L0 = 0.2759064963; /* Mean Meridian in RT90 (rad) */
$EarthRadius = 6366742.5194; /* [m] Mean Earth Radius */
$sin_lat = sin($lat_gd);
$sin2_lat = pow($sin_lat, 2);
$Lat_gc = $lat_gd -
$sin_lat*cos($lat_gd) *
($a1+$sin2_lat*($a2+$sin2_lat*($a3+$sin2_lat*$a4)));
$Xsi = atan(tan($Lat_gc) / cos($long_gd-$RTH_L0));
$Eta = atanh(cos($Lat_gc) * sin($long_gd-$RTH_L0));
$res["X"] = $EarthRadius*($Xsi+$b1*sin(2.0*$Xsi)*cosh(2.0*$Eta) +
$b2*sin(4.0*$Xsi)*cosh(4.0*$Eta) +
$b3*sin(6.0*$Xsi)*cosh(6.0*$Eta));
$res["Y"] = $Y0 + $EarthRadius*($Eta+$b1*cos(2.0*$Xsi)*sinh(2.0*$Eta) +
$b2*cos(4.0*$Xsi) * sinh(4.0*$Eta) +
$b3*cos(6.0*$Xsi) * sinh(6.0*$Eta));
return $res;
} /* RT90 */
?>
Det går säkert åt någonting med motsvarande komplexitetsnivå för att konvertera från Sweref99 till WGS84.
m v h Henrik