Svar till: IGO 8

Start Forum Diverse forumfrågor IGO 8 Svar till: IGO 8

#220
henca
Moderator

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