BLUEnLIVE

구라제거기의 고향입니다

 
 

아래 사진은 SciTechDaily에 올라온 사진이다.
인터넷에서는 이처럼 ISS가 달 위를 지나가는 장면을 촬영한 사진을 가끔 볼 수 있다.

image

파이썬의 라이브러리들을 활용하면 ISS의 위치도, 달의 위치도 꽤 쉽게 알 수 있는데…
그렇다면 ISS가 달 위를 언제 지나가는지 알 수 있을까?

들어가기에 앞서

이전 포스팅들에서 ISS 및 달의 위치를 파악하는 방법을 간단히 알아보았다.
그 때 사용했던 라이브러리들을 활용하면 이 문제도 비교적 쉽게 풀 수 있을 것 같다.

달의 시직경은 약 0.5도이다.
ISS가 달 위를 지나간다는 건 관측자 위치에서 둘의 내각이 0.25도 이내라는 것을 의미한다.
그렇다면 달이 뜰 때부터 질 때까지의 모든 시간에 대해 이를 계산하면 되겠다.

그런데, 문제는 ISS의 위치다.
ISS는 TLE를 통해 위치를 알 수 있는데, TLE는 특성상 며칠 지나면 정확도가 떨어진다.
즉, 계산을 통해 알 수 있는 ISS의 위치는 하루 이상 지나면 신뢰도가 떨어진다는 뜻.

기존 코드들의 활용과 확장

ISS의 TLE를 얻는 방식은 이전 포스팅에서 다룬 것과 같다.

그 외에 달에 대한 정보를 알아내려면 astral 라이브러리를 사용하면 된다.
이 라이브러리를 사용하면 관측자의 위치를 기준으로 해와 달의 위치는 물론 뜨고 지는 시간도 알 수 있다.
날짜에 따라선 달이 뜨거나 지는 시간이 없는 경우가 있으므로 이를 고려해야 한다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

from skyfield.api import load
from datetime import timedelta
from astral import moon


def get_iss_object_skyfield():
    satellites = load.tle_file(url='https://live.ariss.org/iss.txt', reload=True)

    iss = satellites[0]

    return iss


def get_moon_rise_moon_set(location, tm):
    try:
        moon_rise = moon.moonrise(location.observer, tm)
    except:
        tm += timedelta(days=1)
        moon_rise = moon.moonrise(location.observer, tm)

    tm = moon_rise
    try:
        moon_set = moon.moonset(location.observer, tm)
    except:
        tm += timedelta(days=1)
        moon_set = moon.moonset(location.observer, tm)

    if moon_rise > moon_set:
        tm += timedelta(days=1)
        moon_set = moon.moonset(location.observer, tm)

    return moon_rise, moon_set

필요한 함수들의 구현

경위도 좌표를 직교좌표로 변환하는 함수가 우선 필요하다.
우리가 사용하는 모든 좌표는 WGS84 좌표계를 사용하는 점을 고려해야 한다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import numpy as np

def wgs84_llh_to_xyz(llh):
    lat_rad = np.deg2rad(llh[0])
    lon_rad = np.deg2rad(llh[1])
    alt = llh[2] if len(llh) > 2 else 0

    clat = np.cos(lat_rad)
    slat = np.sin(lat_rad)
    clon = np.cos(lon_rad)
    slon = np.sin(lon_rad)

    wgs84_a = 6378137.0
    wgs84_finv = 298.257223563
    wgs84_f = 1 / wgs84_finv
    wgs84_e2 = 1 - (1 - wgs84_f) * (1 - wgs84_f)

    rn = wgs84_a / np.sqrt(1 - wgs84_e2 * slat * slat)

    x = (rn + alt) * clat * clon
    y = (rn + alt) * clat * slon
    z = (rn * (1 - wgs84_e2) + alt) * slat

    return np.array([x, y, z])

우리는 ISS의 위치를 직교좌표로 알 수 있다.
그런데, 달의 위치는 관측자의 위치를 기준으로 방위와 천정각을 알 수 있다.
둘을 비교하려면 ISS의 위치를 관측자의 위치 기준의 방위와 천정각으로 변환해야 한다.

