2018년 12월 5일에 발사된 천리안위성 2A호는 천리안위성 1호의 기상관측 역할을 승계하는 차세대 정지궤도 기상위성으로 기상 및 우주기상 관측임무를 수행합니다. 천리안위성 2A호는 16개의 기본 채널 관측을 통해 52종의 기상산출물을 생산하며, 국내에서 개발한 우주기상 탑재체로 우주기상 감시 및 위성 운영과 관련되는 현상을 관측합니다. 이로써 우리나라의 기상관측, 수치예보, 기후감시, 우주 기상감시의 역량을 강화했습니다.

전구(FD), 동아시아(EA), 확장 국지 영역(ELA), 한반도(KO), 태풍 감시 영역(TP), 지역관측(LA)을 제공하고 있습니다.

영역별_도법_및_영상정보.pdf

가시_0.47_파랑

가시_0.51_초록

가시_0.64_빨강

가시_0.86_식생

근적외_1.37_권운

근적외_1.6_눈얼음

단파적외_3.8_야간안개하층운

수증기_6.3_상층수증기

수증기_6.9_중층수증기

수증기_7.3_하층수증기

적외_8.7_구름상

적외_9.6_오존

적외_10.5_깨끗한대기창

적외_11.2_대기창

적외_12.3_오염된대기창

적외_13.3_이산화탄소

• 사이트 : 기상자료개방포털 (data.kma.go.kr)​
• 메뉴 : 데이터 → 기상위성 → 천리안위성 2A호

• 조건 : 자료형태, 자료구분, 기간, 요소

• 사이트 : 기상청 API허브 (apihub.kma.go.kr)
• 메뉴 : 위성 → 천리안 위성 2A호

• 절차: 회원가입 및 로그인 → API 메뉴 이동 → API 활용신청 → API호출
• 호출 방법: 기상청_api허브_사용_안내서.pdf

소스코드: gk2a_l1b.ipynb

# API 요청을 위한 requests 라이브러리와
# netCDF4 포맷 파일을 읽기 위한 netCDF4 라이브러리,
# 분포도 표출을 위한 matplotlib 라이브러리와
# 한반도 지형 표출을 위한 geopandas 라이브러리를 설치합니다.
%pip install requests netCDF4 matplotlib geopandas

기상청 API허브에 회원가입 후 자신의 인증키를 사용하여 API 요청을 할 수 있습니다.

API 인증키를 코드상에 직접 남겨 사용하는 방식은 다른 사람에게 코드를 공유할때 인증키가 노출되는 문제가 생길 수 있습니다.

따라서 인증키를 환경변수로 등록하여 코드상에 노출되지 않도록 사용하는 것이 바람직합니다.

import os # 환경변수를 불러오는 표준 라이브러리

# 바람직하지 않은 방법
# my_api_key = 'MY_PRIVATE_API_AUTHENTICATE_KEY'

# 인증키를 'myapikey'라는 이름의 환경변수로 미리 등록해놓은 뒤 이를 불러와서 사용합니다.
my_api_key = os.environ['myapikey']

import requests # API 요청을 보내고 받는 라이브러리

# 시간을 다루는 파이썬 기본 패키지
from datetime import datetime, timezone, timedelta

# 천리안위성 2A호 파일 다운로드 API url
api_url = 'https://apihub.kma.go.kr/api/typ01/url/sat_file_down2.php'

# 요청인자로 자료 형태, 관측 시간, 레벨, 채널번호 및 산출물 종류, 영역, 도법이 있습니다.

# 1. 자료 형태를 입력합니다.
# 이진파일(bin)과 이미지(img)를 입력할 수 있습니다.
data_type = 'bin'

