[转载]：经纬度与WGS84坐标转换

时间：2021-11-15 15:28:07

本代码实现在WGS84系统的大地坐标（BLH）和空间直角坐标（XYZ）的互相转换，符合标准语法，可直接使用

如下代码，输出为：
 WGS84:  -2175790.73969891    4461032.11207734     3992337.79032463
 BLH:   38.9999999999998    116.000000000000    33.0000069718808

Module CorrTrans

  !// WGS84 系统BLH坐标与空间直角坐标转换

  !// Fortran Coder http://fcode.cn

  !// Adapted from Acheng's C++ code

  Implicit None

  Integer , parameter :: DP = Selected_Real_Kind( 12 )

  Real(kind=DP) , parameter :: PI = 3.14159265358979323846_DP

  Real(kind=DP) , parameter , private :: a = 6378137._DP

  Real(kind=DP) , parameter , private :: f = 1.0_DP / 298.257223563_DP

  Real(kind=DP) , parameter , private :: t = a * ( 1.0_DP - f )

  Real(kind=DP) , parameter , private :: e = sqrt( (a*a - t*t) / (a*a) )

contains

  Subroutine XYZ2BLH( x , y , z , B , L , H )

    Real(Kind=DP) , Intent( IN )  :: x , y , z

    Real(Kind=DP) , Intent( OUT ) :: B , L , H

    real(kind=DP) :: N

    integer :: i

    L = atan( abs(y/x) )

    B = atan( abs(z/sqrt(x*x+y*y) ) )

    do i = 1 , 5

      N = a / sqrt( 1.0_DP - e*e*sin(B)*sin(B) )

      H = z / sin(B) - N * ( 1.0_DP - e*e )

      B = atan( z*(N+H) / ( sqrt(x*x+y*y)*(N*(1.0_DP-e*e)+H) ) )

    end do

    if ( x < 0._DP .and. y > 0._DP ) L = PI - L

    if ( x > 0._DP .and. y < 0._DP ) L = -1.0_DP * L

    if ( x < 0._DP .and. y < 0._DP ) L = -1.0 * PI + L

    B = B * 180._DP / PI

    L = L * 180._DP / PI

  End Subroutine XYZ2BLH

  Subroutine BLH2XYZ( B , L , H , x , y , z )

    Real(Kind=DP) , Intent( IN )  :: B , L , H

    Real(Kind=DP) , Intent( OUT ) :: x , y , z

    real(kind=DP) :: N , Br , Lr

    Br = B * PI / 180._DP

    Lr = L * PI / 180._DP

    N = a / sqrt( 1.0_DP - e*e*sin(Br)*sin(Br) )

    x = ( N + H ) * cos(Br)*cos(Lr)

    y = ( N + H ) * cos(Br)*sin(Lr)

    z = ( N*(1.0_DP - e*e ) + H ) * sin(Br);

  End Subroutine BLH2XYZ

End Module CorrTrans

Program www_fcode_cn

  Use CorrTrans

  Implicit None

  Real(Kind=DP)  :: x , y , z

  Real(Kind=DP)  :: B , L , H

  B = 39._DP ; L = 116._DP ; H = 33._DP

  call BLH2XYZ( B , L , H , x , y , z )

  write(*,*) 'WGS84:' , x , y , z

  call XYZ2BLH( x , y , z , B , L , H )

  write(*,*) 'BLH:' , B , L , H

End Program www_fcode_cn

转自：http://fcode.cn/code_prof-89-1.html

相关文章