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

수치모델은 한정된 컴퓨터 자원을 이용하여야 하고 현재 관측시스템의 관측 한계 때문에 무제한으로 작게 분해하여 모의할 수 없습니다.

여기서 공간해상도란 대기를 컴퓨터 계산에 적합하게 나타내는 격자점과 격자점 사이의 거리 즉, 격자점이 대표하는 공간규모(수백m~수백㎞)로 이해할 수 있으며 이를 수평해상도라 합니다. 격자점 사이의 거리가 짧을수록 수치모델이 모의할 수 있는 기상현상의 규모도 작아집니다.

일반적으로 모델의 해상도는 사용가능한 컴퓨터의 성능, 모델이 모의할 영역의 크기, 모델이 모의하고자 하는 가장 작은 기상현상의 규모 등에 의해 결정됩니다. 즉, 컴퓨터의 성능이 우수할수록, 모의영역이 작을수록 그리고 모의하고자 하는 현상의 크기가 작을수록 고해상도로 구성하게 됩니다.

모델대기는 역직방향 수십 개 층으로 나누어지는데 각 층은 지구표면에 대해 고정된 개수의 격자점을 포함합니다. 모델대기의 최상단 고도와 함께 연직 층수는 모델 대기의 연직 해상도를 결정합니다. 각 격자점에는 각 격자점이 포함하는 대기요소들의 평균 상태를 나타내는 일련의 값들을 갖습니다(아래 그림 1).

수치모델에서의 적분은 일정한 시간간격으로 연속적으로 진행되어야 하며 이 시간간격은 다시 컴퓨터의 성능에 의해 좌우됩니다. 시간간격(적분시간)이 짧을수록 보다 많은 수의 계산이 수행됩니다. 미래의 새로운 대기상태를 결정하기 위해서 각 공기덩이와 관계된 전체 방정식 계들이 이 적분시간 간격으로 동시에 계산되어야 합니다. 이로부터 새로운 힘의 장이 계산되고 다음 적분에서 공기덩이들이 어떻게 이동할지를 계산합니다.

시간분해능은 주어진 적분시간내에 공기덩이가 모델 격자 거리보다 먼 거리를 이동하지 않도록 공간분해능과 조화를 이루어야 합니다. 만약 v가 가장 빠른 속도이고 d(㎞)가 격자거리라 하면 적분시간은 d/v보다 길어서는 안됩니다. 적분이 진행되는 과정에서 각 공기덩이의 이동을 따라 대기상태를 예측하는 것이 아니라 매 적분시간마다 각 격자내 공기의 새로운 상태(기압, 온도, 습도, 바람 등)들이 계산됩니다.

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

* 기준점 위치(x, y) : UM지역(2953.737876km, 2522.911728km), UM국지(389.478893km, 618.363966km)

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

g512_v070_pres

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

* 기준점 위치(x, y) : UM지역(2953.737876km, 2522.911728km), UM국지(389.478893km, 618.363966km)

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

* 기준점 위치(x, y) : UM지역(2953.737876km, 2522.911728km), UM국지(389.478893km, 618.363966km)

# 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 100000 Pa:fcst time 24 hrs:from 202405141200

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)
)

korea_shp.zip 파일 다운로드 완료

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

(400, 428)

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()