회전 변환을 위해 쿼터니언을 사용하면 편리하다.
쿼터니언을 사용하면 행렬로 표현하는 것보다 훨씬 간단하게 할 수 있다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

from pyquaternion import Quaternion
from astropy import units as u

def calc_azimuth_zenith(my_llh, iss_xyz):
    count = len(iss_xyz)
    ret = np.zeros((count, 2), dtype=np.float64)
    my_xyz = wgs84_llh_to_xyz(my_llh)
    lat = np.deg2rad(my_llh[0])
    lon = np.deg2rad(my_llh[1])

    vectors = iss_xyz.to(u.m).value - my_xyz
    # 이 vectors는 lat, lon에서 출발하는 벡터임
    # 이를 경위도 원점(0,0)에서 출발하는 벡터로 회전 변환 필요

    # iss를 경위도 0,0으로 회전시키는 쿼터니언 생성
    q = Quaternion(axis=[0, 1, 0], angle=lat) * Quaternion(axis=[0, 0, 1], angle=-lon)

    for i, v in enumerate(vectors):
        v = q.rotate(v)             # vector 회전
        v = v / np.linalg.norm(v)   # normalizing

        # azimuth, zenith 계산
        az = 180 - np.degrees(np.arctan2(v[1], -v[2]))
        ze = 90 - np.degrees(np.arcsin(v[0]))
        ret[i] = [az, ze]

    return ret

마지막으로 필요한 함수는 방위와 천정각이 주어졌을 때 두 점 사이의 각을 계산하는 함수다.
처음엔 크기가 1인 벡터를 사용해 만들었는데, copilot 님이 간단한 방식을 제안해주심.

1
2
3
4
5
6
7
8
9
10
11

def calc_angle_between(azimuth1, zenith1, azimuth2, zenith2):
    # 모든 각을 radian으로 변환
    azimuth1 = np.deg2rad(azimuth1)
    zenith1 = np.deg2rad(zenith1)
    azimuth2 = np.deg2rad(azimuth2)
    zenith2 = np.deg2rad(zenith2)

    cos_theta = np.sin(zenith1) * np.sin(zenith2) * np.cos(azimuth1 - azimuth2) + np.cos(zenith1) * np.cos(zenith2)
    theta = np.arccos(cos_theta)

    return np.rad2deg(theta)

본체 구현

실제로 구현해야 할 내용은 다음과 같다.

  • 주어진 날짜의 달이 뜨고 지는 시간 계산
  • 이 시간에 대해 달의 방위 계산
  • 이 시간에 대해 ISS의 위치 계산
  • ISS의 위치를 관측자의 위치 기준의 방위로 변환
  • 두 점 사이의 각 계산

즉, 시간 간격이 짧을수록 정확도는 높아지나 시간은 그만큼 더 소요되는 것이다.

일단 달이 뜨고 지는 시간을 계산하는 코드는 다음과 같다.
시간 간격은 0.5초로 하였다. 계산이 너무 오래 걸리면 이 값을 늘리면 된다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

from datetime import datetime
from pytz import timezone, utc
import matplotlib.pyplot as plt
from astral import LocationInfo

observer_lat = 37.46707105652381
observer_lon = 127.07507550716402
observer_alt = 73.0
step_in_second = 0.5    # 0.5초 간격으로 달이 뜨는 시각부터 지는 시각까지의 시각을 구한다.
days_after = 0          # 며칠 뒤를 계산할 것인지

location = LocationInfo(name="Seoul", region="Korea", timezone="Asia/Seoul", latitude=observer_lat, longitude=observer_lon)
tm = datetime.now(tz=utc) + timedelta(days=days_after)

print('현재 시각:     ', tm.astimezone(timezone(location.timezone)))

moon_rise, moon_set = get_moon_rise_moon_set(location, tm)
print('달 뜨는 시각:  ', moon_rise.astimezone(timezone(location.timezone)))
print('달 지는 시각:  ', moon_set.astimezone(timezone(location.timezone)))

duration = int((moon_set - moon_rise).total_seconds())
print('달 떠있는 시간:', duration, '초')

print('---------------------------------')

다음으로 주어진 시간에 대해 달의 방위각, 천정각을 계산할 차례.
참고로, 천정각은 90도에서 고도를 뺀 값의 개념이다.