# 2. 관측 시간을 년월일시분(UTC)으로 입력합니다.
# 자료의 기준시간이 UTC이므로 대한민국 시간(KST)의 9시간 전이 UTC 시간이 됩니다.
# 여기선 KST기준으로 값을 입력하고 UTC로 변환하는 과정을 거칩니다.
# 이후 입력할 영역에 따라 전 지구영역의 경우 10분 간격으로,
# 아시아/태평양 지역과 한반도의 경우 2분 간격으로 입력할 수 있습니다.
# 이때, 가시영역의 채널의 경우 낮 시간대로 설정해야 영상을 확인할 수 있습니다.
kst_time = '202404011200'

# 입력한 대한민국 시간으로부터 9시간 전인 UTC 시간을 계산합니다.
# 문자로 입력한 시간을 datetime으로 변환하고 시간대를 변환한 뒤, 다시 문자열로 변환합니다.
utc_time = datetime.strptime(
    kst_time, '%Y%m%d%H%M'
).replace(
    tzinfo=timezone(timedelta(hours=9))
).astimezone(
    timezone.utc
).strftime(
    '%Y%m%d%H%M'
)

# 3. 자료의 레벨을 입력합니다.
# 레벨의 경우 기본관측자료(L1B), 기상산출물(L2)을 입력할 수 있습니다.
level = 'L1B'

# 4. 레벨에 따른 채널번호 및 산출물 종류를 입력합니다.
# 레벨 L1B의 경우, 16개의 관측 채널 중 하나를 입력할 수 있습니다.
# 가시영역의 채널의 경우, VI004(파랑), VI005(초록), VI006(빨강), VI008(식생)이 존재합니다.
# 파랑의 경우 지표 및 해안 에어로솔,
# 초록은 안개, 연무 탐지 및 트루컬러 합성,
# 빨강은 구름, 안개, 일사량,
# 식생은 구름, 안개, 일사량 관측에 사용되며
# 빨강, 초록, 파랑을 합쳐 트루컬러로 표출할 수 있습니다.
channel = 'VI004'

# 5. 자료가 나타내는 영역을 입력합니다.
# 영역의 경우 전구(FD), 동아시아(EA, ELA), 한반도(KO), 지역관측(LA), 태풍(TP)을 선택할 수 있습니다.
area = 'KO'

# 6. 자료에 투영된 도법을 입력합니다.
# 도법의 경우 한반도 영역의 경우 Lambert Conic(LC) 도법을 사용합니다.
projection = 'LC'

# API 요청인자들을 묶어 dictionary로 정의합니다.
# 이때, 각 인자는 소문자로 입력하는 것에 유의합니다.
api_parameters = {
    'typ': data_type,
    'tm': utc_time,
    'lvl': level.lower(),
    'chn': channel.lower(),
    'are': area.lower(),
    'map': projection.lower(),
    'authKey': my_api_key
}

파일의 명명규칙과 설명은 다음과 같습니다.

gk2a_ami_le1b_BAND_scnRESpr_yyyymmddhhmn.nc

• gk2a - 천리안위성 2A호(GEO-KOMPSAT-2A:GEOstationary KOrea Multi-Purpose SATellite-2A)

• ami - 고성능 기상센서(Advanced Meteorological Imager)

• le1b - 원시자료(Raw data)에서 분해(LV0), 복사보정(LV1A), 기하보정(LV1B)의 절차를 거친 전처리된 자료(이후 기상요소를 산출하는 LV2 과정)

• BAND - 기본관측자료별 채널과 파장( 예 - VI004는 16개 관측 채널중 첫 번째 채널인 가시영역(파랑) 채널로 파장대는 0.470㎛를 사용합니다.)

• scn - 관측 영역으로 전구 영역(Full Disk), 확장 국지 영역(Extended Local Area), 국지 영역(Local Area), 한반도 영역(KOrea) 등이 있습니다.

• RES - 공간 해상도로 100m를 단위로 세 자리로 표현합니다. (0.5km → 005, 1km → 010, 2km → 020)

• pr - 도법으로 GEOS(ge), Lambert Conformal Conic(lc) 두 가지가 존재합니다.

• yyyymmddhhmn - 관측 시간으로 UTC 시간을 기준으로 합니다.

