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

전구, 동아시아, 한반도

존데 자료 비교(2019. 8. 1. ~ 8. 31. 1개월)

MODIS/AQUA, TERRA 위성의 결과와 비교(2020. 4. 1. ~ 4. 30.)

레이더 CAPPI 1.5 km 및 낙뢰 자료와 비교(2019. 8. 1. ~ 8. 20., 총 7개 대류운 사례)

MODIS/AQUA, TERRA 위성의 결과와 비교(2020. 4. 1. ~ 4. 30.)

MODIS 위성의 운형결과와 비교(2019.8.1. ~ 8.31.)

MODIS 센서의운상결과와비교(2020.1.1.~ 1.31.)

MODIS 센서의운상결과와비교(2020.1.1.~ 1.31.)

(COT) : 주간→ bias 20%, 야간 → bias 30 %

(CER) : 수적→ bias 4㎛, 빙정 → bias 10 ㎛

(LWP) : bias 25 or 15 %

(IWP) : bias 25 or 30 %

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

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

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

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

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

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

# API 요청을 위한 requests 라이브러리와
# netCDF4 포맷 파일을 읽기 위한 netCDF4 라이브러리,
# 분포도 표출을 위한 matplotlib 라이브러리와
# 이미지 및 점의 좌표계 변환을 위한 rasterio, pyproj 라이브러리,
# 지도를 표출하기 위한 folium 라이브러리를 설치합니다.
%pip install requests netCDF4 matplotlib rasterio pyproj folium

기상청 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분 간격으로 생산됩니다.
kst_time = '202404171200'
 
# 입력한 대한민국 시간으로부터 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 = 'L2'
 
# 4. 레벨에 따른 채널번호 및 산출물 종류를 입력합니다.
# 레벨 L2의 경우, 기상산출물 중 하나를 입력합니다.
# 여기선 기상산출물 중 예시로 안개(FOG)를 입력합니다.
data = 'FOG'
 
# 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(),
    'dat': data.lower(),
    'are': area.lower(),
    'map': projection.lower(),
    'authKey': my_api_key
}

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

gk2a_ami_le2_DATA_scnRESpr_yyyymmddhhmn.nc

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

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

• le2 - 기본채널 관측자료를 통해 산출된 기상산출물 레벨2

• DATA - 기상산출물 종류( 예 - 구름탐지(CLD), 안개(FOG), 강우강도(RR) 등 )

• 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',
}
 
# 파일 이름은 명명규칙에 따라 정해지며 API 응답을 저장할 파일 형식은 netCDF4파일(nc)로 저장합니다.
file_name = (
    f'gk2a_ami_{params_to_filename[level]}_{data.lower()}_'
    f'{area.lower()}020{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} 파일 다운로드 완료')

gk2a_ami_le2_fog_ko020lc_202404170300.nc 파일 다운로드 완료

다운로드 받은 이진파일은 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):
  file_creation_time: 2024-04-17T03:19:08.871Z
  origianl_sourece_file: gk2a_ami_le2_fog_ea020ge_202404170300.nc
  algorithm_version: GK2A_FOG_v1.2023.01.1
  comment: TBD
  conventions: CF-1.7
  institution: KMA/NMSC> Korea Meteorological Administration, National Meteorological Satellite Center
  instrument: GK-2A Advanced Meteorological Imager
  license: Access is restricted to approved users only
  metadata_conventions: Unidata Dataset Discovery v1.0
  names_of_product: FOG, DQF_FOG
  number_of_product: 2
  observation_mode: FD
  processing_area: EA
  processing_environment: operation
  references: RAP03 FOG v1.0
  title: AMI L2 FOG
  dimensions(sizes): dim_y(900), dim_x(900)
  variables(dimensions): uint8 DQF_FOG(dim_y, dim_x), int16 Del_Fta(dim_y, dim_x), uint16 FOG(dim_y, dim_x), int32 gk2a_imager_projection()
  groups:

실제 데이터는 netCDF의 variables에 저장되어있으며 순서대로 DQF_FOG, Del_Fta, FOG, gk2a_imager_projection이 존재합니다.

각 데이터를 살펴보면 다음과 같습다.

# 파일을 읽기모드로 열어 파일내의 관련 정보들을 출력합니다.
with netCDF4.Dataset(file_name) as ncfile:
    # 실제데이터 FOG 정보들을 출력합니다.
    # 이외에도 DQF_FOG, Del_Fta, gk2a_imager_projection을 살펴볼 수 있습니다.
    var_name = 'FOG'
    print(ncfile[var_name])

