영국 통합모델 기반으로 되어 있는 UM 전지구 예보모델은 수평으로 2560 × 1920개 격자(N1280)로 구성되어 약 10㎞의 해상도를 가지며, 연직으로 약 80㎞ 고도까지 70층으로 구성되어있습니다. 수평으로 Arakawa-C 격자를, 연직으로 Charney-Phillips 방식을 채택하고 있습니다. UM 전지구 예보모델의 경우, 00UTC, 12UTC에는 288시간 예측을 수행하고, 06UTC, 18UTC에는 87시간 예측을 수행합니다.

위 그림은 기상청이 운영하고 있는 주요 수치예보시스템 구성을 요약 한 것으로 2020년 4월부터 현업 운영 중인 한국형수치예보모델(KIM)의 안정적 성능 확보까지 기존 통합모델 기반의 수치예보시스템과 병행 운영 중입니다. 또한, 전지구예보모델인 한국형수치예보모델의 현업 운영에 따른 연계모델(앙상블, 국지, 응용모델 등)들의 순차 개발이 진행 중에 있습니다.

* 저장순서: 남→북, 서→ 동 으로 저장

g128_v070_ergl_{pres|unis|isen|etas}_h{000~288}.yyyymmddhh.gb2

• UM 변수목록: UM 변수 목록

* 저장순서: 남→북, 서→ 동 으로 저장

l015_v070_erlo_{pres|unis|isen|etas}_h{000~48}.yyyymmddhh.gb2

• UM 변수목록: UM 변수 목록

• 사이트 : 기상자료개방포털 (data.kma.go.kr)​
• 메뉴 : 데이터 → 수치모델 → 단중기예측 → 지역예보모델, 국지예보모델

• 조건 : 기간, 자료구분, 분석시간, 예측시간

• 사이트 : 기상청 API허브 (apihub.kma.go.kr)
• 메뉴 : 수치모델 → 수치예보모델 → 수치모델 경량화 다운로드(예측시간+변수+고도별), 통합모델(UM) 자료 조회

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

소스코드: skew_t_log_p.ipynb

# API 요청을 위한 requests 라이브러리와 
# grib 파일을 다루기 위한 pygrib 라이브러리
# 한반도 지형 표출을 위한 geopandas 라이브러리를 설치합니다.
# 또, 단열선도 표출을 위한 matplotlib, metpy 라이브러리,
# 좌표계 변환을 위한 pyproj 라이브러리를 설치합니다.
%pip install requests pygrib geopandas matplotlib metpy pyproj

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

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

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

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

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

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

import requests     # API 요청을 보내고 받는 라이브러리
import numpy as np  # 배열을 다루는 라이브러리

# 수치모델 경량화 파일 다운로드 API url
api_url = 'https://apihub.kma.go.kr/api/typ06/url/nwp_vars_down.php'

# 요청인자로 모델종류, 파일구분, 기준시각, 예측시각, 저장포맷이 있습니다.

# 1. 모델종류로 전구(Global, g), 지역(Region, r), 국지(Local, l) 중 하나를 입력합니다.
# 여기선 예시로 국지모델을 사용하며 1.5km의 해상도를 갖고 있어 l015로 명명합니다.
nwp = 'l015'

# 2. 파일구분으로 등압면(pres), 단일면(unis), 단일온면(isen) 중 하나를 입력하고 원하는 요소를 입력합니다.
# 여기선 등압면(pres)을 선택합니다.
# 등압면의 경우 24개 레벨로 구분하여 제공되며 여기선 모든 레벨의 등압면을 사용합니다.
file_seperator = {
    'sub': 'pres',
    'vars': 'tmpr'
}
# 제공하는 24개 레벨의 등압면 배열
pres_levels = [
    1000, 975, 950, 925, 900, 875, 850, 800, 750, 700, 650,
    600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100, 70, 50
]

# 3. 기준시각을 년월일시(UTC) 형태로 입력합니다.
# 수치모델은 일 4회 6시간 간격으로 파일을 생성하므로 00, 06, 12, 18시 중 하나를 입력합니다.
tmfc_utc = '2024051412'

# 4. 예측시각을 입력합니다.
# 국지모델은 기준시간으로부터 1시간 간격으로 최대 48시간을 예측합니다.
ef = 24

