zoukankan      html  css  js  c++  java
  • [转载]:经纬度与WGS84坐标转换

    本代码实现在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

  • 相关阅读:
    webpack
    react 原理
    jest
    input 自动获取焦点
    taro
    html5标签
    webpack
    每日日报
    每日日报
    每日日报
  • 原文地址:https://www.cnblogs.com/China3S/p/4692690.html
Copyright © 2011-2022 走看看