<class 'netCDF4._netCDF4.Variable'>
uint16 FOG(dim_y, dim_x)
  _FillValue: 65535
  ancillary_variables: DQF_FOG
  grid_mapping: gk2a_imager_projection
  long_name: AMI L2 FOG 2km
  product_meanings: 1 : Clear, 2 : Middle or High Cloud, 3 : Unkown, 4 : Probably Fog, 5 : Fog, 6 : Snow, 7 : Desert or semi-Desert
  units: 
  valid_max: 7
  valid_min: 1
unlimited dimensions: 
current shape = (900, 900)
filling on

살펴본 결과는 다음과 같습니다.

FOG는 각 픽셀별로 안개 여부(1 ~ 7)를 저장한 변수이고,

DQF_FOG는 FOG 변수에 대한 데이터 품질 표식(Data Quality Flag)(0 ~ 15)을 나타냅니다.

Del_Fta는 안개로 판명된 지점에 대해 이미지 표출을 위한 변수이며

gk2a_imager_projection는 데이터의 도법에 관한 설명을 담고있습니다.

따라서 Del_Fta 변수를 사용해 안개이미지를 표출합다.

import matplotlib.pyplot as plt # 분포도 표출을 위한 라이브러리
import numpy as np              # 배열을 다루는 라이브러리
 
# 분포도 표출을 위한 컬러맵
from matplotlib.colors import ListedColormap, BoundaryNorm, Normalize
 
# 안개 배열을 분포도로 표출하기 위해 6가지 색의 색상표를 RGB로 정의합니다.
colormap_fog = ListedColormap(np.array([
    [255, 0, 0], [192, 0, 0], [255, 192, 0],
    [255, 255, 0], [0, 176, 80], [146, 208, 80]
]) / 255)
 
# 색상표에서 벗어난 값은 투명한 색상(RGBA 중 A가 0)을 나타내도록 합니다.
colormap_fog.set_bad([0, 0, 0, 0])
 
# 색상표에서 각각의 색이 나타내는 범위를 정합니다.
bounds = np.array([
    0, 1, 2, 3, 4, 5, 6
])
 
# 색상 개수와 값 범위를 서로 맞춥니다.
norm = BoundaryNorm(boundaries=bounds, ncolors=len(colormap_fog.colors))
 
# 색상표를 표출합니다.
colormap_fog

# 파일을 다시 읽기모드로 열어 파일 내의 관련 정보들을 변수에 저장합니다.
with netCDF4.Dataset(file_name) as ncfile:
    # 이미지 표출을 위해 배열과 배열의 비어있는 영역을 나타내는 마스크를 정의합니다.
    arr_fog = np.asarray(ncfile['Del_Fta'])
    null_mask = (arr_fog <= -30000)
 
    # 데이터의 단위를 저장합니다.
    unit = ncfile['Del_Fta'].units
 
    # Lambert Conic 도법에 관련된 정보를 변수에 저장합니다.
    projection_var = ncfile['gk2a_imager_projection']
 
    lat_1 = projection_var.standard_parallel1
    lat_2 = projection_var.standard_parallel2
    lat_0 = projection_var.origin_latitude
    lon_0 = projection_var.central_meridian
    false_easting = projection_var.false_easting
    false_northing = projection_var.false_northing
    source_width = projection_var.image_width
    source_height = projection_var.image_height
    source_resolution = projection_var.pixel_size
 
    lower_bottom_x = projection_var.lower_left_easting
    lower_bottom_y = projection_var.lower_left_northing
    upper_top_x = projection_var.upper_right_easting
    upper_top_y = projection_var.upper_right_northing

# 분포도에서 색상표 위치, 크기를 조절하는 함수
from mpl_toolkits.axes_grid1 import make_axes_locatable
 
# 안개 배열 각각의 값을 색상표의 색이 나타내는 인덱스로 매핑시킵니다.
# 매핑된 각 인덱스 범위(0 ~ 6)를 정규화(0 ~ 1)합니다.
# 정규화를 통해 잃게된 투명색의 범위를 다시 지정합니다.
# 마지막으로 배열의 형태를 (행 개수, 열 개수) 에서 (행 개수, 열 개수, RGBA)로 변환합니다.
colored_array = norm(arr_fog)
colored_array = Normalize(
    0, len(colormap_fog.colors)
)(colored_array)
colored_array[null_mask] = np.nan
colored_array = (colormap_fog(colored_array) * 255).astype(np.uint8)
 
