컴퓨터에 지구를 위치시키는 법 (feat. Geographic Coordinates)

컴퓨터로 지구 또는 지구안에 있는 무언가를 표현하려면 경도, 위도같은 좌표가 필요하다.

• Cartesian Coordinates (데카르트 좌표계)
• Geographic Coordinates (지리적 좌표계)

---

 

Geographic Coordinates (지리적 좌표계)

• Cesium에서는 Cartographic 객체
• spherical coordinate: (azimuth, inclination, radius) 로 표현 - (방위각, 경사도, 반경)
• geographic coordinates : (longitude, latitude, height) 로 표현 - (경도, 위도, 높이)
  • 경도: 서쪽에서 동쪽으로의 각도 (기준선 = 본초 자오선)
    • 본초 자오선 = 0º
    • 서쪽 = -180º
    • 동쪽 = +180º
  • 위도: 남쪽에서 북쪽으로의 각도 (기준선 = 적도)
    • 적도 = 0º
    • 남반구 = -90º
    • 북반구 = +90º
  • 높이: 표면 위 또는 아래로의 선형거리
• 대부분의 벡터 데이터는 Geographic 으로 표현
  • GPS 등에 Geographic 이용

---

 
경위도 $(0^{◦},{80}^{◦})$에서 $({10}^{◦},{90}^{◦})$의 면적과 $(0^{◦},0^{◦})$에서 $({10}^{◦},{10}^{◦})$을 비교해보자. 얼핏보면 각도의 제곱이 같기 때문에 면적이 같을 것 같다.

하지만 지구는 3차원 타원체다. 즉, 위치에 따라 편평률이 다르기 때문에 경위도 $(0^{◦},{80}^{◦})$에서 $({10}^{◦},{90}^{◦})$ 면적과 $(0^{◦},0^{◦})$에서 $({10}^{◦},{10}^{◦})$ 면적은 서로 다르다.
 
따라서 지리적 그리드 테셀레이션에서와 같이 극 근처의 균일한 경도/위도 그리드 오버샘플을 기반으로 하는 알고리즘이 필요하다. (나중에 설명)
 

사람들은 경, 위도에 익숙하다. 예를 들어, 대한민국 경위도는 127º, 37.5º 이다.
하지만 지구의 중심을 기준으로 대한민국의 데카르트 좌표는 다음과 같다.

 

X ≈ −3,044,804 m,
Y ≈ 4,043,821 m,
Z ≈ 3,884,000 m
 
사람들은 이런 표현에 익숙하지 않다. 하지만 컴퓨터는 이런 방식을 통해 무언가를 나타낸다.
대한민국 경위도를 데카르트좌표로 변환하는 방법을 알아보자.

---

 

원점을 중심으로 구의 방정식: $x_s^2+y_s^2+z_s^2=r^2$ (2.1)

 
1. 방정식 (2.1)을 만족하는 점($x_s,y_s,z_s)$는 구의 표면(surface)에 있다.
2. 구(sphere)에 있는 점의 경우, 점을 원점을 기준으로 벡터로 취급하고 정규화 하면 표면 법선을 구할 수 있다.
 
지구를 구로 표현하면 계산이 용이하겠지만, 지구는 타원체이기에 타원체 방정식을 사용하는 것이 더 정밀하고 유연하다.
 

타원체 방정식: $\frac{x_s^2}{a^2}+\frac{y_s^2}{b^2}+\frac{z_s^2}{c^2}=1$ (2.2)

 
1. 타원체는 각각 X축, Y축, Z축을 따라 세개의 반지름(a, b, c)로 정의된다.
2. 방정식 (2.2)를 만족하는 점($x_s,y_s,z_s)$는 타원체의 표면(surface)에 있다.
 

타원체의 한 점에서 지구 중심을 기준으로 법선을 그리면 이를 Geocentric이라고 한다.
하지만 타원체는 곡률이 존재하기 때문에, Geocentric은 대부분 점에 대해 표면에 수직이 아니다.
 
따라서 타원체의 한 점에서 지구 중심이 아닌 적도$(X,Y)$ 평면에 법선을 그려야 한다. 이를 Geodetic이라고한다.
Geodetic은 타원에 한 점에 대해 실제 법선이다.
 
구의 경우, Geodetic, Geocentric은 동일하다.
타원체는 대부분 점에 대해 Geodetic, Geocentric이 크게 다른데, 45º에서 가장 크게 다르다.

 
표면 법선 벡터(surface normal vector) = $\widehat{n_s}$ 구하기

$m=(\frac{x_s}{a^2},\frac{y_s}{b^2},\frac{z_s}{c^2})$
$\widehat{n_s}=\frac{m}{||m||}$

 

Geocentric 오차

