가끔 folium의 지역위주의 지도표시를 넘어 국가적, 지구적인 지도(그림파일)에 데이터 표시가 필요할 때가 있습니다.
특히나 지구본 처럼 구 형태에 시각화를 하고 싶을 때, 바로 그때 유용한 패키지가 basemap 인데요.
저는 이번에 지구본에서 최단경로로 움직이는 비행경로(측지선 geodesic)를 시각화해 볼게요!!!
물론 자유 소프트웨어 파이썬을 이용합니다! 충성!! 충성!!!
먼저 좋은 아이디어와 수학적 증명은 아래 블로그에서 힌트를 많이 얻었습니다. 감사드립니다!
C/C++ 반복-이완 계산법으로 알아보는 비행경로 (tistory.com)
여기서는 핵심만 요약해서 간략히 정리할 예정입니다.
1. 지구 위에 2 지점(위경도) 사이에 최소거리는 ?
다양한 해법이 있겠지만, 저는 이 방법이 제일 간단하고 이해도 쉬운 것 같아요.
바로 '구좌표계에서 직각좌표계 변환', '벡터의 내적' 2가지를 이용해 최단거리를 계산할거에요!
물론 이때 계산의 편의를 위해 지구는 완전한 구형이라고 전제 합니다.
구형이 아니라면 계산이 많이 복잡해질테니까요.
아마 일반적인 측지선방정식으로 확장하면 풀이할 수 있을 것 같은데 이건 정말 제 능력밖이네요.ㅠㅠ
2. 2 지점 사이 최단 비행경로의 궤적은 ??
(th=위도, ph=경도 로 두면,) 먼저 두 지점사이의 거리식은 L = R * ∫√(th'**2 + (cos(th)*ph')**2)dt
정상적이 되어야 하는 정적분으로 부터 오일러-라그랑주 방정식을 구하고
이렇게 구한 상미분방정식(th'' = -cos(th)*sin(th)*ph'**2 , ph'' = 2*tan(th)*th'*ph')을
적절한 경계조건의 초기값을 주어 Runge-Kutta 방법을 사용해서 풀 수 있습니다.
이렇게 구한 해는 우리가 직관적으로 이해하고 있는 이른바 대원(great circle)이 됩니다.
그럼,,, 이제 지구본에서 (인천공항-런던히드로공항 간) 최단 비행경로를 시각화하는 소스코드를 볼게요.
from numpy import *
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
R = 6400
def ang2rad(x) : return x*pi/180
def rad2ang(x) : return x*180/pi
def s2r(r, th, ph) : return r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)
def ina(v0, v1) : # included angle
return arccos(dot(v0, v1)/sqrt(dot(v0, v0)*dot(v1, v1)))
ICN = [37.46, 126.44]
LHR = [51.47, -0.45]
v0 = array([0, 0, 1])
v1 = array(s2r(1, ang2rad(90 - ICN[0]), ang2rad(ICN[1])))
v2 = array(s2r(1, ang2rad(90 - LHR[0]), ang2rad(LHR[1])))
print(ina(v1, v2)*R)
fig=plt.figure()
ax=fig.add_axes([0.1,0.1,0.8,0.8])
m = Basemap(llcrnrlon=-30.,llcrnrlat=25.,urcrnrlon=150.,urcrnrlat=75.,
rsphere=(6378137.00,6356752.3142),
resolution='l',projection='merc',
lat_0=30.,lon_0=-20.,lat_ts=20.)
gcc = m.drawgreatcircle(ICN[1],ICN[0],LHR[1],LHR[0], linewidth=2, color='None')
xys = gcc[0].get_xydata()
m.plot(xys[:,0], xys[:,1], 'g')
m.drawcoastlines()
m.fillcontinents()
m.drawparallels(arange(10,90,10),labels=[1,1,0,1])
m.drawmeridians(arange(-180,180,20),labels=[1,1,0,1])
ax.set_title('Great Circle from Incheon to London')
fig=plt.figure()
map = Basemap(projection='ortho',lat_0=40,lon_0=80,resolution='l')
map.drawgreatcircle(ICN[1],ICN[0],LHR[1],LHR[0], linewidth=2, color='b')
map.drawcoastlines(linewidth=0.25)
map.drawcountries(linewidth=0.25)
map.fillcontinents(color='coral',lake_color='aqua')
map.drawmapboundary(fill_color='aqua')
map.drawmeridians(arange(0,360,30))
map.drawparallels(arange(-90,90,30))
plt.title('Great Circle from Incheon to London')
plt.show(block=False)
지구의 반지름이 6400km 인데, ICH-LHR 간 최단경로 거리는 약 8900km 정도네요. (자동차 6개월 주행거리 정도구나!)
항공기가 이 최단경로를 시속 약 1000km 속도로 구름 위에서 하늘을 가르며 날고 있다고 상상해 보세요!
크~ 멋집니다ㅎㅎ
좋고 멋있는 것들이 거의 전부 영국에서 유래된 것이 많은 것 같아요.
영어, 셰익스피어, 그리니치천문대, 뉴턴, 시민혁명, 산업혁명 등등,,, 하다못해 게임에도 최애유닛 장궁병...
그래서인지 더,,, 오래전 여행 갔었던 영국에 아련하고 좋은 추억이 떠오르네요.
이대로 끝내기 아쉬워 몇장의 베스트 컷을 공유합니다^^;
----- 2023.4.17. 추가 -----
3. 허블 우주망원경 위성궤도는 어떻게 그릴까??
미국은 과거에 특정한 궤도경사각(극궤도 90˚, 정지궤도 0˚ 등)이 필요한 경우를 제외하곤, 일반적인 우주선은 케이프커너배럴(N28.45˚, W80.68)에서 정동쪽으로 로켓을 발사했습니다.(궤도경사각이 28.45˚) 아직 현역인 허블 우주망원경도 마찬가지고요. 왜냐면,,, 적도에 가까이에서 동쪽으로 로켓을 발사하면 지구 자전속도를 추가로 얻을 수 있어 위성궤도 진입을 위한 에너지를 더 낮출 수 있고, 동쪽이 망망대해 대서양이기 때문입니다.
참고로 허블 우주망원경은 지국자전방향과 동일한 태양동기궤도로 하루에 지구를 15회 돈다고 합니다.
기사 - [우주]우주선과 인공위성 궤도경사각, 그리고 발사장의 입지 (ddanzi.com) 를 참고해서 작성했습니다!!!
'파이썬' 카테고리의 다른 글
| 로또번호 분석 & 제조기 (1) | 2023.01.14 |
|---|---|
| 파이썬으로 별 이름표 붙이기(별주석) (0) | 2022.06.14 |
| GPIB를 이용한 네트워크 분석기 활용 (0) | 2022.05.08 |
| 파이썬 GUI 실행파일 만들기(전력변환기) (0) | 2021.12.03 |
| 파이썬으로 별자리판 만들기(1) (0) | 2021.01.15 |