# 5. 데이터 저장포맷으로 GRIB, BIN, TEXT 중 하나를 입력합니다.
# 여기선 GRIB 파일로 지정합니다.
data_type = 'GRIB' 

# 경량화 파일 명명 규칙에 따라 24개 레벨의 파일 이름을 저장할 배열을 정의합니다.
# API 응답을 저장할 파일 형식은 gb2 파일로 저장합니다.
file_name_list = np.array([
    f"{nwp}_v070_erlo_{file_seperator['sub']}_{file_seperator['vars']}_"
    f"{file_seperator['sub'][0]}{(level % 1000):03d}_"
    f"h{ef:03d}.{tmfc_utc[:10]}.gb2" for level in pres_levels]
)

# 파일 목록을 출력합니다.
file_name_list

# 각 파일의 존재 여부를 확인합니다.
file_exists_list = np.vectorize(os.path.exists)(file_name_list)

# 저장된 파일이 없는 경우에만 API를 요청해 파일을 다운로드합니다.
if not file_exists_list.all():
    # API 요청인자들을 묶어 dictionary로 정의합니다.
    api_parameters = {
        'nwp': nwp,
        'tmfc': tmfc_utc,
        'ef': ef,
        'dataType': data_type,
        'authKey': my_api_key,
        **file_seperator
    }
    # 파일 다운로드를 위한 세션을 만듭니다.
    with requests.Session() as session:
        # 존재하지 않는 파일에 대해서 반복합니다.
        for idx in np.where(~file_exists_list)[0]:
            # API 요청인자에 레벨을 입력합니다.
            api_parameters[file_seperator['sub']] = pres_levels[idx]

            # API 요청인자와 함께 API 요청
            response = session.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_list[idx], 'wb') as f:
                    f.write(response.content)
                print(f'{file_name_list[idx]} 파일 다운로드 완료')

단열선도는 고층의 기상 자료를 분석하는데 사용되는 그래프입니다. 대기의 건조상태, 구름의 분포, 대기의 불안정 여부 등을 분석할 수 있습니다.

먼저, 데이터를 읽고 한반도 분포도를 표출하고, 수치 모델이 예측한 대기 상태를 분석하기 위해 단열선도를 표출합니다.

# grib 파일을 다루는 라이브러리
import pygrib

# 저장된 경량화 수치모델을 읽습니다.
# 분포도 예시로 고도가 낮은 등압면 1000hPa과 고도가 높은 등압면 50hPa를 읽습니다.
with pygrib.open(file_name_list[0]) as grbs:
    # 파일 내의 저장된 메세지들을 전부 출력합니다.
    for grb_message_1000 in grbs:
        print(grb_message_1000)

with pygrib.open(file_name_list[-1]) as grbs:
    # 파일 내의 저장된 메세지들을 전부 출력합니다.
    for grb_message_50 in grbs:
        print(grb_message_50)

1:Temperature:K (instant):lambert:isobaricInhPa:level 5000 Pa:fcst time 24 hrs:from 202405141200

grib 파일은 각각의 메세지를 파일 내에 저장하며 메세지 별로 키(Key)가 존재합니다.

인덱스, 이름, 단위, 격자종류, 기준시각, 예측시각 등이 key로 저장됩니다.

메세지의 data에는 실제 데이터, 격자의 위도와 경도 값이 저장됩니다.

실제 데이터를 이용해 기온 분포도를 표출합니다.

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

또, 데이터에 대한 참고자료를 다운로드 받습니다.

import geopandas as gpd # shape 파일을 다루는 라이브러리
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)

# 국가 경계를 구분한 shp파일은 ARD_NAION_AS으로
# 이 파일을 읽어 수치모델이 갖는 좌표계인 LCC 도법으로 변환합니다.
lcc_crs =(
    "+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=126 "
    "+x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
)
gdf = gpd.read_file(
    os.path.join(extract_directory, 'ARD_NAION_AS.shp')
).to_crs(lcc_crs)

# 1m 단위의 도법에서 수치모델의 해상도 단위(1.5km)에 맞춥니다.
# 참고자료에 따르면 기준점 위치의 경우 x: 389.478893km, y: 618.363966km이며
# 각 격자마다 1.5km의 간격을 가집니다.
resolution = 1500
origin_x = 389.478893/1.5
origin_y = 618.363966/1.5
gdf.geometry = gdf.geometry.scale(
    xfact=1/resolution, yfact=1/resolution, zfact=1.0,
    origin=(origin_x, origin_y)
)

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

