위도, 경도를 격자 좌표로 변환하기

개발

위도, 경도를 격자 좌표로 변환하기

Beretta 2024. 11. 24. 09:43

개요

https://gist.github.com/fronteer-kr/14d7f779d52a21ac2f16

기상청 격자 <-> 위경도 변환

GitHub Gist: instantly share code, notes, and snippets.

gist.github.com

위 코드를 참조하면 되지만, 그냥 가져다 쓰기에는 모르는 내용이 너무 많아 이것저것 찾아보았다

상세

1. LCC DFS 좌표 변환의 의미

LCC는 람베르트 정각원추도법이라고 한다.

https://ko.wikipedia.org/wiki/%EB%9E%8C%EB%B2%A0%EB%A5%B4%ED%8A%B8_%EC%A0%95%EA%B0%81%EC%9B%90%EC%B6%94%EB%8F%84%EB%B2%95

람베르트 정각원추도법 - 위키백과, 우리 모두의 백과사전

위키백과, 우리 모두의 백과사전. 람베르트 정각원추도법(람베르트正角圓錐圖法, Lambert conformal conic projection)은 정각도법인 원추도법의 일종이다. 독일 사람 요한 하인리히 람베르트가 고안한

ko.wikipedia.org

DFS는 digital forecast system의 약자로 보인다.

2. 각 상수들의 의미

    var RE = 6371.00877; // 지구 반경(km)
    var GRID = 5.0; // 격자 간격(km)
    var SLAT1 = 30.0; // 투영 위도1(degree)
    var SLAT2 = 60.0; // 투영 위도2(degree)
    var OLON = 126.0; // 기준점 경도(degree)
    var OLAT = 38.0; // 기준점 위도(degree)
    var XO = 43; // 기준점 X좌표(GRID)
    var YO = 136; // 기1준점 Y좌표(GRID)

GRID: 격자 간격. 치환된 좌표의 1단위가 몇 km을 표현하는지 나타낸다. 만약 x=1과 x=2의 간격을 1m로 두고 싶으면 GRID 값을 0.001로 두면 된다.
SLAT1, SLAT2: 람베르트 정각원추도법에서 사용하는 두 개의 표준 위도. 이 위도 사이에서 정확도가 제일 높다고 한다. 따라서 어떤 지역에서 사용하느냐에 따라 값을 변경하면 되며, 기상청에서는 30도~60도의 중위도 지역을 넓게 잡았다.
OLON, OLAT: 기준점 위경도. 이 위경도 근처에서 가장 정확하다.
XO, YO: 기준점 XY좌표. 그런데 OLON, OLAT만 있으면 사실상 캐싱의 의미밖에 없을 듯.

ChatGPT를 이용해서 코드를 C#으로 변환하기

링크에도 다양한 예제가 있지만, 딱 원하는 코드가 없어서 ChatGPT에게 JavaScript를 주고 컨버팅을 시켰다.

이 코드를 C#으로 재작성해줘. 이때, 다음과 같은 지시사항을 따라줘
1. 변수명을 줄여쓰지 않고 있는 그대로 써줘. 예를들면 lat는 latitude로 바꿔야겠지.
2. public static class로 만들고 위경도를 XY로 변환하는 함수와, XY를 위경도로 변환하는 static 함수 두 개를 만들어줘
3. 기준점 경도, 기준점 위도는 static 함수의 인자로 포함시켜줘. 내 위경도를 기준으로 좌표를 펼치려고 해.

1단위좌표를 1미터로 두기 위해 GridSpacing을 0.001로 수정했다.

public static class LCCConverter
{
    // 상수 정의
    private const double EarthRadius = 6371.00877; // 지구 반경 (km)
    private const double GridSpacing = 0.001; // 격자 간격 (km)
    private const double StandardLatitude1 = 30.0; // 투영 위도 1 (degree)
    private const double StandardLatitude2 = 60.0; // 투영 위도 2 (degree)

