728x90
OPINET의 API를 이용하여 유가 정보를 조회하는 기능을 추가하던중
주유소의 좌표와 인근 주유소를 조회하는 API의 좌표가 KATEC이라고 하여
근 일주일 정도를 좌표계가 뭔지, 상호 변환은 어떻게 하는지 머리를 싸맸네요...
네이버나 카카오, SK 등에서 좌표 변환 API를 지원하기는 하는데,
사용자별로 Key를 신청해서 써야 하는지라 범용으로 사용하기 어려운 점이 있습니다.
여기저기 검색하다가 아래 페이지를 검색하여 좌표 변환 소스를 얻었습니다.
https://www.androidpub.com/android_dev_info/1318647
https://javaexpert.tistory.com/142
하지만 VBA로 포팅을 해 봤는데, JAVA 관련 지식이 미천하여 시행착오를 일주일,,,
제가 사용할 수 있는 좌표변환 API를 사용하여 전국 좌표를 많이 샘플링하여 좌표를 변환하고,
계산 결과의 차이 평균이 최소가 되도록 일부 상수를 조금 수정했습니다.
코드를 조금 더 효율적으로 수정하고, 변환 가능한 CRS를 몇 개 추가하였습니다.
혹시 잘아시는 분이 계시면 제가 잘 못 수정했거나 틀린 정보가 있으면 알려 주세요.
원본에서는 GeoTransPoint 클래스를 사용하는데, X,Y,Z 좌표값만 사용하는 것이라 변수로 변경했습니다.
사용 가능한 좌표계는 6개입니다.
아래 설명부분은 제가 찾아 보면서 정리한 것이라 공식적인 이름과 다를 수도 있으며 적용사례도 틀릴 수 있습니다.
WGS84 타원체 사용 좌표계
- 0 : EPSG:4326, WGS84위경도, GPS 위경도 좌표계
- 1 : EPSG:32652, UTM52N
GRS80 타원체 사용 좌표계
- 2 : EPSG:5179, UTM-K, NHN:2048, 네이버 지도
- 3 : EPSG:5181, 타원체를 바꾼 중부원점, 다음지도
Bessel 타원체 사용 좌표계
- 4 : EPSG:5174, 2002년 이전의 중부원점, 2002년 이전 연속지적도
- 5 : EPSG없음, KATEC, 네비게이션용
실제 사용 전에 기준 변환 자료와 대조하여 오류나 오차값을 확인하시기 바랍니다.
먼저 사용할 변수들입니다. 이해할 수 있는 것만 코메트를 달았습니다.
Option Explicit
Option Base 0
'//
'// 원저자(행복한 수지아빠) : https://javaexpert.tistory.com/142
'// 보완자(김한국) : https://www.androidpub.com/android_dev_info/1318647
'//
'// EGTools : https://egtools.tistory.com/entry/VBAGEOConverter
'// egexcelvba@gmail.com 2024-02-07, VBA로포팅, KATEC 점검
'// 2024-02-10, 기본 좌표계 상수 및 SetVariables, 좌표계 정리 및 추가
'//
'
' * @author aquilegia
' *
' * The code based on hyosang(http://hyosang.kr/tc/96) and aero's blog ((http://aero.sarang.net/map/analysis.html)
' * License: LGPL : http://www.gnu.org/copyleft/lesser.html
' */
'// WGS84 타원체 사용 좌표계
Private Const WGS84 = 0 '// EPSG:4326, WGS84, GPS 위경도 좌표계
Private Const UTM52N = 1 '// EPSG:32652, UTM52N
'// GRS80 타원체 사용 좌표계
Private Const UTM_K = 2 '// EPSG:5179, UTM-K, NHN:2048, 네이버 지도
Private Const GRS80 = 3 '// EPSG:5181, 타원체를 바꾼 중부원점, 다음지도
'// Bessel 타원체 사용 좌표계
Private Const B5174 = 4 '// EPSG:5174, 다음TM, 2002년 이전의 중부원점, 2002년 이전 연속지적도
Private Const KATEC = 5 '// EPSG없음, KATEC, 네비게이션용
Private Const gX = 0
Private Const gY = 1
Private Const gZ = 2
'// 지구타원체
Private Const WGS84a = 6378137# '// WGS84타원체
Private Const WGS84b = 6356752.31424518 '// 298.257223563
Private Const GRS80a = 6378137# '// GRS80타원체
Private Const GRS80b = 6356752.31414035 '// 298.257222101, GRS80은 WGS84와 아주 조금 차이남
Private Const Bessel1841a = 6377397.155 '// Bessel1841 타원체
Private Const Bessel1841b = 6356078.96282175 '// 299.15281285
Private m_ind(5) As Double
Private m_Es(5) As Double
Private m_Esp(5) As Double
Private src_m(5) As Double '// 경도원점의 적도 이격 거리(m)
Private dst_m(5) As Double '// 경도원점의 적도 이격 거리(m)
Private Const EPSLN = 0.0000000001
Private m_arMajor(5) As Double '// 타원체 최대 지름
Private m_arMinor(5) As Double '// 타원체 최소 지름
Private m_arScaleFactor(5) As Double '//
Private m_arLonCenter(5) As Double '// 기준점 경도 기준의 Radian
Private m_arLatCenter(5) As Double '// 기준점 위도의 Radian
Private m_arFalseEasting(5) As Double '// 경도 좌표 원점 이격
Private m_arFalseNorthing(5) As Double '// 위도 좌표 원점 이격
Private datum_params(2) As Double
Private Const PI = 3.14159265358979 '// 3.14159265358979323846...
Private Const HALF_PI = 1.5707963267949 '// 1.57079632679489661923...
Private Const COS_67P5 = 0.38268343236509 '/* cosine of 67.5 degrees */
Private Const AD_C = 1.0026
공통으로 사용할 수치변환 함수들입니다.
Private Function D2R(ByVal degree As Double)
D2R = degree * PI / 180
End Function
Private Function R2D(ByVal radian As Double)
R2D = radian * 180 / PI
End Function
Private Function e0fn(ByVal x As Double)
e0fn = 1 - 0.25 * x * (1 + x / 16 * (3 + 1.25 * x))
End Function
Private Function e1fn(ByVal x As Double)
e1fn = 0.375 * x * (1 + 0.25 * x * (1 + 0.46875 * x))
End Function
Private Function e2fn(ByVal x As Double)
e2fn = 0.05859375 * x * x * (1 + 0.75 * x)
End Function
Private Function e3fn(ByVal x As Double)
e3fn = x * x * x * (35 / 3072)
End Function
Private Function mlfn(ByVal e0 As Double, ByVal e1 As Double, ByVal e2 As Double, ByVal e3 As Double, ByVal phi As Double)
mlfn = e0 * phi - e1 * Sin(2 * phi) + e2 * Sin(4 * phi) - e3 * Sin(6 * phi)
End Function
Private Function asinz(ByVal value As Double)
If (Abs(value) > 1) Then value = IIf(value > 0, 1, -1)
asinz = Application.Asin(value)
End Function
많은 변수들의 초기값을 설정합니다.
Private Sub SetVariables()
m_arScaleFactor(WGS84) = 1# '// EPSG:4326, WGS84, GPS 위경도 좌표계
m_arLonCenter(WGS84) = 0
m_arLatCenter(WGS84) = 0
m_arFalseNorthing(WGS84) = 0
m_arFalseEasting(WGS84) = 0
m_arMajor(WGS84) = WGS84a
m_arMinor(WGS84) = WGS84b
m_arScaleFactor(UTM52N) = 0.9996 '// EPSG:32652, UTM52N
m_arLonCenter(UTM52N) = 2.25147473507268 '// 129
m_arLatCenter(UTM52N) = 0 '// 적도
m_arFalseEasting(UTM52N) = 500000#
m_arFalseNorthing(UTM52N) = 0
m_arMajor(UTM52N) = WGS84a
m_arMinor(UTM52N) = WGS84b
m_arScaleFactor(UTM_K) = 0.9996 '// EPSG:5179, UTM-K, NHN:2048, 네이버 지도
m_arLonCenter(UTM_K) = 2.22529479629277 '// 127.5
m_arLatCenter(UTM_K) = 0.663225115757845 '// 38
m_arFalseEasting(UTM_K) = 1000000# '//
m_arFalseNorthing(UTM_K) = 2000000# '//
m_arMajor(UTM_K) = GRS80a
m_arMinor(UTM_K) = GRS80b
m_arScaleFactor(GRS80) = 1# '// EPSG:5181, 타원체를 바꾼 중부원점, 다음 지도
m_arLonCenter(GRS80) = 2.2165681500328 '// 127
m_arLatCenter(GRS80) = 0.663225115757845 '// 38
m_arFalseEasting(GRS80) = 200000# '//
m_arFalseNorthing(GRS80) = 500000# '//
m_arMajor(GRS80) = GRS80a
m_arMinor(GRS80) = GRS80b
m_arScaleFactor(B5174) = 1# '// EPSG:5174, 다음TM, 2002년 이전의 중부원점, 2002년 이전 연속지적도 -> Tmap에 최적화
m_arLonCenter(B5174) = 2.21661859489632 '// 2.21661859489671(127+10.485 minute)
m_arLatCenter(B5174) = 0.663225115757845 '// 38
m_arFalseEasting(B5174) = 200000# '//
m_arFalseNorthing(B5174) = 500000# '//
m_arMajor(B5174) = Bessel1841a
m_arMinor(B5174) = Bessel1841b
m_arScaleFactor(KATEC) = 0.9999 '// EPSG없음, KATEC, 네비게이션용 -> Tmap에 최적화
m_arLonCenter(KATEC) = 2.23402144255274 '// 128
m_arLatCenter(KATEC) = 0.663225115757845 '// 38
m_arFalseEasting(KATEC) = 400000# '//
m_arFalseNorthing(KATEC) = 600000# '//
m_arMajor(KATEC) = Bessel1841a
m_arMinor(KATEC) = Bessel1841b
'// Datum Transformation Parameters Local Geodetic(Bessel 1841) to WGS84
'// 원본 : -146.43, 507.89, 681.46
'// 여러 웹 소개 자료에 나오는 Bursa-Wolf 모델의 Proj4용 파라미터
'// +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43
'// https://epsg.io/5178에 수록되 정보
'// TOWGS84[-145.907,505.034,685.756,-1.162,2.347,1.592,6.342]], Revision date: 2011-01-25
'// 3개 API에 대조한 가장 근사치는 epsg.io/5178에 수록된 정보임
datum_params(gX) = -145.907
datum_params(gY) = 505.034
datum_params(gZ) = 685.756
Dim i&, tmp#
For i = LBound(m_Es) To UBound(m_Es)
tmp = m_arMinor(i) / m_arMajor(i)
m_Es(i) = 1# - tmp * tmp
m_Esp(i) = m_Es(i) / (1# - m_Es(i))
If (m_Es(i) < 0.00001) Then m_ind(i) = 1# Else m_ind(i) = 0
src_m(i) = m_arMajor(i) * mlfn(e0fn(m_Es(i)), e1fn(m_Es(i)), e2fn(m_Es(i)), e3fn(m_Es(i)), m_arLatCenter(i))
dst_m(i) = m_arMajor(i) * mlfn(e0fn(m_Es(i)), e1fn(m_Es(i)), e2fn(m_Es(i)), e3fn(m_Es(i)), m_arLatCenter(i))
Next i
End Sub
좌표 변환을 수행하는 함수들입니다.
Main 함수는 convert() 입니다.
원본의 코멘트 및 원본이 참조하는 자료의 코멘트는 가능한 그대로 두었습니다.
Private Function convert(ByVal srctype As Integer, ByVal dsttype As Integer, ByVal in_pt As Variant)
Dim tmpPt, out_pt
If srctype = WGS84 Then '// WGS84, GRS80 타원체, 기준원점 그대로
ReDim tmpPt(2) As Double
tmpPt(gX) = D2R(in_pt(gX))
tmpPt(gY) = D2R(in_pt(gY))
Else
tmpPt = tm2geo(srctype, in_pt)
End If
If dsttype = WGS84 Then '// WGS84, GRS80 타원체, 기준원점 그대로
ReDim out_pt(2) As Double
out_pt(gX) = R2D(tmpPt(gX))
out_pt(gY) = R2D(tmpPt(gY))
Else
out_pt = geo2tm(dsttype, tmpPt)
End If
convert = out_pt
End Function
Private Function geo2tm(ByVal dsttype As Integer, ByVal in_pt As Variant)
Dim x#, y#, delta_lon#, sin_phi#, cos_phi#, b#, con#, al#, als#, c#, t#, tq#, n#, ml#, out_pt#(2), tmpx#
in_pt = transform(WGS84, dsttype, in_pt)
delta_lon = in_pt(gX) - m_arLonCenter(dsttype)
sin_phi = Sin(in_pt(gY))
cos_phi = Cos(in_pt(gY))
If (m_ind(dsttype) <> 0) Then
b = cos_phi * Sin(delta_lon)
If ((Abs(Math.Abs(b) - 1#)) < EPSLN) Then
Debug.Print "무한대 에러"
End If
Else
b = 0
x = 0.5 * m_arMajor(dsttype) * m_arScaleFactor(dsttype) * Log((1 + b) / (1 - b))
con = Application.Acos(cos_phi * Cos(delta_lon) / Sqr(1 - b * b))
If (in_pt(gY) < 0) Then
con = con * -1
y = m_arMajor(dsttype) * m_arScaleFactor(dsttype) * (con - m_arLatCenter(dsttype))
End If
End If
al = cos_phi * delta_lon
als = al * al
c = m_Esp(dsttype) * cos_phi * cos_phi
tq = Tan(in_pt(gY))
t = tq * tq
con = 1# - m_Es(dsttype) * sin_phi * sin_phi
n = m_arMajor(dsttype) / Sqr(con)
ml = m_arMajor(dsttype) * mlfn(e0fn(m_Es(dsttype)), e1fn(m_Es(dsttype)), e2fn(m_Es(dsttype)), e3fn(m_Es(dsttype)), in_pt(gY))
out_pt(gX) = m_arScaleFactor(dsttype) * n * al * (1 + als / 6 * (1 - t + c + als / 20 * (5 - 18 * t + t * t + 72 * c - 58 * m_Esp(dsttype)))) + m_arFalseEasting(dsttype)
tmpx = (0.5 + als / 24 * (5 - t + 9 * c + 4 * c * c + als / 30 * (61 - 58 * t + t * t + 600 * c - 330 * m_Esp(dsttype))))
out_pt(gY) = m_arScaleFactor(dsttype) * (ml - dst_m(dsttype) + n * tq * (als * tmpx)) + m_arFalseNorthing(dsttype)
geo2tm = out_pt
End Function
Private Function tm2geo(ByVal srctype As Integer, ByVal in_pt As Variant)
Dim max_iter%, f#, g#, temp#, h#, con#, phi#, i%, delta_Phi#, tmpPt#(2), out_pt, tmpx#
Dim sin_phi#, cos_phi#, tan_phi#, c#, cs#, t#, ts#, cont#, n#, R#, d#, ds#
ReDim out_pt(2)
tmpPt(gX) = in_pt(gX)
tmpPt(gY) = in_pt(gY)
max_iter = 6
If (m_ind(srctype) <> 0) Then
f = Exp(in_pt(gX) / (m_arMajor(srctype) * m_arScaleFactor(srctype)))
g = 0.5 * (f - 1 / f)
temp = m_arLatCenter(srctype) + tmpPt(gY) / (m_arMajor(srctype) * m_arScaleFactor(srctype))
h = Cos(temp)
con = Sqr((1 - h * h) / (1 + g * g))
out_pt(gY) = asinz(con)
If (temp < 0) Then out_pt(gY) = out_pt(gY) * (-1)
If ((g = 0) & (h = 0)) Then
out_pt(gX) = m_arLonCenter(srctype)
Else
out_pt(gX) = Application.atan(g / h) + m_arLonCenter(srctype)
End If
End If
tmpPt(gX) = tmpPt(gX) - m_arFalseEasting(srctype)
tmpPt(gY) = tmpPt(gY) - m_arFalseNorthing(srctype)
con = (src_m(srctype) + tmpPt(gY) / m_arScaleFactor(srctype)) / m_arMajor(srctype)
phi = con
i = 0
Do While (True)
delta_Phi = ((con + e1fn(m_Es(srctype)) * Sin(2 * phi) - e2fn(m_Es(srctype)) * Sin(4 * phi) + e3fn(m_Es(srctype)) * Sin(6 * phi)) / e0fn(m_Es(srctype))) - phi
phi = phi + delta_Phi
If (Math.Abs(delta_Phi) <= EPSLN) Then Exit Do
If (i >= max_iter) Then
Debug.Print "무한대 에러"
Exit Function
End If
i = i + 1
Loop
If Abs(phi) < (PI / 2) Then
sin_phi = Sin(phi)
cos_phi = Cos(phi)
tan_phi = Tan(phi)
c = m_Esp(srctype) * cos_phi * cos_phi
cs = c * c
t = tan_phi * tan_phi
ts = t * t
cont = 1 - m_Es(srctype) * sin_phi * sin_phi
n = m_arMajor(srctype) / Sqr(cont)
R = n * (1 - m_Es(srctype)) / cont
d = tmpPt(gX) / (n * m_arScaleFactor(srctype))
ds = d * d
tmpx = (5 + 3 * t + 10 * c - 4 * cs - 9 * m_Esp(srctype) - ds / 30 * (61 + 90 * t + 298 * c + 45 * ts - 252 * m_Esp(srctype) - 3 * cs))
out_pt(gY) = phi - (n * tan_phi * ds / R) * (0.5 - ds / 24 * tmpx)
tmpx = (5 - 2 * c + 28 * t - 3 * cs + 8 * m_Esp(srctype) + 24 * ts)
out_pt(gX) = m_arLonCenter(srctype) + (d * (1 - ds / 6 * (1 + 2 * t + c - ds / 20 * tmpx)) / cos_phi)
Else
out_pt(gY) = PI * 0.5 * Sin(tmpPt(gY))
out_pt(gX) = m_arLonCenter(srctype)
End If
out_pt = transform(srctype, WGS84, out_pt)
tm2geo = out_pt
End Function
사용하는 타원체에 따라 상호 변환을 해 주는 함수들입니다.
'//
'// Author: Richard Greenwood rich@greenwoodmap.com
'// License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
'//
'// * convert between geodetic coordinates (longitude, latitude, height)
'// * and gecentric coordinates (X, Y, Z)
'// * ported from Proj 4.9.9 geocent.c
'//
'// following constants from geocent.c
'// private static final double HALF_PI = 0.5 * Math.PI;
'// private static final double COS_67P5 = 0.38268343236508977; /* cosine of 67.5 degrees */
'// private static final double AD_C = 1.0026000 ;
'// /* Toms region 1 constant */
'//
Private Function transform(ByVal srctype As Integer, ByVal dsttype As Integer, ByVal point As Variant)
If srctype = dsttype Then Exit Function
If srctype > GRS80 Or dsttype > GRS80 Then
'// Convert to geocentric coordinates.
point = geodetic_to_geocentric(srctype, point)
'// Convert between datums
If srctype > GRS80 Then point = geocentric_to_wgs84(point)
If dsttype > GRS80 Then point = geocentric_from_wgs84(point)
'// Convert back to geodetic coordinates
point = geocentric_to_geodetic(dsttype, point)
End If
transform = point
End Function
Private Function geodetic_to_geocentric(ByVal iType As Integer, ByVal p As Variant)
'// The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
'// (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
'// according to the current ellipsoid parameters.
'//
'// Latitude : Geodetic latitude in radians (input)
'// Longitude : Geodetic longitude in radians (input)
'// Height : Geodetic height, in meters (input)
'// X : Calculated Geocentric X coordinate, in meters (output)
'// Y : Calculated Geocentric Y coordinate, in meters (output)
'// Z : Calculated Geocentric Z coordinate, in meters (output)
Dim Longitude#, Latitude#, Height#, x#, y#, Z#, Rn#, Sin_Lat#, Sin2_Lat#, Cos_Lat#
Longitude = p(gX)
Latitude = p(gY)
Height = p(gZ)
'Rn '// Earth radius at location
'Sin_Lat '// Sin(Latitude)
'Sin2_Lat '// Square of Sin(Latitude)
'Cos_Lat '// Cos(Latitude)
'// Don't blow up if Latitude is just a little out of the value
'// range as it may just be a rounding issue.
'// Also removed longitude test, it should be wrapped by Cos() and Sin().
'// NFW for PROJ.4, Sep/2001.
If (Latitude < -HALF_PI And Latitude > -1.001 * HALF_PI) Then
Latitude = -HALF_PI
ElseIf (Latitude > HALF_PI And Latitude < 1.001 * HALF_PI) Then
Latitude = HALF_PI
ElseIf ((Latitude < -HALF_PI) Or (Latitude > HALF_PI)) Then '/* Latitude out of range */
geodetic_to_geocentric = True
End If
'/* no errors */
If (Longitude > PI) Then Longitude = Longitude - (2 * PI)
Sin_Lat = Sin(Latitude)
Cos_Lat = Cos(Latitude)
Sin2_Lat = Sin_Lat * Sin_Lat
Rn = m_arMajor(iType) / (Sqr(1 - m_Es(iType) * Sin2_Lat))
x = (Rn + Height) * Cos_Lat * Cos(Longitude)
y = (Rn + Height) * Cos_Lat * Sin(Longitude)
Z = ((Rn * (1 - m_Es(iType))) + Height) * Sin_Lat
p(gX) = x
p(gY) = y
p(gZ) = Z
geodetic_to_geocentric = p
End Function
Private Function geocentric_to_geodetic(ByVal iType As Integer, ByVal p As Variant)
'// Convert_Geocentric_To_Geodetic
'// The method used here is derived from 'An Improved Algorithm for
'// Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
Dim x#, y#, Z#, Longitude#, Latitude#, Height#, W#, W2#, T0#, T1#, S0#, S1#, Sin_B0#, Sin3_B0#, Cos_B0#, Sin_P1#, Cos_P1#, Rn#, Sum#, At_Pole As Boolean
x = p(gX)
y = p(gY)
Z = p(gZ)
Longitude = 0
Latitude = 0
Height = 0
'W '* distance from Z axis */
'W2 '* square of distance from Z axis */
'T0 '* initial estimate of vertical component */
'T1 '* corrected estimate of vertical component */
'S0 '* initial estimate of horizontal component */
'S1 '* corrected estimate of horizontal component */
'Sin_B0 '* Math.sin(B0), B0 is estimate of Bowring aux doubleiable */
'Sin3_B0 '* cube of Math.sin(B0) */
'Cos_B0 '* Math.cos(B0) */
'Sin_P1 '* Math.sin(phi1), phi1 is estimated latitude */
'Cos_P1 '* Math.cos(phi1) */
'Rn '* Earth radius at location */
'Sum '* numerator of Math.cos(phi1) */
'At_Pole '* indicates location is in polar region */
At_Pole = False
If (x <> 0) Then
Longitude = Application.Atan2(x, y) ' Application.Atan2(y, x)
Else
If (y > 0) Then
Longitude = HALF_PI
ElseIf (y < 0) Then
Longitude = -HALF_PI
Else
At_Pole = True
Longitude = 0
If (Z > 0) Then '/* north pole */
Latitude = HALF_PI
ElseIf (Z < 0) Then '/* south pole */
Latitude = -HALF_PI
Else '/* center of earth */
Latitude = HALF_PI
Height = -m_arMinor(iType)
Exit Function
End If
End If
End If
W2 = x * x + y * y
W = Sqr(W2)
T0 = Z * AD_C
S0 = Sqr(T0 * T0 + W2)
Sin_B0 = T0 / S0
Cos_B0 = W / S0
Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0
T1 = Z + m_arMinor(iType) * m_Esp(iType) * Sin3_B0
Sum = W - m_arMajor(iType) * m_Es(iType) * Cos_B0 * Cos_B0 * Cos_B0
S1 = Sqr(T1 * T1 + Sum * Sum)
Sin_P1 = T1 / S1
Cos_P1 = Sum / S1
Rn = m_arMajor(iType) / Sqr(1 - m_Es(iType) * Sin_P1 * Sin_P1)
If (Cos_P1 >= COS_67P5) Then
Height = W / Cos_P1 - Rn
ElseIf (Cos_P1 <= -COS_67P5) Then
Height = W / -Cos_P1 - Rn
Else
Height = Z / Sin_P1 + Rn * (m_Es(iType) - 1)
End If
If (At_Pole = False) Then Latitude = Atn(Sin_P1 / Cos_P1)
p(gX) = Longitude
p(gY) = Latitude
p(gZ) = Height
geocentric_to_geodetic = p
End Function
Private Function geocentric_to_wgs84(ByVal p As Variant)
'// geocentic_to_wgs84(defn, p )
'// defn = coordinate system definition,
'// p = point to transform in geocentric coordinates (x,y,z)
p(gX) = p(gX) + datum_params(gX)
p(gY) = p(gY) + datum_params(gY)
p(gZ) = p(gZ) + datum_params(gZ)
geocentric_to_wgs84 = p
End Function
Private Function geocentric_from_wgs84(ByVal p As Variant)
'// geocentic_from_wgs84()
'// coordinate system definition,
'// point to transform in geocentric coordinates (x,y,z)
p(gX) = p(gX) - datum_params(gX)
p(gY) = p(gY) - datum_params(gY)
p(gZ) = p(gZ) - datum_params(gZ)
geocentric_from_wgs84 = p
End Function
실제 Excel에서 호출하여 사용할 수 있는 함수를 별도로 추가하였습니다.
Function GEOConvert(ByVal Src As Integer, ByVal Dst As Integer, ByVal x As Double, ByVal y As Double)
Dim in_pt#(2), out_pt
If Src < 0 Or Src > 5 Or Dst < 0 Or Dst > 5 Then GEOConvert = CVErr(xlErrValue): Exit Function
in_pt(gX) = x: in_pt(gY) = y
SetVariables
out_pt = convert(Src, Dst, in_pt)
'// 사용하는 버전에 따라 아래 출력 방식 선택
GEOConvert = Array(out_pt(gX), out_pt(gY)) '//Excel2021 이상에서는 2개 셀에 동시 출력
'GEOConvert = out_pt(gX) & "," & out_pt(gY) '//Excel2019 이하 - 값을 쉼표로 구분하여 텍스트로 출력
End Function
728x90
'Excel - VBA' 카테고리의 다른 글
| Google Calendar를 Excel에 연결하기 (2) | 2024.03.18 |
|---|---|
| vWorld.kr의 Open API를 이용한 Geocoding (0) | 2024.02.12 |
| [Excel] Windows 11 에서 한글 입력 이상할 때 (1) | 2023.12.26 |
| [VBA] Module이나 Sub의 크기 한계 (3) | 2023.12.10 |
| [Excel] 기존 리본 메뉴에 Custom 그룹 삽입하기 (0) | 2023.11.27 |