# 각 레벨별, 채널별 파일 명명법을 정의한 dictionary
params_to_filename = {
    # 레벨별 명명법
    'L1B': 'le1b',
    'L2': 'le2',

    # 채널별 해상도
    'VI004': '010',
    'VI005': '010',
    'VI006': '005',
    'VI008': '010',
    'NR013': '020',
    'NR016': '020',
    'SW038': '020',
    'WV063': '020',
    'WV069': '020',
    'WV073': '020',
    'IR087': '020',
    'IR096': '020',
    'IR105': '020',
    'IR112': '020',
    'IR123': '020',
    'IR133': '020',
}

# 파일 이름은 명명규칙에 따라 정해지며 API 응답을 저장할 파일 형식은 netCDF4파일(nc)로 저장합니다.
file_name = (
    f'gk2a_ami_{params_to_filename[level]}_{channel.lower()}_'
    f'{area.lower()}{params_to_filename[channel]}{projection.lower()}_{utc_time}.nc'
)

# 저장된 파일이 없는 경우에만 API를 요청해 파일을 다운로드합니다.
if not os.path.exists(file_name):
    # API 요청인자와 함께 API 요청
    response = requests.get(api_url, params=api_parameters)

    # 잘못된 응답을 받거나 짧은 에러메세지를 응답으로 받은 경우 에러 메세지 출력
    if response.status_code != 200 or len(response.content) < 500:
        raise requests.RequestException(
            response.content.decode('euc-kr')
        )
    else:
        # 그 외의 올바른 응답에 대해서만 파일로 저장합니다.
        with open(file_name, 'wb') as f:
            f.write(response.content)
        print(f'{file_name} 파일 다운로드 완료')

가시영상을 트루컬러로 표출하기 위해 RGB 3개 채널의 파일을 다운로드 합니다. (VI004, VI005, VI006)

다운로드 받은 이진파일은 netCDF4 포맷이며 netCDF4 라이브러리를 통해 파일을 읽을 수 있습니다.

netCDF 포맷의 가장 큰 특징으로 파일 내에 데이터뿐만 아니라 데이터를 설명하는 관련 정보들도 함께 담겨있습니다.

import netCDF4  # netCDF4 파일을 읽기 위한 라이브러리

# 파일을 읽기모드로 열어 파일내의 관련 정보들을 출력합다.
with netCDF4.Dataset(file_name, mode='r', format='netcdf4') as ncfile:
    print(ncfile)

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    projection_type: lambert_conformal_conic
    standard_parallel1: 30.0
    standard_parallel2: 60.0
    ...
    satellite_name: GK-2A
    dimensions(sizes): dim_y(1800), dim_x(1800)
    variables(dimensions): uint16 image_pixel_values(dim_y, dim_x)
    groups: 

관련 정보들 중 projection_type, standard_parallel1, standard_parallel2, origin_latitude, central_meridian, false_easting, false_northing은 Lambert Conic 도법에 관련된 정보인 것을 알 수 있습니다.

또, channel_spatial_resolution이 공간해상도를 의미하는 것도 알 수 있습니다.

실제 데이터(variables)는 image_pixel_values에 uint16형으로 저장되어 있으며 데이터의 차원(dimensions)은 dim_y(1800), dim_x(1800)로 정의되어 있습니다.

# 파일을 다시 읽기모드로 열어 파일 내의 관련 정보들을 변수에 저장합니다.
with netCDF4.Dataset(file_name) as ncfile:
    # Lambert Conic 도법에 관련된 정보를 저장합니다.
    lat_1 = ncfile.standard_parallel1
    lat_2 = ncfile.standard_parallel2
    lat_0 = ncfile.origin_latitude
    lon_0 = ncfile.central_meridian
    false_easting = ncfile.false_easting
    false_northing = ncfile.false_northing

    # 데이터의 공간해상도 정보를 저장합니다.
    resolution = ncfile.pixel_size

    # 격자 자료의 모서리 좌표 정보를 저장합니다.
    lower_bottom_x = ncfile.lower_left_easting / resolution
    lower_bottom_y = ncfile.lower_left_northing / resolution
    upper_top_x = ncfile.upper_right_easting / resolution
    upper_top_y = ncfile.upper_right_northing / resolution

    # 실제데이터(variables)인 image_pixel_values 정보들을 출력합니다.
    print(ncfile['image_pixel_values'])

