Go: WGS84를 WCONGNAMUL로 변환 함수

컴퓨터/Go language 2024. 3. 24. 19:44

728x90

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

Cut-off ratio for small singular values of a. For the purposes of rank determination, singular values are treated as zero if they are smaller than rcond times the largest singular value of a. Changed in version 1.14.0: If not set, a FutureWarning is given.

numpy.org

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
}

테스트를 잘 통과한다.

번외

WCONGNAMUL을 다시 WGS84로 변환하기 (margin error ~= 0.0001)

const (
	// Constants related to the WGS84 ellipsoid.
	aWGS84           float64 = 6378137 // Semi-major axis.
	flatteningFactor float64 = 0.0033528106647474805

	// Constants for Korea TM projection.
	k0               float64 = 1      // Scale factor.
	dx               float64 = 500000 // False Easting.
	dy               float64 = 200000 // False Northing.
	lat0             float64 = 38     // Latitude of origin.
	lon0             float64 = 127    // Longitude of origin.
	radiansPerDegree float64 = math.Pi / 180
)

// ConvertWCONGToWGS84 translates WCONGNAMUL coordinates to WGS84.
func ConvertWCONGToWGS84(x, y float64) (float64, float64) {
	return transformKoreaTMToWGS84(aWGS84, flatteningFactor, dx, dy, k0, lat0, lon0, x/2.5, y/2.5)
}

// transformKoreaTMToWGS84 transforms coordinates from Korea TM to WGS84.
// a, f, falseEasting, falseNorthing, scaleFactor, latOrigin, lonOrigin
func transformKoreaTMToWGS84(d, e, h, f, c, l, m, x, y float64) (float64, float64) {
	u := e
	if u > 1 {
		u = 1 / u
	}
	w := math.Atan(1) / 45 // Conversion factor from degrees to radians
	o := l * w
	D := m * w
	u = 1 / u
	B := d * (u - 1) / u
	z := (d*d - B*B) / (d * d)
	u = (d*d - B*B) / (B * B)
	B = (d - B) / (d + B)

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

	o = G*o - E*math.Sin(2*o) + I*math.Sin(4*o) - J*math.Sin(6*o) + L*math.Sin(8*o)
	o *= c
	o = y + o - h
	M := o / c
	H := d * (1 - z) / math.Pow(math.Sqrt(1-z*math.Pow(math.Sin(0), 2)), 3)
	o = M / H
	for i := 0; i < 5; i++ {
		B = G*o - E*math.Sin(2*o) + I*math.Sin(4*o) - J*math.Sin(6*o) + L*math.Sin(8*o)
		H = d * (1 - z) / math.Pow(math.Sqrt(1-z*math.Pow(math.Sin(o), 2)), 3)
		o += (M - B) / H
	}
	H = d * (1 - z) / math.Pow(math.Sqrt(1-z*math.Pow(math.Sin(o), 2)), 3)
	G = d / math.Sqrt(1-z*math.Pow(math.Sin(o), 2))
	B = math.Sin(o)
	z = math.Cos(o)
	E = B / z
	u *= z * z
	A := x - f
	B = E / (2 * H * G * math.Pow(c, 2))
	I = E * (5 + 3*E*E + u - 4*u*u - 9*E*E*u) / (24 * H * G * G * G * math.Pow(c, 4))
	J = E * (61 + 90*E*E + 46*u + 45*E*E*E*E - 252*E*E*u - 3*u*u + 100*u*u*u - 66*E*E*u*u - 90*E*E*E*E*u + 88*u*u*u*u + 225*E*E*E*E*u*u + 84*E*E*u*u*u - 192*E*E*u*u*u*u) / (720 * H * G * G * G * G * G * math.Pow(c, 6))
	H = E * (1385 + 3633*E*E + 4095*E*E*E*E + 1575*E*E*E*E*E*E) / (40320 * H * G * G * G * G * G * G * G * math.Pow(c, 8))
	o = o - math.Pow(A, 2)*B + math.Pow(A, 4)*I - math.Pow(A, 6)*J + math.Pow(A, 8)*H
	B = 1 / (G * z * c)
	H = (1 + 2*E*E + u) / (6 * G * G * G * z * z * z * math.Pow(c, 3))
	u = (5 + 6*u + 28*E*E - 3*u*u + 8*E*E*u + 24*E*E*E*E - 4*u*u*u + 4*E*E*u*u + 24*E*E*u*u*u) / (120 * G * G * G * G * G * z * z * z * z * z * math.Pow(c, 5))
	z = (61 + 662*E*E + 1320*E*E*E*E + 720*E*E*E*E*E*E) / (5040 * G * G * G * G * G * G * G * z * z * z * z * z * z * z * math.Pow(c, 7))
	A = A*B - math.Pow(A, 3)*H + math.Pow(A, 5)*u - math.Pow(A, 7)*z
	D += A

	return o / w, D / w // LATITUDE, LONGITUDE
}

Python

pyproj를 쓰면 쉽게 변환할 수 있다.

from pyproj import CRS, Transformer

# EPSG:4326 (WGS84)와 EPSG:5181 (카카오 좌표계) 정의
crs_wgs84 = CRS.from_epsg(4326)
crs_kakao = CRS.from_proj4("+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=GRS80 +units=m +no_defs")

# Transformer 객체 생성
transformer = Transformer.from_crs(crs_wgs84, crs_kakao, always_xy=True)

# 테스트 좌표 (위도, 경도)
lat, lon = 35.73294563400083, 127.37264182214031

# 변환 수행
x_kakao, y_kakao = transformer.transform(lon, lat)
x_kakao, y_kakao = x_kakao * 2.5, y_kakao * 2.5
print("WGS84 to Kakao EPSG:5181 -> X:", x_kakao, "Y:", y_kakao)

# 예상 결과와 비교
expected_x_kakao = 584279.0
expected_y_kakao = 621193.0
print("Expected -> X:", expected_x_kakao, "Y:", expected_y_kakao)

# 결과 비교
if abs(x_kakao - expected_x_kakao) < 1 and abs(y_kakao - expected_y_kakao) < 1:
    print("변환 결과가 예상 값과 일치합니다.")
else:
    print("변환 결과가 예상 값과 일치하지 않습니다.")

728x90

저작자표시 비영리 변경금지

'컴퓨터 > Go language' 카테고리의 다른 글

Go: 백엔드 웹소켓 채팅 (fiber v2 + websocket + rabbitmq) (0)	2024.03.31
Go: fiber v2 백엔드에서 토스 페이 API 사용하기 (0)	2024.03.04
dismember: 메모리 탐색기 (0)	2022.07.06

인기포스트 MORE POST

ABOUT ME

STUDY BITS STUDY BITS

계수 테스트 (Go언어)

번외

Python

'컴퓨터 > Go language' 카테고리의 다른 글

티스토리툴바

인기포스트 MORE POST

ABOUT ME

계수 테스트 (Go언어)

번외

Python

'컴퓨터 > Go language' 카테고리의 다른 글

관련글 관련글 더보기

티스토리툴바