    public static (double X, double Y) ConvertToXY(double userLatitude, double userLongitude, double targetLatitude, double targetLongitude)
    {
        const double DegreeToRadian = Math.PI / 180.0;

        var re = EarthRadius / GridSpacing;
        var standardLat1Radian = StandardLatitude1 * DegreeToRadian;
        var standardLat2Radian = StandardLatitude2 * DegreeToRadian;
        var originLongitudeRadian = userLongitude * DegreeToRadian;
        var originLatitudeRadian = userLatitude * DegreeToRadian;

        var scaleFactor = Math.Tan(Math.PI * 0.25 + standardLat2Radian * 0.5) / Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactor = Math.Log(Math.Cos(standardLat1Radian) / Math.Cos(standardLat2Radian)) / Math.Log(scaleFactor);

        var scaleFactorAtOrigin = Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactorAtOrigin = Math.Pow(scaleFactorAtOrigin, scaleFactor) * Math.Cos(standardLat1Radian) / scaleFactor;

        var radiusAtOrigin = Math.Tan(Math.PI * 0.25 + originLatitudeRadian * 0.5);
        radiusAtOrigin = re * scaleFactorAtOrigin / Math.Pow(radiusAtOrigin, scaleFactor);

        var radius = Math.Tan(Math.PI * 0.25 + targetLatitude * DegreeToRadian * 0.5);
        radius = re * scaleFactorAtOrigin / Math.Pow(radius, scaleFactor);

        var angle = targetLongitude * DegreeToRadian - originLongitudeRadian;
        if (angle > Math.PI) angle -= 2.0 * Math.PI;
        if (angle < -Math.PI) angle += 2.0 * Math.PI;
        angle *= scaleFactor;

        var x = radius * Math.Sin(angle);
        var y = radiusAtOrigin - radius * Math.Cos(angle);

        return (X: x, Y: y);
    }

    public static (double Latitude, double Longitude) ConvertToLL(double userLatitude, double userLongitude, double x, double y)
    {
        const double DegreeToRadian = Math.PI / 180.0;
        const double RadianToDegree = 180.0 / Math.PI;

        var re = EarthRadius / GridSpacing;
        var standardLat1Radian = StandardLatitude1 * DegreeToRadian;
        var standardLat2Radian = StandardLatitude2 * DegreeToRadian;
        var originLongitudeRadian = userLongitude * DegreeToRadian;
        var originLatitudeRadian = userLatitude * DegreeToRadian;

        var scaleFactor = Math.Tan(Math.PI * 0.25 + standardLat2Radian * 0.5) / Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactor = Math.Log(Math.Cos(standardLat1Radian) / Math.Cos(standardLat2Radian)) / Math.Log(scaleFactor);

        var scaleFactorAtOrigin = Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactorAtOrigin = Math.Pow(scaleFactorAtOrigin, scaleFactor) * Math.Cos(standardLat1Radian) / scaleFactor;

        var radiusAtOrigin = Math.Tan(Math.PI * 0.25 + originLatitudeRadian * 0.5);
        radiusAtOrigin = re * scaleFactorAtOrigin / Math.Pow(radiusAtOrigin, scaleFactor);

        var deltaX = x;
        var deltaY = radiusAtOrigin - y;
        var radius = Math.Sqrt(deltaX * deltaX + deltaY * deltaY);
        if (scaleFactor < 0.0) radius = -radius;

        var latitude = Math.Pow((re * scaleFactorAtOrigin / radius), (1.0 / scaleFactor));
        latitude = 2.0 * Math.Atan(latitude) - Math.PI * 0.5;

        double angle;
        if (Math.Abs(deltaX) <= 0.0)
        {
            angle = 0.0;
        }
        else
        {
            angle = Math.Atan2(deltaX, deltaY);
        }

        var longitude = angle / scaleFactor + originLongitudeRadian;
        return (Latitude: latitude * RadianToDegree, Longitude: longitude * RadianToDegree);
    }