# 분포도를 그릴 크기를 지정합니다.
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
 
# 분포도의 제목을 지정합니다.
ax.set_title('Colored array')
 
# 배경색을 회색으로 지정합니다.
ax.set_facecolor('#cccccc')
 
# 색상표의 크기와 위치를 조절합니다.
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0)
 
# 배열을 해당 크기에 맞춰 그립니다.
im = ax.imshow(colored_array, cmap=colormap_fog, norm=norm)
 
# 색상표에 표시될 글자 크기 및 제목을 설정합니다.
cbar = fig.colorbar(im, cax=cax)
cbar.ax.tick_params(labelsize=8)
cbar.ax.set_title(unit, fontsize=8)
 
 
# 그린 이미지를 표출합니다.
plt.show()

위와 같은 이미지를 동적인 지도 위에 표출하기 전, 지도와 위성 자료(이미지)의 좌표계(Coordinate Reference System, CRS)를 맞추어야 합니다.

현재 동적인 지도가 사용하는 좌표계는 웹 상에서 지도를 표출할 때 많이 사용되는 EPSG:3857입니다.

위성 자료가 사용하는 좌표계는 데이터의 gk2a_imager_projection 변수를 통해 확인 가능합니다.

위성 자료가 사용하는 좌표계는 Lambert Conformal Conic Projection(LCC)으로 한반도를 중심으로 표현하기 적합해 여러 데이터에서도 사용된다고 합니다.

위성 자료를 EPSG:3857로 변환하여 지도 위에 표출합니다.

# 좌표계 변환을 위한 변환 행렬
from rasterio.transform import Affine
 
# 좌표계 변환을 위해 필요한 함수
from rasterio.warp import calculate_default_transform, reproject, Resampling
 
# 위성 자료의 좌표계를 proj.4 형태의 문자열로 정의합니다.
# 투영법(Projection)은 LCC이며 LCC 좌표계 정의에 필요한 위도 선 2개(lat_1, lat_2)와
# 중심 위도, 경도(lat_0, lon_0)을 정의합니다.
# 좌표계의 x, y는 각각 오른쪽, 위로 증가하는 방향을 가지고(False easting, False northing)
# 지구 타원체를 WGS84로 정의합니다.
# 마지막으로 좌표계가 사용하는 단위는 미터(m)로 정의합니다.
source_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"
)
 
# 이미지(배열)의 행과 열이 LCC 좌표계에서의 좌표로 변환되기 위한 변환 행렬(Affine Transform Matrix)을 정의합니다.
# 먼저 각 이미지 픽셀이 실제로 2000m를 나타내므로
# 변환 행렬은 크기를 2000배 키워주는 동시에 x축 대칭을 하고, (Scale, Reflection)
# 이미지의 중심점을 (0, 0) 위치로 이동시키는(Translation) 역할을 합니다.
# 이때, 행렬의 곱하는 순서에 유의합니다.(AB != BA)
source_transform = Affine.translation(-upper_top_x, upper_top_y) * Affine.scale(source_resolution, -source_resolution)
 
# 변환 행렬을 거친 이미지가 나타내는 경계를 정의합니다.
source_bounds = {
    'left': lower_bottom_x,
    'bottom': lower_bottom_y,
    'right': upper_top_x,
    'top': upper_top_y
}
 
# 변환 후 이미지의 변환 행렬과 너비와 높이를 계산합니다.
dest_transform, dest_width, dest_height = calculate_default_transform(
    src_crs=source_crs,
    dst_crs='EPSG:3857',
    width=source_width,
    height=source_height,
    **source_bounds,
)
 
# 변환 후의 이미지가 담길 비어있는 배열을 정의합니다.
converted_array = np.ones((dest_height, dest_width, 4), dtype=np.uint8)
 
# RGBA 각 채널에 대해 좌표계 변환을 수행합니다.
# 사용하는 resampling 기법으로 가까운 값을 선택하는 nearest를 선택합니다.
for i in range(4):
    reproject(
        source=colored_array[:, :, i],
        destination=converted_array[:, :, i],
        src_transform=source_transform,
        src_crs=source_crs,
        dst_transform=dest_transform,
        dst_crs='EPSG:3857',
        resampling=Resampling.nearest,
    )
 
