wgs84 gc02 bd09 坐标系转换excel宏

Private Const PI As Double = 3.14159265358979
Private Const x_pi = 3.14159265358979 * 3000# / 180#
Private Const X_EARTH_RADIUS = 6378245
Private Const EE = 6.69342162296594E-03
 
Private Type Point
  lng     As Double
  lat     As Double
End Type
Private Function transform(lng As Double, lat As Double) As Point
   dlat = transformLat(lng - 105#, lat - 35#)
   dlng = transformlng(lng - 105#, lat - 35#)
   radlat = lat / 180# * PI
   magic = Math.Sin(radlat)
   magic = 1 - EE * magic * magic
   sqrtmagic = Sqr(magic)
   dlat = (dlat * 180#) / ((X_EARTH_RADIUS * (1 - EE)) / (magic * sqrtmagic) * PI)
   dlng = (dlng * 180#) / (X_EARTH_RADIUS / sqrtmagic * Math.Cos(radlat) * PI)
   mgLat = lat + dlat
   mglng = lng + dlng
   Dim ret As Point
   ret.lng = mglng
   ret.lat = mgLat
   transform = ret
 End Function
Private Function transformlng(x As Double, y As Double) As Double
   ret = 300# + x + 2# * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Sqr(Math.Abs(x))
   ret = ret + (20# * Math.Sin(6# * x * PI) + 20# * Math.Sin(2# * x * PI)) * 2# / 3#
   ret = ret + (20# * Math.Sin(x * PI) + 40# * Math.Sin(x / 3# * PI)) * 2# / 3#
   ret = ret + (150# * Math.Sin(x / 12# * PI) + 300# * Math.Sin(x / 30# * PI)) * 2# / 3#
   transformlng = ret
 End Function
Private Function transformLat(x As Double, y As Double) As Double
   ret = -100# + 2# * x + 3# * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Sqr(Math.Abs(x))
   ret = ret + (20# * Math.Sin(6# * x * PI) + 20# * Math.Sin(2# * x * PI)) * 2# / 3#
   ret = ret + (20# * Math.Sin(y * PI) + 40# * Math.Sin(y / 3# * PI)) * 2# / 3#
   ret = ret + (160# * Math.Sin(y / 12# * PI) + 320 * Math.Sin(y * PI / 30#)) * 2# / 3#
   transformLat = ret
 End Function
Private Function BD09toGCJ02(lng As Double, lat As Double) As Point
    Dim p As Point
    x = lng - 0.0065
    y = lat - 0.006
    Z = Sqr(x * x + y * y) - 0.00002 * Math.Sin(y * x_pi)
    theta = WorksheetFunction.Atan2(x, y) - 0.000003 * Math.Cos(x * x_pi)
    p.lng = Z * Math.Cos(theta)
    p.lat = Z * Math.Sin(theta)
    BD09toGCJ02 = p
End Function
Private Function GCJ02toBD09(lng As Double, lat As Double) As Point
    Dim p As Point
    Z = Sqr(lng * lng + lat * lat) + 0.00002 * Math.Sin(lat * x_pi)
    theta = WorksheetFunction.Atan2(lng, lat) + 0.000003 * Math.Cos(lng * x_pi)
    p.lng = Z * Math.Cos(theta) + 0.0065
    p.lat = Z * Math.Sin(theta) + 0.006
    GCJ02toBD09 = p
End Function
 
Private Function WGS84toGCJ02(lng As Double, lat As Double) As Point
    Dim p As Point
    dlat = transformLat(lng - 105#, lat - 35#)
    dlng = transformlng(lng - 105#, lat - 35#)
    radlat = lat / 180# * PI
    magic = Math.Sin(radlat)
    magic = 1 - EE * magic * magic
    sqrtmagic = Sqr(magic)
    dlat = (dlat * 180#) / ((X_EARTH_RADIUS * (1 - EE)) / (magic * sqrtmagic) * PI)
    dlng = (dlng * 180#) / (X_EARTH_RADIUS / sqrtmagic * Math.Cos(radlat) * PI)
    p.lat = lat + dlat
    p.lng = lng + dlng
    WGS84toGCJ02 = p
End Function
 
Public Function BD09toGCJ02Lng(lng As Double, lat As Double) As Double
    Dim p As Point
    p = BD09toGCJ02(lng, lat)
    BD09toGCJ02Lng = p.lng
End Function
Public Function BD09toGCJ02Lat(lng As Double, lat As Double) As Double
    Dim p As Point
    p = BD09toGCJ02(lng, lat)
    BD09toGCJ02Lat = p.lat
End Function
    
Public Function GCJ02toBD09Lat(lng As Double, lat As Double) As Double
    Dim p As Point
    p = GCJ02toBD09(lng, lat)
    GCJ02toBD09Lat = p.lat
End Function
    
Public Function GCJ02toBD09Lng(lng As Double, lat As Double) As Double
    Dim p As Point
    p = GCJ02toBD09(lng, lat)
    GCJ02toBD09Lng = p.lng
End Function
    
    
Public Function WGS84toGCJ02Lng(lng As Double, lat As Double) As Double
    Dim p As Point
    p = WGS84toGCJ02(lng, lat)
    WGS84toGCJ02Lng = p.lng
End Function
     
      
Public Function WGS84toGCJ02Lat(lng As Double, lat As Double) As Double
    Dim p As Point
    p = WGS84toGCJ02(lng, lat)
    WGS84toGCJ02Lat = p.lat
End Function
     
 
Public Function GCJ02toWGS84Lng(lng As Double, lat As Double) As Double
    Dim p As Point
    p = transform(lng, lat)
    mglng = lng * 2 - p.lng
    mgLat = lat * 2 - p.lat
    GCJ02toWGS84Lng = mglng
End Function
Public Function GCJ02toWGS84Lat(lng As Double, lat As Double) As Double
    Dim p As Point
    p = transform(lng, lat)
    mglng = lng * 2 - p.lng
    mgLat = lat * 2 - p.lat
    GCJ02toWGS84Lat = mgLat
End Function
 
Public Function WGS84toBD09Lat(lng As Double, lat As Double) As Double
    Dim p As Point
    p = WGS84toGCJ02(lng, lat)
    WGS84toBD09Lat = GCJ02toBD09Lat(p.lng, p.lat)
End Function
      
Public Function WGS84toBD09Lng(lng As Double, lat As Double) As Double
   Dim p As Point
    p = WGS84toGCJ02(lng, lat)
    WGS84toBD09Lng = GCJ02toBD09Lng(p.lng, p.lat)
End Function
 
Public Function BD09toWGS84Lat(lng As Double, lat As Double) As Double
    Dim p As Point
    p = BD09toGCJ02(lng, lat)
    BD09toWGS84Lat = GCJ02toWGS84Lat(p.lng, p.lat)
End Function
      
Public Function BD09toWGS84Lng(lng As Double, lat As Double) As Double
   Dim p As Point
    p = BD09toGCJ02(lng, lat)
    BD09toWGS84Lng = GCJ02toWGS84Lng(p.lng, p.lat)
End Function