# 분포도에서 색상표 위치, 크기를 조절하는 함수
from mpl_toolkits.axes_grid1 import make_axes_locatable

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

# 이미지의 제목을 지정합니다.
fig.suptitle(f"Forecast after {ef} hours from {tmfc_utc}")
axes[0].set_title(f"Temperature at 1000hPa")
axes[1].set_title(f"Temperature at 50hPa")

# 국가간의 경계선을 검정색으로 그립니다.
gdf.boundary.plot(ax=axes[0], color='black', linewidth=0.5)
gdf.boundary.plot(ax=axes[1], color='black', linewidth=0.5)

# 이미지를 각 영역에 그립니다.
im1 = axes[0].imshow(grb_message_1000.data()[0], cmap='turbo', origin='lower')
im2 = axes[1].imshow(grb_message_50.data()[0], cmap='turbo', origin='lower')

# 색상표의 크기와 위치를 조절합니다.
divider1 = make_axes_locatable(axes[0])
cax1 = divider1.append_axes('right', size='5%', pad=0)
divider2 = make_axes_locatable(axes[1])
cax2 = divider2.append_axes('right', size='5%', pad=0)

# 색상표에 표시될 글자 크기 및 제목을 설정합니다.
cbar1 = fig.colorbar(im1, cax=cax1)
cbar2 = fig.colorbar(im2, cax=cax2)
cbar1.ax.tick_params(labelsize=8)
cbar1.ax.set_title('K', fontsize=8)
cbar2.ax.tick_params(labelsize=8)
cbar2.ax.set_title('K', fontsize=8)

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

임의의 지역에서의 단열선도를 그리기 위해 지역의 위·경도 값을 격자 좌표로 변환합니다.

그리고 격자 좌표에서의 기압-기온 값을 이용해 단열선도를 표출합니다.

# 지점의 좌표계 변환을 위한 transformer
from pyproj.transformer import Transformer

# 위·경도 좌표계에서 LCC 좌표계로의 변환하는 transformer를 정의합니다.
transformer = Transformer.from_crs('EPSG:4326', lcc_crs)

# 원하는 지역의 위도 경도 값을 입력합니다.
# 예시로 강원도의 한 지점의 값을 입력합니다.
target_lat, target_lon = 38.189989, 128.474987

# 데이터의 격자 좌표로 변환합니다.
transformed_x, transformed_y = transformer.transform(target_lat, target_lon)
transformed_x = round(transformed_x / resolution + origin_x)
transformed_y = round(transformed_y / resolution + origin_y)

# 변환된 격자 좌표를 표출합니다.
transformed_x, transformed_y

from metpy.plots.skewt import SkewT # 단열선도 표출을 위한 SkewT
from metpy.units import units       # 단열선도 표출을 위한 단위 지정

# 기압 레벨에 따른 각 파일들을 읽고 격자 좌표에 해당하는 기온 값을 배열에 저장합니다.
# 단위는 절대온도(K)입니다.
temperature_array = np.array(
    [
        pygrib.open(f).message(1).data()[0][transformed_y, transformed_x]
        for f in file_name_list
    ]
) * units.degree_Kelvin

# 단열선도를 그릴 크기를 지정합니다.
fig = plt.figure(figsize=(9, 9))

# 단열선도를 정의합니다.
skew = SkewT(fig)

# 제목을 지정합니다.
skew.ax.set_title(f"Skew T forecast {ef} hour after from {tmfc_utc}(UTC) in {target_lat}N, {target_lon}E")

# y축(기압축)의 범위를 레벨에 맞게 지정합니다.
skew.ax.set_ylim(pres_levels[0], pres_levels[-1])

# 단열선도에서 건조단열선, 포화단열선(습윤단열선), 포화혼합비선을 그립니다.
skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()

# 기압과 기온을 이용해 단열선도를 그립니다.
skew.plot(pres_levels, temperature_array, color='red')

# 그래프의 x, y축 라벨을 입력합니다.
plt.xlabel('Temperature (℃)')
plt.ylabel('Pressure (hPa)')

# 단열선도를 표출합니다.
plt.show()