# 변환된 이미지를 분포도로 표출합니다.
# 이전 분포도와 다르게 너비, 높이가 달라졌음을 확인합니다.
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
 
# 분포도의 제목을 지정합니다.
ax.set_title('Converted array')
 
# 배경색을 회색으로 지정합니다.
ax.set_facecolor('#cccccc')
 
# 색상표의 크기와 위치를 조절합니다.
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0)
 
# 배열을 해당 크기에 맞춰 그립니다.
im = ax.imshow(converted_array, cmap=colormap_fog, norm=norm)
 
# 색상표에 표시될 글자 크기 및 제목을 설정합니다.
cbar = fig.colorbar(im, cax=cax)
cbar.ax.tick_params(labelsize=8)
cbar.ax.set_title(unit, fontsize=8)
 
# 그린 이미지를 표출합니다.
plt.show()

마지막으로 변환된 이미지를 동적인 지도 위에 표출합니다.

import folium   # 지도로 표출하기 위한 라이브러리
import branca.colormap as cm    # 지도에 색상표를 표출하기 위한 색상표 모듈
 
# 지도에 마우스 위치를 표출하는 플러그인
from folium.plugins import MousePosition
 
# 점끼리의 좌표계 변환을 위한 transformer
from pyproj.transformer import Transformer
 
# folium을 이용한 좌표는 경도, 위도만 입력 가능하므로
# EPSG:3857 좌표계에서 EPSG:4326 좌표계로의 변환하는 transformer를 정의합니다.
degree_tranformer = Transformer.from_crs('EPSG:3857', 'EPSG:4326')
 
# 너비, 높이 크기가 800인 지도 생성
fig = folium.Figure(width=800, height=800)
 
# 지도 중심의 경도, 위도와 지도가 나타낼 경도, 위도의 범위를 정한 맵을 생성합니다.
map = folium.Map(
    location=[38, 126],
    zoom_start=6,
    min_zoom=5,
    min_lat=29,
    max_lat=46,
    min_lon=113,
    max_lon=139,
    max_bounds=True,
    tiles='Cartodb Positron',
).add_to(fig)
 
# 지도에 색상표를 추가합니다.
map.add_child(cm.StepColormap(
    [tuple(i) for i in colormap_fog.colors],
    vmin=bounds[0], vmax=bounds[-1], tick_labels=bounds.tolist(), caption=unit
))
 
# 변환된 이미지를 지도 위에 겹쳐 그립니다.
# 불투명도를 0.4로 정하고
# 이미지가 그려질 범위를 경도, 위도로 입력합니다.
folium.raster_layers.ImageOverlay(
    image=converted_array,
    name='fog',
    opacity=0.4,
    bounds=[
        degree_tranformer.transform(*dest_transform.__mul__((0, dest_height))),
        degree_tranformer.transform(*dest_transform.__mul__((dest_width, 0)))
    ],
    zindex=1,
).add_to(map)
 
# 지도에 마우스 위치를 표출하는 플러그인과 layer를 조절하는 버튼을 추가합니다.
MousePosition().add_to(map)
folium.LayerControl().add_to(map)
 
# 지도를 표출합니다.
fig

# 지도 상에서 원하는 점의 경도, 위도를 알아낸뒤, 그 값을 인덱싱하여 가져올 수 있습니다.
# 예를 들어, 서해의 어느 한 지점의 경도 위도가 37.62293, 125.79964 라면
# 그 지점의 좌표를 다시 LCC 좌표계로 변환하여 그 값을 가져올 수 있습니다.
 
# 위·경도 좌표계에서 LCC 좌표계로 변환하는 transformer를 정의합니다.
lcc_transformer = Transformer.from_crs('EPSG:4326', source_crs)
 
# 원하는 지점의 위·경도를 입력합니다.
target_lat = 37.62293
target_lon = 125.79964
 
# 위·경도 값을 LCC 좌표계로 변환합니다.
# source_transform 연산의 반대이므로 역행렬을 곱합니다.
index_col, index_row = source_transform.__invert__().__mul__(
    lcc_transformer.transform(target_lat, target_lon)
)
 
# 지점의 안개값을 출력합니다.
print(
    f"위도 {target_lat}, 경도 {target_lon} 지점의 {kst_time} 안개값은 "
    f"{arr_fog[round(index_row-0.5), round(index_col-0.5)]:.02f}K 입니다."
)

위도 37.62293, 경도 125.79964 지점의 202404171200 안개값은 4.00K 입니다.