<class 'netCDF4._netCDF4.Variable'>
uint16 image_pixel_values(dim_y, dim_x)
    number_of_total_pixels: 3240000
    average_pixel_value: 515.8856228070175
    channel_name: VI004
    max_pixel_value: 2044
    min_pixel_value: 178
    number_of_total_bits_per_pixel: 16
    number_of_valid_bits_per_pixel: 11
unlimited dimensions: 
current shape = (1800, 1800)
filling on, default _FillValue of 65535 used

관련 정보들 중 값이 가질 수 있는 최대값(max_pixel_value), 최소값(min_pixel_value)과 자료의 유효비트수(number_of_valid_bits_per_pixel)를 확인할 수 있습니다.

uint16 형의 16개 비트 수 중 유효비트가 아닌 앞의 2자리 비트는 데이터 품질 표식(Quality Flag)으로써

00은 좋은 품질을, 01은 조건에 따라 사용가능(conditionally usable),

10은 관측 영역 외부를, 11은 오류 존재를 의미합니다.

여기선 품질 표식 값에 상관없이 유효비트가 나타내는 값을 사용해 가시영상을 표출합니다.

각 채널의 값을 트루컬러로 합성하기 전, 각 값들은 0에서 1사이의 float 값들로 변환되어야 합니.

국가기상위성센터(NMSC)에서 보정표(Calibration Table)를 이용해 각 채널의 값을 복사휘도(Radiance)와 알베도(Albedo)로 변환하여 사용할 수 있습니다.

보정표에서 명시한 각 가시영역 채널별 복사휘도, 알베도 변환 계수는 다음과 같습니다.

여기에선 각 채널의 값들을 복사휘도를 거쳐 알베도 값으로 변환하여 사용합니다.

import numpy as np  # 배열을 다루는 라이브러리

# 파랑 채널의 배열을 파일로부터 읽어들입니다.
with netCDF4.Dataset(f'gk2a_ami_le1b_vi004_ko010lc_{utc_time}.nc') as ncfile:
    valid_bits = ncfile['image_pixel_values'].number_of_valid_bits_per_pixel
    arr_b = np.asarray(ncfile['image_pixel_values'])
# 유효비트수만 사용해 uint16에서 float64로 형변환하고, 복사 휘도를 거쳐 알베도로 변환합니다.
arr_b = np.bitwise_and(arr_b, (1 << valid_bits) - 1).astype(np.float64)
arr_b = ((0.363545805215835 * arr_b - 7.27090454101562) * 0.0015582450)

# 초록 채널의 배열을 파일로부터 읽어들입니다.
with netCDF4.Dataset(f'gk2a_ami_le1b_vi005_ko010lc_{utc_time}.nc') as ncfile:
    valid_bits = ncfile['image_pixel_values'].number_of_valid_bits_per_pixel
    arr_g = np.asarray(ncfile['image_pixel_values'])
# 유효비트수만 사용해 uint16에서 float64로 형변환하고, 복사 휘도를 거쳐 알베도로 변환합니다.
arr_g = np.bitwise_and(arr_g, (1 << valid_bits) - 1).astype(np.float64)
arr_g = ((0.343625485897064 * arr_g - 6.87249755859375) * 0.0016595767)

# 빨강 채널의 배열을 파일로부터 읽어들입니다.
# 빨강 채널은 다른 채널에 비해 공간해상도가 2배 세밀하므로
# 배열을 슬라이싱하여 다른 배열과의 크기를 맞추어 사용합니다.
with netCDF4.Dataset(f'gk2a_ami_le1b_vi006_ko005lc_{utc_time}.nc') as ncfile:
    valid_bits = ncfile['image_pixel_values'].number_of_valid_bits_per_pixel
    # 배열의 x, y 좌표를 두 칸씩 건너뛰어 슬라이싱합니다.
    arr_r = np.asarray(ncfile['image_pixel_values'])[::2, ::2]