1
2
3
4
5
6
7
8
9
10
11
12

# 주어진 시간에 대해 달의 방위각, 천정각을 구한다.
moon_times = []

for i in np.arange(0, duration, step_in_second):
    moon_times.append(moon_rise + timedelta(seconds=i))

moon_az = np.zeros(len(moon_times))
moon_ze = np.zeros(len(moon_times))

for i, mt in enumerate(moon_times):
    moon_az[i] = moon.azimuth(location.observer, mt)
    moon_ze[i] = moon.zenith(location.observer, mt)

그 다음은 ISS의 방위각, 천정각을 계산.
ISS의 XYZ 좌표를 구하고, 이를 관측자의 위치를 기준으로 방위각, 천정각으로 변환한다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 주어진 시간에 대해 ISS의 XYZ 좌표를 구한다.
ts = load.timescale()
moon_times = ts.from_datetimes(moon_times)

iss = get_iss_object_skyfield()
iss_geocentric = iss.at(moon_times)
iss_xyz = iss_geocentric.itrf_xyz().to(u.m)
iss_xyz = iss_xyz.transpose()

# ISS의 XYZ 좌표를 이용하여 azimuth, zenith를 구한다.
iss_az_ze = calc_azimuth_zenith((observer_lat, observer_lon, observer_alt), iss_xyz)

iss_az = iss_az_ze[:, 0]
iss_ze = iss_az_ze[:, 1]

마지막으로 두 각들의 사잇각을 계산하면 된다.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# 달과 ISS의 사잇각을 구한다.
angles = calc_angle_between(moon_az, moon_ze, iss_az, iss_ze)

# moon_times 값들을 local 시간으로 변환한다.
moon_times = [mt.astimezone(timezone(location.timezone)) for mt in moon_times]

# angles가 최솟값인 시각을 찾는다.
min_angle = np.min(angles)
min_index = np.argmin(angles)
min_time = moon_times[min_index]

print('달과 ISS 사이가 가장 가까운 각도:', min_angle)
print('달과 ISS 사이가 가장 가까운 시각:', min_time)

my_xyz = wgs84_llh_to_xyz((observer_lat, observer_lon, observer_alt))
min_iss_pos = iss_xyz[min_index].to(u.m).value
min_iss_xyz = min_iss_pos - my_xyz
min_distance = np.sqrt(min_iss_xyz[0] ** 2 + min_iss_xyz[1] ** 2 + min_iss_xyz[2] ** 2)
print('그 때 ISS 까지의 거리:          ', min_distance / 1000, 'km')
print('그 때 ISS의 방위각:             ', iss_az[min_index], '도, 천정각:', iss_ze[min_index], '도')


## 그래프 그리기
plt.plot(moon_times, angles)
plt.xlabel('Date/Time')
plt.ylabel('Angle between Moon and ISS')
plt.title('Angle between Moon and ISS')

plt.show()

계산 결과

실행 결과는 다음과 같다.
달과 ISS 사잇각이 가장 작은 시각은 4월 8일 09시 42분 32.5초로, 이 때의 각도는 4.69도이다.
백주대낮이니 달이 제대로 보이진 않을 것이고…

현재 시각:      2024-04-07 18:23:46.015400+09:00
달 뜨는 시각:   2024-04-08 05:45:00+09:00
달 지는 시각:   2024-04-08 18:30:00+09:00
달 떠있는 시간: 45900 초
---------------------------------
[#################################] 100% iss.txt
달과 ISS 사이가 가장 가까운 각도: 4.694649913207852
달과 ISS 사이가 가장 가까운 시각: 2024-04-08 09:42:32.500000+09:00
그 때 ISS 까지의 거리:           543.9669704560661 km
그 때 ISS의 방위각:              129.85901531930875 도, 천정각: 41.67761609274922 도

시간에 따른 사잇각의 그래프는 아래와 같다.

image 04-08 10 바로 왼쪽의 최저점이 위에 적은 4.69도

관측이 가능한 시간이 나오더라도 유지되는 시간은 1초 미만으로 극히 짧다.
즉, 이렇게 찍힌 사진들은 그 자체만으로 기적과도 같은 것이다.

카테고리:

업데이트:

참고