지구 중심에서 특정 벡터 사이의 각도 : ${\varphi }_c$ = geocentric

지구 표면에서 특정 벡터 사이의 각도 : ${\varphi }_d$ = geodetic

위경도를 geocentric 또는 geodetic으로 구한 후, 높이(h)는 geodetic을 이용해서 구해야 한다.
 
높이(h)를 geocentric으로 구할 경우 우주와 같이 높이가 높은 곳에서 오차가 발생한다.
특히 geocentric, geodetic 각도 차이가 클 수록 오차가 커진다.

---

 

좌표 변환

Geographic coordinates → WGS84

실제 지리 데이터가 저장될 때는, 보통 경위도(지리적 좌표)로 저장된다.
 
컴퓨터에서 지구는 보통 WGS84 좌표계로 렌더링하기 때문에 경위도를 WGS84 좌표로 변환하는 작업은 필수적이다.
 
주의: WGS84 또한 Geographic coordinates 중 하나다.

• 전 세계에서 가장 널리 사용되는 지리 좌표계 (GPS 표준)
• 즉, WGS84는 Geographic Coordinates 중 하나일 뿐
• GCS(Geographic Coordinate System) ≠ WGS84
• WGS84 ⊂ GCS
• 원점: 지구 중심(Geocentric)
• 단위: 도(º)
• 좌표값: (위도, 경도, 타원체 고도)

---

 

 

$\widehat{n_s}=\cos \varphi \cos \lambda \hat{i}+\cos \varphi \sin \lambda \hat{j}+\sin \varphi \hat{k}$ (2.3)

구(Sphere)형 모델에서 유도된 형태이며, 완전한 타원체 모델을 표현하기 위해 추가 변환이 필요하다.
 

$\overrightarrow{n_s}=\frac{x_s}{a^2}\hat{i}+\frac{y_s}{b^2}\hat{j}+\frac{z_s}{c^2}\hat{k}$ (2.4)

타원체 표면에서의 법선 벡터를 나타내는 식  (단위벡터 아님)
 

$\widehat{n_s}=\gamma \overrightarrow{n_s}$ (2.5)

임의의 $\gamma$를 곱해 단위벡터를 찾을 수 있다.
 

$\widehat{n_s}=\gamma (\frac{x_s}{a^2}\hat{i}+\frac{y_s}{b^2}\hat{j}+\frac{z_s}{c^2}\hat{k})$ (2.6)

(2.4)의 법선 벡터를 (2.5)의 관계를 통해 단위 벡터로 변환하는 식
 

$\widehat{n_x}=\frac{\gamma x_s}{a^2}$
$\widehat{n_y}=\frac{\gamma y_s}{b^2}$ (2.7)
$\widehat{n_z}=\frac{\gamma z_s}{c^2}$

단위 법선 벡터의 성분을 타원체 좌표 $(x_s,y_s,z_s)$와 타원체 반지름 $(a,b,c)$로 나타낸 식
 

$x_s=\frac{a^2\widehat{n_x}}{\gamma }$
$y_s=\frac{b^2\widehat{n_y}}{\gamma }$ (2.8)
$z_s=\frac{c^2\widehat{n_z}}{\gamma }$

단위 법선 벡터 성분을 이용하여 타원체 표면의 좌표 $(x_s,y_s,z_s)$로 나타낸 식
 

타원체 방정식 $\frac{x_s^2}{a^2}+\frac{y_s^2}{b^2}+\frac{z_s^2}{c^2}=1$에 $(x_s,y_s,z_s)$ 대입

$\frac{(\frac{a^2\widehat{n_x}}{\gamma }{)}^2}{a^2}+\frac{(\frac{b^2\widehat{n_y}}{\gamma }{)}^2}{b^2}+\frac{(\frac{c^2\widehat{n_z}}{\gamma }{)}^2}{c^2}=1$

 
정리하면

$a^2{\widehat{n_x}}^2+b^2{\widehat{n_y}}^2+c^2{\widehat{n_z}}^2={\gamma }^2$ (2.9)
$\gamma =\sqrt{a^2{\widehat{n_x}}^2+b^2{\widehat{n_y}}^2+c^2{\widehat{n_z}}^2}$

• 정규화 상수 $\gamma$를 구하는 식

 
또한

$\overrightarrow{h}=h\widehat{n_s}$

$\overrightarrow{h}$ 는 타원체 표면에서 법선 방향으로의 변위 벡터
즉, $\overrightarrow{h}$ 는 타원체 표면의 한점에서 법선 벡터 방향으로의 거리 $h$ 만큼 떨어진 점의 위치 변화$h$ 는 타원체 표면에서 해당 점까지의 거리

 

$\because$ $r$ = $r_s$ + $h\widehat{n_s}$ (2.10)