# 유효비트수만 사용해 uint16에서 float64로 형변환하고, 복사 휘도를 거쳐 알베도로 변환합니다.
arr_r = np.bitwise_and(arr_r, (1 << valid_bits) - 1).astype(np.float64)
arr_r = ((0.154856294393539 * arr_r - 6.19424438476562) * 0.0019244840)

마지막으로 위성 가시영상에서 지형의 경계를 명확히 표시하기 위해 국토지리정보원에서 제공하는 대한민국주변도 shape 파일을 다운로드 받아 압축을 해제합니다.

import requests # API 요청을 보내고 받는 라이브러리
import zipfile  # zip 압축 파일을 다루는 모듈

# 국토지리정보원의 대한민국주변도 shape 파일 다운로드 url
shp_url = 'https://www.ngii.go.kr/other/file_down.do?sq=107951'

# shape 파일은 압축파일 형태로 제공됩니다.
shape_file_name = 'korea_shp.zip'
extract_directory = 'korea_shape'

# shape 압축파일이 없는 경우에만 API를 요청해 파일을 다운로드합니다.
if not os.path.exists(shape_file_name):
    # url에 대해 GET 요청을 보냅니다.
    response = requests.get(shp_url)

    # 받은 응답을 zip 파일 형태로 저장합니다.
    with open(shape_file_name, 'wb') as f:
        f.write(response.content)
    print(f'{shape_file_name} 파일 다운로드 완료')

    # 다운로드 받은 압축 파일을 korea_shape 폴더 아래에 압축 해제합니다.
    with zipfile.ZipFile(shape_file_name, 'r') as zip_ref:
        for file_info in zip_ref.infolist():
            file_info.filename = file_info.filename.encode('cp437').decode('euc-kr')
            zip_ref.extract(file_info, extract_directory)

import geopandas as gpd # shape 파일을 다루는 라이브러리

# 국가 경계를 구분한 shp파일은 ARD_NAION_AS으로
# 이 파일을 읽어 LCC 도법으로 변환합니다.
gdf = gpd.read_file(
    os.path.join(extract_directory, 'ARD_NAION_AS.shp')
).to_crs(
    f"+proj=lcc +lat_1={lat_1} +lat_2={lat_2} +lat_0={lat_0} +lon_0={lon_0} "
    f"+x_0={false_easting} +y_0={false_northing} +ellps=WGS84 +units=m +no_defs"
)

# 1m 단위의 도법에서 위성 해상도 단위(1000m)에 맞춥니다.
gdf.geometry = gdf.geometry.scale(
    xfact=1/resolution, yfact=1/resolution, zfact=1.0, origin=(0, 0)
)

import matplotlib.pyplot as plt # 이미지 표출에 사용되는 라이브러리

# (1800 x 1800) 크기의 단일 채널의 이미지를 rgb 순서로 합쳐
# (1800 x 1800 x 3)의 트루컬러로 합성합니다.
# 0-1 사이의 float64 값을 0-255 사이의 uint8로 형변환합니다.
img = (np.dstack([arr_r, arr_g, arr_b]) * 255).astype(np.uint8)

# 이미지를 그릴 크기를 지정합니다.
fig, ax = plt.subplots(1, 1, figsize=(10, 10))

# 이미지의 제목을 지정합니다.
ax.set_title(f"GK2A RGB {datetime.strptime(kst_time, '%Y%m%d%H%M')}(KST)")

# 국가간의 경계선을 노란색으로 그립니다.
gdf.boundary.plot(ax=ax, color='y', linewidth=0.5)

# 이미지를 해당 크기에 맞춰 그립니다.
ax.imshow(img, extent=[lower_bottom_x, upper_top_x, lower_bottom_y, upper_top_y])

# 그린 이미지를 표출합니다.
plt.show()