    public static double CalculateDistance(double userLatitude, double userLongitude, double lat1, double lon1, double lat2, double lon2)
    {
        // 두 위경도를 LCC 격자 좌표로 변환
        var (x1, y1) = ConvertToXY(userLatitude, userLongitude, lat1, lon1);
        var (x2, y2) = ConvertToXY(userLatitude, userLongitude, lat2, lon2);

        // 두 격자 좌표 간의 거리 계산
        var distance = Math.Sqrt(Math.Pow(x2 - x1, 2) + Math.Pow(y2 - y1, 2)) * GridSpacing * 1000; // km -> m
        return distance;
    }
}

(한번에 딱 되지는 않았고 몇 번 지적해줘서 만들었다.)

ChatGPT 코드 검산하기

이 코드가 제대로 된 코드인지 알아보기 위해서, 위경도를 네이버 지도에 입력해보고, 거리를 재본다.

ChatGPT에게 명동 성당 중심으로 서쪽으로 150m, 동쪽으로 150m 되는 지점 2개를 잡아달라고 했다.

네이버 지도 검색에서 위경도를 넣으면 그 지점을 정확히 집어준다.

두 좌표를 네이버에 넣고 거리를 재보면...

정확히 300m 차이가 나는 좌표를 집어줬다. 이걸 LCCConverter에 집어넣고 제대로 된 거리를 주는지 보자.

약간의 오차가 있다.

표준 위도의 설정에 따라 오차가 발생할 수 있으므로, 한국에서만 사용한다고 가정하여 위경도를 바꿔본다.

ChatGPT가 추천해준 값을 적용해서 다시 계산하면 다음과 같이 나온다.

그런데, 이미 코드에서 userLatitude, userLongitude를 받고 있으므로, 이것을 기준으로 ±2 정도의 마진을 주면 될 것 같다.

내가 사용할 시스템은 나의 위치를 중심으로 한 좌표계를 구성할 것이기 때문이다.

다음은 userLatitude로부터 standardLatitude를 생성하도록 바꾼 코드이다.

public static class LCCConverter
{
    // 상수 정의
    private const double EarthRadius = 6371.00877; // 지구 반경 (km)
    private const double GridSpacing = 0.001; // 격자 간격 (km)
    private const double StandardLatitudePadding = 2;

    public static (double X, double Y) ConvertToXY(double userLatitude, double userLongitude, double targetLatitude, double targetLongitude)
    {
        const double DegreeToRadian = Math.PI / 180.0;

        var re = EarthRadius / GridSpacing;
        var standardLat1 = userLatitude - StandardLatitudePadding;
        var standardLat2 = userLatitude + StandardLatitudePadding;
        var standardLat1Radian = standardLat1 * DegreeToRadian;
        var standardLat2Radian = standardLat2 * DegreeToRadian;
        var originLongitudeRadian = userLongitude * DegreeToRadian;
        var originLatitudeRadian = userLatitude * DegreeToRadian;

        var scaleFactor = Math.Tan(Math.PI * 0.25 + standardLat2Radian * 0.5) / Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactor = Math.Log(Math.Cos(standardLat1Radian) / Math.Cos(standardLat2Radian)) / Math.Log(scaleFactor);

        var scaleFactorAtOrigin = Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactorAtOrigin = Math.Pow(scaleFactorAtOrigin, scaleFactor) * Math.Cos(standardLat1Radian) / scaleFactor;

        var radiusAtOrigin = Math.Tan(Math.PI * 0.25 + originLatitudeRadian * 0.5);
        radiusAtOrigin = re * scaleFactorAtOrigin / Math.Pow(radiusAtOrigin, scaleFactor);

        var radius = Math.Tan(Math.PI * 0.25 + targetLatitude * DegreeToRadian * 0.5);
        radius = re * scaleFactorAtOrigin / Math.Pow(radius, scaleFactor);

        var angle = targetLongitude * DegreeToRadian - originLongitudeRadian;
        if (angle > Math.PI) angle -= 2.0 * Math.PI;
        if (angle < -Math.PI) angle += 2.0 * Math.PI;
        angle *= scaleFactor;

        var x = radius * Math.Sin(angle);
        var y = radiusAtOrigin - radius * Math.Cos(angle);

        return (X: x, Y: y);
    }