• 최종적으로, 타원체 표면 $r_s$으로부터 높이 $h$만큼 떨어진 점의 위치를 결정하는 식

---

 

WGS84 → Geographic coordinates

일반적으로 WGS84에서 지리적 좌표로 변환하는 것은 반대 방향으로 변환하는 것 보다 더 많은 작업이 필요하므로 여러 단계로 나눈다.

$\widehat{n_s}=\frac{n_s}{||n_s||}$ (2.4 참고)
$\lambda =\arctan \frac{\widehat{n_y}}{\widehat{n_x}}$,
$\varphi =arcsin\frac{\widehat{n_z}}{||n_s||}$

 

두 점 r0, r1을 Geocetnric, Geodetic 두 방식으로 나눠 크기를 조정할 수 있다.
 
(a) : Geocentric Nomrmal을 따라 타원체 중심의 벡터가 타원체와 교차하여 표면지점을 결정
(b) : Geocentric Nomrmal을 따라 반복 작업(뉴턴-랩슨법 사용)을 통해 스케일링

뉴턴-랩슨법 특징
1. 점진적으로 정답을 찾아가는 근사 해법이라서 빠르게 수렴함.
2. 삼각함수 없이 계산 가능 ⇒ 효율적이다.

 
 
Geocentric Nomrmal (a) 특징

• 계산이 간단함 (직선 방정식과 타원 방정식으로 해결 가능)
• 벡터 연산만 사용하므로 빠르게 계산 가능
• WGS84에서 지리 좌표계로 변환 했을 때 의 정확도를 보장할 수 없다.
  • 지구는 완벽한 구가 아니라 타원체 이기 때문.
• 따라서, Geocentric Nomrmal이 아닌 Geodetic Normal(b)을 사용해야 한다.

 
Geodetic Normal (b) 특징

• 정확도가 높음 (지구 타원체의 실제 표면을 고려함)
• WGS84 좌표 변환에서 가장 많이 사용됨
• 반복 연산이 필요함 (하지만 뉴턴 랩슨법을 통해 빠르게 수렴함)

$\overrightarrow{r}=r-\overrightarrow{0}$
$r=(x,y,z)$ 지구 중심에서의 위치를 의미
$\overrightarrow{r_s}=\beta \overrightarrow{r}$
$\beta$ = 크기 조정(scaling) 계수로, $\overrightarrow{r}$을 조정하여 표면에 맞춤
$\beta =\frac{1}{\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}}}$ (2.11)

• $\overrightarrow{r_s}$ 의미
  • 주어진 점 $(x,y,z)$이 지구 중심에서 어디에 있든, 이 점을 타원체 표면으로 옮겨야 함
  • 즉, $(x,y,z)$ 를 적절한 비율로 축소(Scaling)하여 $(x_s,y_s,z_s)$를 만들고 싶음
  • 이때, 스케일링 계수를 β 로 두고 아래처럼 표현
    $x_s=\beta x$
    $y_s=\beta y$ (2.12)
    $z_s=\beta z$
  • 즉, β값을 곱해서 현재 점을 표면으로 투영하는 방식
• $\beta$ 의 역할
  • 현재 점을 지구 중심 방향으로 투영하여, 타원체 표면 위로 스케일링하는 역할
  • 만약 주어진 점이 지구 표면과 가깝다면, $\beta \approx 1$에 가까워짐
  • 만약 주어진 점이 지구 표면보다 훨씬 멀다면, $\beta$가 작아져서 점을 더 많이 축소함
• $\beta$ 값 유도 과정 (2.11)
  • $(x_s,y_s,z_s)$는 타원체 표면에 있으므로 타원체 방정식을 만족해야한다.${\beta }^2$ 값에 대해 정리하면 정리하면$\beta$ 값은 양수여야 하므로
  • $\beta =\frac{1}{\sqrt{\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}}}$ (2.11)
  • → ${\beta }^2(\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2})=1$
  • → $\frac{(\beta x{)}^2}{a^2}+\frac{(\beta y{)}^2}{b^2}+\frac{(\beta z{)}^2}{c^2}=1$

$\overrightarrow{r}=\overrightarrow{r_s}+\overrightarrow{h}$ (2.13)
$\overrightarrow{r}:$ 임의의 점 벡터 (구하고자 하는 점 위치)
$r_s:$ 표면에 있는 점 벡터
$\overrightarrow{h}:$ 표면 점으로부터 임의의 점 높이 벡터

• $\overrightarrow{h}$는 $\overrightarrow{n_s}$ 와 같은 방향을 가지지만 크기(magnitude)가 다르다.
  ⇒ $\overrightarrow{h}=\alpha \overrightarrow{n_s}$
  $\overrightarrow{n_s}:$ 지구 중심에서 표면점에 대해 정규화되지 않은 벡터
