Go: WGS84를 WCONGNAMUL로 변환 함수

Go언어 백엔드에서 카카오맵 API를 이용하다 보니, WCONGNAMUL로 전환하는 일도 꽤 생겼다.

근데 외부 API를 부르자니 뭔가 싫고 공부겸 값을 찾아보기로 했다.

(사실 마음 편하게 카카오 API를 이용하자: https://developers.kakao.com/docs/latest/ko/local/dev-guide#trans-coord)

 

우선 변환된 값들을 이용해서 python numpy의 lstsq로 linear 하게 계수를 찾아보는 방법을 택했다.

numpy.linalg.lstsq — NumPy v1.26 Manual

from numpy import array
from numpy.linalg import lstsq

# WGS84 좌표 (longitude, latitude)
wgs84_coords = array([
    [126.99116337285824, 37.248098895147216, 1],
    [127.28434974725708, 36.762436923118536, 1],
    [127.37264182214031, 35.73294563400083, 1],
    [127.37264182214031, 35.7328, 1]
])

# WCONGNAMUL 좌표 (x, y)
wcongnamul_x = array([498040.0, 563473.0, 584279.0, 584280.0])
wcongnamul_y = array([1041367.0, 906718.0, 621193.0, 621153.0])

# x coefficients (A, B, C)
x_coeffs = lstsq(wgs84_coords, wcongnamul_x, rcond=None)[0]

# y coefficients (D, E, F)
y_coeffs = lstsq(wgs84_coords, wcongnamul_y, rcond=None)[0]

print(x_coeffs, y_coeffs)

"""
아래처럼 나옴
(array([ 2.21113038e+05, -1.24710175e+03, -2.75349098e+07]),
 array([ 1.87549269e+02,  2.77361611e+05, -9.31364281e+06]))
"""

 

계수 테스트 (Go언어)

func ConvertWGS84ToWCONGNAMUL(lat, long float64) WCONGNAMULCoord {
	x := 226244.957*long + 32.0744737*lat - 28234282.6
	y := 242.017919*long + 277367.474*lat - 9320788.94
	return WCONGNAMULCoord{X: x, Y: y}
}

func TestWCONGNAMUL(t *testing.T) {
	tests := []struct {
		name                 string
		lat, long            float64 // Starting point
		expectedX, expectedY float64 // Expected ending point
	}{
		{
			name:      "Test 1",
			lat:       37.248098895147216,
			long:      126.99116337285824,
			expectedX: 498040.0,
			expectedY: 1041367.0,
		},
		{
			name:      "Test Not korea 1",
			lat:       33.248098895147216,
			long:      126.99116337285824,
			expectedX: 497941.0,
			expectedY: -68085.0,
		},
		{
			name:      "Test 3",
			lat:       35.73294563400083,
			long:      127.37264182214031,
			expectedX: 584279.0,
			expectedY: 621193.0,
		},
		{
			name:      "Test 4",
			lat:       35.7328,
			long:      127.37264182214031,
			expectedX: 584280.0,
			expectedY: 621153.0,
		},
		{
			name:      "Test 5",
			lat:       34.248098895147216,
			long:      126.96666,
			expectedX: 492322.0,
			expectedY: 209211.0,
		},
	}

	for _, tt := range tests {
		t.Run(tt.name, func(t *testing.T) {
			result := ConvertWGS84ToWCONGNAMUL(tt.lat, tt.long)

			t.Logf("Calculated WCONGNAMUL: %f, %f", result.X, result.Y)
			if math.Abs(result.X-tt.expectedX) > 1 || math.Abs(result.Y-tt.expectedY) > 1 { // Allowing a margin of error of 1
				t.Errorf("ConvertWGS84ToWCONGNAMUL(%v, %v) = (%v, %v), want (%v, %v)",
					tt.lat, tt.long, result.X, result.Y, tt.expectedX, tt.expectedY)
			}
		})
	}
}

 

카카오 좌표 API를 이용해 1000개 정도의 샘플을 만들고 값을 돌려도 점점 더 차이가 나서 쓰기 힘들었다.

이런저런 소스를 찾다가 값도 정확히 나오는 함수가 완성되었다.

 

• d: 지구의 평균 반경
• e: 타원체의 이심률
• h, f: 변환 결과 영향 주는 보정 값
• c: 스케일 팩터
• l: 기준 경도
• m: 기준 위도
• lat, lon: 변환하려는 지점

 

위도와 경도를 라디안 단위로 변환하고, 타원체 관련 계산을 통해 중간 값을 얻어 새로운 좌표로 변환한다.

타원체라 이심률 (e)가 1보다 크면 역수를 써서 타원체를 고려해야 한다.

 

수학 잘하시는 분이 계시다면 자세한 설명 추가 부탁드립니다.

(TM은 Transverse Mercator 좌표계를 의미함, W콩나물에서 WGS84 변환은 반대로 하면 됨)

// ConvertWGS84ToWCONGNAMUL converts coordinates from WGS84 to WCONGNAMUL.
func ConvertWGS84ToWCONGNAMUL(lat, long float64) WCONGNAMULCoord {
	x, y := transformWGS84ToKoreaTM(6378137, 0.0033528106647474805, 500000, 200000, 1, 38, 127, lat, long)
	x = math.Round(x * 2.5)
	y = math.Round(y * 2.5)
	return WCONGNAMULCoord{X: x, Y: y}
}

// transformWGS84ToKoreaTM conversion
func transformWGS84ToKoreaTM(d, e, h, f, c, l, m, lat, lon float64) (float64, float64) {
	A := math.Pi / 180
	latRad := lat * A
	lonRad := lon * A
	lRad := l * A
	mRad := m * A

	w := 1 / e
	if e > 1 {
		w = e
	}

	z := d * (w - 1) / w
	G := 1 - (z*z)/(d*d)
	w = (d*d - z*z) / (z * z)
	z = (d - z) / (d + z)

	E := d * (1 - z + 5*(z*z-z*z*z)/4 + 81*(z*z*z*z-z*z*z*z*z)/64)
	I := 3 * d * (z - z*z + 7*(z*z*z-z*z*z*z)/8 + 55*z*z*z*z*z/64) / 2
	J := 15 * d * (z*z - z*z*z + 3*(z*z*z*z-z*z*z*z*z)/4) / 16
	L := 35 * d * (z*z*z - z*z*z*z + 11*z*z*z*z*z/16) / 48
	M := 315 * d * (z*z*z*z - z*z*z*z*z) / 512

	D := lonRad - mRad
	u := E*lRad - I*math.Sin(2*lRad) + J*math.Sin(4*lRad) - L*math.Sin(6*lRad) + M*math.Sin(8*lRad)
	z = u * c
	sinLat := math.Sin(latRad)
	cosLat := math.Cos(latRad)
	t := sinLat / cosLat
	G = d / math.Sqrt(1-G*sinLat*sinLat)

	u = E*latRad - I*math.Sin(2*latRad) + J*math.Sin(4*latRad) - L*math.Sin(6*latRad) + M*math.Sin(8*latRad)
	o := u * c

	E = G * sinLat * cosLat * c / 2
	I = G * sinLat * math.Pow(cosLat, 3) * c * (5 - t*t + 9*w + 4*w*w) / 24
	J = G * sinLat * math.Pow(cosLat, 5) * c * (61 - 58*t*t + t*t*t*t + 270*w - 330*t*t*w + 445*w*w + 324*w*w*w - 680*t*t*w*w + 88*w*w*w*w - 600*t*t*w*w*w - 192*t*t*w*w*w*w) / 720
	H := G * sinLat * math.Pow(cosLat, 7) * c * (1385 - 3111*t*t + 543*t*t*t*t - t*t*t*t*t*t) / 40320
	o += D*D*E + D*D*D*D*I + D*D*D*D*D*D*J + D*D*D*D*D*D*D*D*H
	y := o - z + h

	o = G * cosLat * c
	z = G * math.Pow(cosLat, 3) * c * (1 - t*t + w) / 6
	w = G * math.Pow(cosLat, 5) * c * (5 - 18*t*t + t*t*t*t + 14*w - 58*t*t*w + 13*w*w + 4*w*w*w - 64*t*t*w*w - 25*t*t*w*w*w) / 120
	u = G * math.Pow(cosLat, 7) * c * (61 - 479*t*t + 179*t*t*t*t - t*t*t*t*t*t) / 5040
	x := f + D*o + D*D*D*z + D*D*D*D*D*w + D*D*D*D*D*D*D*u

	return x, y
}

 

테스트를 잘 통과한다.