    public static (double Latitude, double Longitude) ConvertToLL(double userLatitude, double userLongitude, double x, double y)
    {
        const double DegreeToRadian = Math.PI / 180.0;
        const double RadianToDegree = 180.0 / Math.PI;

        var re = EarthRadius / GridSpacing;
        var standardLat1 = userLatitude - StandardLatitudePadding;
        var standardLat2 = userLatitude + StandardLatitudePadding;
        var standardLat1Radian = standardLat1 * DegreeToRadian;
        var standardLat2Radian = standardLat2 * DegreeToRadian;
        var originLongitudeRadian = userLongitude * DegreeToRadian;
        var originLatitudeRadian = userLatitude * DegreeToRadian;

        var scaleFactor = Math.Tan(Math.PI * 0.25 + standardLat2Radian * 0.5) / Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactor = Math.Log(Math.Cos(standardLat1Radian) / Math.Cos(standardLat2Radian)) / Math.Log(scaleFactor);

        var scaleFactorAtOrigin = Math.Tan(Math.PI * 0.25 + standardLat1Radian * 0.5);
        scaleFactorAtOrigin = Math.Pow(scaleFactorAtOrigin, scaleFactor) * Math.Cos(standardLat1Radian) / scaleFactor;

        var radiusAtOrigin = Math.Tan(Math.PI * 0.25 + originLatitudeRadian * 0.5);
        radiusAtOrigin = re * scaleFactorAtOrigin / Math.Pow(radiusAtOrigin, scaleFactor);

        var deltaX = x;
        var deltaY = radiusAtOrigin - y;
        var radius = Math.Sqrt(deltaX * deltaX + deltaY * deltaY);
        if (scaleFactor < 0.0) radius = -radius;

        var latitude = Math.Pow((re * scaleFactorAtOrigin / radius), (1.0 / scaleFactor));
        latitude = 2.0 * Math.Atan(latitude) - Math.PI * 0.5;

        double angle;
        if (Math.Abs(deltaX) <= 0.0)
        {
            angle = 0.0;
        }
        else
        {
            angle = Math.Atan2(deltaX, deltaY);
        }

        var longitude = angle / scaleFactor + originLongitudeRadian;
        return (Latitude: latitude * RadianToDegree, Longitude: longitude * RadianToDegree);
    }

    public static double CalculateDistance(double userLatitude, double userLongitude, double lat1, double lon1, double lat2, double lon2)
    {
        // 두 위경도를 LCC 격자 좌표로 변환
        var (x1, y1) = ConvertToXY(userLatitude, userLongitude, lat1, lon1);
        var (x2, y2) = ConvertToXY(userLatitude, userLongitude, lat2, lon2);

        // 두 격자 좌표 간의 거리 계산
        var distance = Math.Sqrt(Math.Pow(x2 - x1, 2) + Math.Pow(y2 - y1, 2)) * GridSpacing * 1000; // km -> m
        return distance;
    }
}

결과는 조금 더 정확해졌다.

'개발' 카테고리의 다른 글

Unity에서 System.Text.Json 사용하기 (0)	2024.11.24
Unity에서 외부 폴더의 코드를 참조하기 (0)	2024.11.24

현재글위도, 경도를 격자 좌표로 변환하기

베레타 블로그

개발 관련 히스토리 정리하는 블로그입니다.

Today :
Yesterday :

베레타 블로그