• $\overrightarrow{n_s}=\frac{x_s}{a^2}\hat{i}+\frac{y_s}{b^2}\hat{j}+\frac{z_s}{c^2}\hat{k}$ (2.14)

이 식을 (2.13)에 대입하면
$\overrightarrow{r}=\overrightarrow{r_s}+\alpha \overrightarrow{n_s}$
$\overrightarrow{r}$을 $\overrightarrow{r_s}$와 $\alpha \overrightarrow{n_s}$를 이용해서, 세 개의 스칼라 방정식( $x,y,z)$ 으로 나타내면

$x=x_s+\alpha \frac{x_s}{a^2},$
$y=y_s+\alpha \frac{y_s}{b^2},$
$z=z_s+\alpha \frac{z_s}{c^2}$

 
이 식을 인수 분해 하면

$x=x_s(1+\frac{\alpha }{a^2}),$
$y=y_s(1+\frac{\alpha }{b^2}),$
$z=z_s(1+\frac{\alpha }{c^2})$

 
이 식을 표면에 있는 점을 기준으로 재배열 하면

$x_s=\frac{x}{1+\frac{\alpha }{a^2}},$
$y_s=\frac{y}{1+\frac{\alpha }{b^2}},$ (2.15)
$z_s=\frac{z}{1+\frac{\alpha }{c^2}}$

즉, 알고 있는 점 $r=(x,y,z)$와 타원체의 반지름 $(a,b,c)$을 통해 $(x_s,y_s,z_s)$를 나타냈다.
식 (2.15) 에서 모르는 미지수는 $\alpha$ 이다.
 
즉, 미지수 $\alpha$만 찾으면 WGS84 좌표에서 임의의 점 $r$에 대해Geographic 좌표로 변환된 $\overrightarrow{r}$을 찾을 수 있다.
미지수 $\alpha$를 찾기 위해 타원체 방정식$(S)$ 에 $(x_s,y_s,z_s)$를 대입한다. (2.17)
 

$S=\frac{x_s^2}{a^2}+\frac{y_s^2}{b^2}+\frac{z_s^2}{c^2}-1=0$ (2.16)

 

$S=\frac{x^2}{a^2(1+\frac{\alpha }{a^2}{)}^2}+\frac{y^2}{b^2(1+\frac{\alpha }{b^2}{)}^2}+\frac{z^2}{c^2(1+\frac{\alpha }{c^2}{)}^2}-1=0$ (2.17)

${\alpha }_0=(1-\beta )\frac{||r||}{||n_s||}$

• $(1-\beta )$ 의미
  • $\beta$ 값은 (2.11 참고)
  • 현재 점이 표면에서 얼마나 떨어져 있는지를 나타내는 비율
  • $\beta$가 크면 현재 점이 표면에 가까워짐
  • $\beta$가 작으면 현재 점이 표면에서 멀리 떨어져짐
• $\frac{||r||}{||n_s||}$ 의미
  • 현재 좌표에서 법선 방향으로 얼마나 이동해야 하는지를 계산하는 기준
  • $||r||:$중심으로부터 임의의 점 $(x,y,z)$벡터 $\overrightarrow{r}$ 크기
    • 원래 주어진 WGS84 좌표의 벡터 크기
  • $||n_s||:$중심으로부터 지구 표면에 있는 점 $(x_s,y_s,z_s)$벡터 $\overrightarrow{n_s}$크기
    • 표면 좌표에서의 법선 벡터 크기
• 뉴턴 랩슨법
  • 초기 값 : ${\alpha }_0$
  • ${\alpha }_{n+1}={\alpha }_n-\frac{S({\alpha }_n)}{S^{\prime }({\alpha }_n)}$해당 작업을 통해 $\alpha$를 찾고, (2.17)에 대입하여 0과 근사해지는 값을 찾음.
    • 코드 예시
      • $\alpha$를 대입해 S 값이 0과 가까운 ${10}^{-10}$보다 작은 경우까지 반복
        
  • ⇒ $\alpha =\alpha -\frac{S}{\frac{\partial S}{\partial \alpha }}$
  • $\frac{\partial S}{\partial \alpha }=-2[\frac{x^2}{a^4(1+\frac{\alpha }{a^2}{)}^3}+\frac{y^2}{b^4(1+\frac{\alpha }{b^2}{)}^3}+\frac{z^2}{c^4(1+\frac{\alpha }{c^2}{)}^3}]$
• 타원체 방정식$(S)$과 뉴턴-랩슨법을 이용
• $S$ 와 $\frac{\partial S}{\partial \alpha }$를 평가하여 0에 가까워지는 $\alpha$를 찾는다.
• 즉 미지수 $\alpha$를 근사해, 현재 위치와 타원체 표면 사이의 거리를 근사적으로 계산한다.