수치예보(Numerical Weather Prediction : NWP)란 대기현상의 역학 및 물리적 원리에 대한 지배방정식들을 컴퓨터를 활용하여 연속적으로 수치적분함으로써 현재의 대기상태를 분석하고 나아가 미래의 대기상태를 정략적으로 예측하는 일련의 과정입니다.

수치예보의 수행과정은 기상관측자료 전처리 → 자료동화 → 수치예보모델 → 후처리 → 검증 순으로 진행됩니다.

기상관측자료 전처리 수치예보모델에 입력할 관측자료를 준비하는 과정으로 전세계에서 관측한 다양한 종류의 관측자료를 수집하여 해독하고(자료수집 및 해독 과정) 수집된 자료 중 품질이 낮은 관측자료를 제거하는 과정(품질검사 과정)으로 구성되어 있습니다. 관측전처리 과정을 통해 관측자료에 포함된 오차를 평가하고 편차를 제거하는 과정은 자료동화 및 수치예보모델의 성능에 중요한 역할을 합니다.

모델은 항상 오차를 가지게 마련이며, 이 오차는 관측자료를 이용해 보정해 주어야 합니다. 전처리과정에서 수집된 관측자료를 모델에 입력해 모델이 가지고 있는 예측오차를 보정하는 과정을 자료동화라고 합니다. 자료동화를 통해 전 시간의 모델 예측결과를 보정하여 다음시간 예측을 위한 고품질의 초기장을 생산할 수 있습니다.

전통적인 자료동화	현대의 자료동화

수치예보모델을 이용해 대기의 상태 및 운동을 설명하는 대기방정식을 풀어 미래의 날씨를 계산하는 과정입니다. 대기방정식은 너무 복잡하여 정확한 해를 직접적으로 구할 수 없기 때문에, 수치해석적인 방법으로 근사적인 해를 구하게 됩니다. 수치예보모델은 대기의 큰 흐름이나 운동을 정의하는 역학과정과 작은 규모의 물리적인 현상들을 정의하는 물리과정으로 구성되어 있습니다.

역학과정	물리과정
지구를 수평· 연직 방향의 격자로 나누어 날씨의 변화를 지배하는 편미분 형태의 방정식계를 계산합니다. 따라서 격자 모양이나 간격, 방정식계를 미적분하는 방식 등에 따라 다양한 형태로 표현될 수 있으며, 이에 따라 수치모델 결과도 다르게 나타납니다.	수치모델의 역학과정에서 명시적으로 표현하지 못하는 모델 격자보다 작은 규모(아격자 규모)의 현상들은 수치모델의 물리과정에서 계산됩니다. 즉, 수치모델의 물리과정에서는 대기복사, 강수, 난류 등과 같이 역학과정에서 표현하기 어려운 지면/해양/대기의 물리적 현상의 격자 평균 효과를 통계적, 물리적 가정들을 사용하여 각 격자지점에서 적절히 모의하며, 이를 모수화라고 합니다.

수치모델이 생산한 디지털자료를 가시화하여 일기도와 가이던스를 생산하는 가시화 과정과 모델의 오차를 통계적으로 보정해 더 정확하고 상세한 예측결과를 제공하는 통계과정 등이 있습니다.

가시화	예보가이던스
수치예보모델 산출 결과는 다양한 변수가 숫자 형태로 예측시간별로 저장되어 있어 바로 이해하기 어렵습니다. 가시화는 이러한 숫자 형태의 자료를 보기 쉽고 이해하기 쉽게 그림이나 표 등의 시각적 형태로 표현하는 과정입니다.	수치예보모델로 예측된 자료를 과학적으로 처리하여 수요자(예보관)가 필요로 하는 정보로 가공한 결과물입니다. 수치예보모델은 다양한 시간, 공간해상도를 가지고 수행되지만 최종적으로는 동네예보(5km×5km간격, 3일) 또는 중기예보(15개 구역, 10일)의 형태로 제공되기 때문에 객관예보의 생산과 효율성 도모를 위하여 예보가이던스의 생산과 제공이 필요합니다. 또한 수치예보모델은 예측기간과 예측간격도 다양하여, 1시간 또는 3시간 간격의 기온과 강수량을 예측하지만 기상청의 예보는 낮 최고기온, 아침 최저기온, 6시간 강수량 등의 형태로 제공되며, 운량은 예측되지만 하늘상태가 예측되지는 않아 기상청 예보 정의에 맞는 예보요소로 가공하는 작업도 이루어집니다.

수치예보모델이 가지는 오차가 얼마나 되는지 수치로 평가하는 과정이다. 모델의 결과를 관측자료를 기준으로 검증하는 것을 관측검증이라 하며, 자료동화를 통해 생산된 고품질의 분석장을 기준으로 모델을 검증하는 것을 분석검증이라 한다.

관측검증	분석검증
세계기상기구(WMO)에서 매년 업데이트해서 제공하는 약 700여개의 고층관측(라디오존데) 지점의 관측자료를 이용하여 표준등압면의 기상요소를 검증하는 것을 말합니다. 관측지점에서 가장 가까운 모델 격자점을 이용하여 검증합니다.	수치예보모델에서 현재의 대기상태에 가장 가깝게 모의하여 생성한 분석장과 예측장을 비교하는 것을 말합니다. 모델 격자에 있는 분석장과 예측장을 일정한 크기 (WMO 기준: 1.5°×1.5°)의 검증격자로 변환하여 검증하며, 이때 해당되는 면적크기에 따라 가중 평균한 값을 적용합니다.

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

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

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

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

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

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

한국형수치예보모델은 2011년부터 개발을 시작하여 2020년 4월부터 날씨정보 생산에 활용하였습니다.

한국형수치예보모델은 육면체구 격자 기반의 전지구모델로 분광요소법 기반 수치계의 비정역학 코어와 고해상도 모델에 적합한 물리과정으로 구성되어 있습니다. 육면체구 격자는 위경도 격자의 극 특이점 문제를 해결하고 격자 거리가 비교적 균일하도록 고안된 구면 격자입니다(Sadourny, 1972). 육면체구 격자체계는 내접하는 정육면체의 각 면에서 직교좌표계를 형성한 다음 구 표면으로 투영시킨 형태의 격자로 6개의 패널을 가지는데, 각 패널에서 x축 또는 y축 한 방향으로의 요소 개수를 ne라고 하고, 각 요소에서 한 방향으로의 가우스 구적점 개수를 np라고 합니다. 아래 그림 2.3.1은 ne=10, np=4인 육면체구 격자와 한반도가 패널의 중심에 위치하도록 회전된 육면체구 격자로, 육면체구의 한 면을 차지하는 패널 안에는 10개의 요 소를 가지고 각 요소는 패널 경계 상의 점을 포함하여 4개의 점으로 구성되어 있음을 볼 수 있습니다. 현재 KIM의 해상도는 ne360np3으로, 하나의 패널에 360개의 요소와 각 요소에서 는 3개의 점으로 있으며, 이는 공간해상도 약 12 km입니다.

(a) ne=10, np=4 인 육면체구 격자, (b) 한반도가 패널의 중심에 위치하도록 회전된 육면체구 격자

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

KIM-지역모델은 2022년 5월 12일부터 현업 운영을 시작하였습니다. 동아시아 영역에 대 해 일 4회 3일(72시간) 예측을 수행하며, 수평해상도는 3km이며 예측 간격은 1시간 입니다.

r030_v040_ne36_{pres|unis|etas}_h{000~072}.yyyymmddhh.gb2

• 사이트 : 기상청 API허브 (apihub.kma.go.kr)
• 메뉴 : 수치모델 → 수치예보모델 → 한국형수치예보모델(KIM) 자료 조회

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

소스코드: vertical_cross_section.ipynb

# API 요청을 위한 requests 라이브러리와
# 지도 표출을 위한 cartopy 라이브러리,
# 연직단면도 표출을 위한 metpy 라이브러리를 설치합니다.
%pip install requests cartopy metpy

기상청 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  # 배열을 다루는 라이브러리

# KIM 자료를 txt파일로 조회하는 API url
api_url = 'https://apihub.kma.go.kr/api/typ06/cgi-bin/url/nph-kim_nc_xy_txt1'

# 요청인자로 모델 구분, 모델 기반 종류, 자료 종류, 변수명, 고도,
# 사용영역 및 격자영역, 분석시간, 예측시각, 제공형태가 있습니다.

# 1. 모델 구분으로 KIMG를 입력합니다.
group = 'KIMG'

# 2. 모델 기반 종류로 NE36을 입력합니다.
nwp = 'NE36'

# 3. 자료 종류로 등압면(P), 단일면(U), 단일면 강수/해면기압(A)이 있습니다.
# 여기선 등압면을 사용합니다.
data = 'P'

# 4. 변수명을 입력합니다.
# 여기선 사용할 변수명으로 기온(Air Temperature, T)을 사용합니다.
name = 'T'

# 5. 고도를 입력합니다.
# 여기선 등압면 24개 레벨을 모두 사용합니다.
pres_levels = np.array([
    1000, 975, 950, 925, 900, 875, 850, 800, 750, 700, 650, 600, 550,
    500, 450, 400, 350, 300, 250, 200, 150, 100, 70, 50
])

# 6. 사용영역을 전체(F), 일부(S), LCC변환 일부(R)중 하나를 입력합니다.
# 여기선 동아시아 일부 영역을 사용하므로 S를 입력합니다.
# 동아시아 영역의 좌표를 입력하기 위해 전체 영역이 사용하는 좌표계를 이해해야 합니다.
# 전체 영역은 위도[-89.9375, 89.9375]를 1440개로 나누고 (0.125간격),
# 경도[0, 360)를 2880개로 나눈 격자 좌표(0.125간격)를 사용합니다.
# 동아시아 영역이 나타내는 범위를 위도 25.8125 ~ 50.3125, 경도 104 ~ 148로 가정합니다.
# 마지막으로 격자영역은 xmin, ymin, xmax, ymax 형태로 입력합니다.
map = 'S'
degree_interval = 0.125
lat_1, lon_1, lat_2, lon_2 = 25.8125, 104, 50.3125, 148

# 격자 좌표를 계산하기위해 전체영역의 시작 위도를 더해줍니다.
# 그리고 첫 시작점을 1부터 시작시키며 정수로 변환합니다.
sub = (
    np.array([lon_1, lat_1 + 89.9375, lon_2, lat_2 + 89.9375]) / degree_interval + 1
).astype(int)

# 7. 분석시간을 년월일시(UTC)로 입력합니다.
# 예시로 2023년 8월 9일 06시(UTC)를 입력합니다.
tmfc_utc = '2023080906'

# 8. 예측시각을 입력합니다.
# 예시로 24시간 후의 예측을 입력합니다.
hf = 24

# 9. 마지막으로 제공형태를 ASCI(A), Binary(B) 중 하나를 입력합니다.
# 활용 편의를 위해 이진형태의 파일을 선택합니다.
disp = 'B'

# 파일 명명 규칙에 따라 24개 레벨의 파일 이름을 저장할 배열을 정의합니다.
# API 응답을 저장할 파일 형식은 bin 파일로 저장합니다.
file_name_list = np.array([
    f"{group.lower()}_{nwp.lower()}_pres_{name}_p{(level % 1000):03d}_h{hf:03d}.{tmfc_utc[:10]}.bin"
    for level in pres_levels
])

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

     array(['kimg_ne36_pres_T_p000_h024.2023080906.bin',
            'kimg_ne36_pres_T_p975_h024.2023080906.bin',
            'kimg_ne36_pres_T_p950_h024.2023080906.bin',
            'kimg_ne36_pres_T_p925_h024.2023080906.bin',
            'kimg_ne36_pres_T_p900_h024.2023080906.bin',
            'kimg_ne36_pres_T_p875_h024.2023080906.bin',
            'kimg_ne36_pres_T_p850_h024.2023080906.bin',
            'kimg_ne36_pres_T_p800_h024.2023080906.bin',
            'kimg_ne36_pres_T_p750_h024.2023080906.bin',
            'kimg_ne36_pres_T_p700_h024.2023080906.bin',
            'kimg_ne36_pres_T_p650_h024.2023080906.bin',
            'kimg_ne36_pres_T_p600_h024.2023080906.bin',
            'kimg_ne36_pres_T_p550_h024.2023080906.bin',
            'kimg_ne36_pres_T_p500_h024.2023080906.bin',
            'kimg_ne36_pres_T_p450_h024.2023080906.bin',
            'kimg_ne36_pres_T_p400_h024.2023080906.bin',
            'kimg_ne36_pres_T_p350_h024.2023080906.bin',
            'kimg_ne36_pres_T_p300_h024.2023080906.bin',
            'kimg_ne36_pres_T_p250_h024.2023080906.bin',
            'kimg_ne36_pres_T_p200_h024.2023080906.bin',
            'kimg_ne36_pres_T_p150_h024.2023080906.bin',
            'kimg_ne36_pres_T_p100_h024.2023080906.bin',
            'kimg_ne36_pres_T_p070_h024.2023080906.bin',
            'kimg_ne36_pres_T_p050_h024.2023080906.bin'], dtype='<U41')

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

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

            # API 요청인자와 함께 API 요청
            response = session.get(api_url, params=api_parameters)

            # 잘못된 응답을 받거나 짧은 에러메세지를 응답으로 받은 경우 에러 메세지 출력
            if response.status_code != 200 or len(response.content) < 300:
                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]} 파일 다운로드 완료')

  kimg_ne36_pres_T_p000_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p975_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p950_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p925_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p900_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p875_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p850_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p800_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p750_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p700_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p650_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p600_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p550_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p500_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p450_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p400_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p350_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p300_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p250_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p200_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p150_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p100_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p070_h024.2023080906.bin 파일 다운로드 완료
  kimg_ne36_pres_T_p050_h024.2023080906.bin 파일 다운로드 완료

다운로드받은 24개 레벨의 파일을 읽어 하나의 xarray 데이터셋으로 합쳐서 다룹니다.

import xarray as xr # 자료를 다루기 위한 라이브러리

# 24개 레벨의 파일을 읽어 xarray 데이터셋을 만듭니다.
# xarray 데이터셋은 크게 variables, coordinates, attributes로 구성되어 있습니다.
# attributes는 데이터에 사용된 각종 정보들을 나타냅니다. 
xr_array = xr.Dataset(
    # variable은 실제 데이터를 담고 있는 변수로 변수이름과 (차원, 데이터)를 dictionary 형태로 입력합니다.
    # 데이터는 파일로부터 읽어들이며 이진 파일의 앞의 4byte, 4byte는 데이터의 가로, 세로 격자수를 나타내며
    # 나머지 데이터는 4byte 소수형으로 읽어들입니다.
    data_vars={
        'temperature': (
            ['isobaric', 'lat', 'lon'], (data_array:=np.stack([
                np.fromfile(f, offset=8, dtype=np.float32).reshape(
                    sub[3] - sub[1] + 1, sub[2] - sub[0] + 1
                )
                for f in file_name_list]))
            )
        },
    # coordinates는 variable의 차원을 나타내는 변수입니다.
    # 앞서 variable의 차원부분에 입력한 이름에 맞게 실제 차원 수를 dictionary 형태로 입력합니다.
    coords={
        'isobaric': pres_levels,
        'lat': np.linspace(lat_1, lat_2, num=data_array.shape[1]),
        'lon': np.linspace(lon_1, lon_2, num=data_array.shape[2]),
    },
    # attributes는 데이터에 사용된 각종 정보들을 나타냅니다.
    # API 요청시 사용한 인자들을 dictionary 형태로 입력합니다.
    attrs={
        'group': group,
        'nwp': nwp,
        'data': data,
        'name': name,
        'map': map,
        'sub': ','.join([str(x) for x in sub]),
        'tmfc': tmfc_utc,
        'hf': hf,
    }
)

# 생성한 xarray 데이터셋을 출력합니다.
xr_array

  <xarray.Dataset> Size: 7MB
  Dimensions:      (isobaric: 24, lat: 197, lon: 353)
  Coordinates:
  * isobaric     (isobaric) int64 192B 1000 975 950 925 900 ... 150 100 70 50
  * lat          (lat) float64 2kB 25.81 25.94 26.06 26.19 ... 50.06 50.19 50.31
  * lon          (lon) float64 3kB 104.0 104.1 104.2 104.4 ... 147.8 147.9 148.0
  Data variables:
      temperature  (isobaric, lat, lon) float32 7MB 305.4 305.2 ... 216.9 216.9
  Attributes:
      group:    KIMG
      nwp:      NE36
      data:     P
      name:     T
      map:      S
      sub:      833,927,1185,1123
      tmfc:     2023080906
      hf:       24

연직단면도를 표출하기 앞서, 기압 레벨에 따른 한반도 지역의 기온 분포도를 표출합니다.

import matplotlib.pyplot as plt     # 이미지 표출에 사용되는 라이브러리
import cartopy.crs as ccrs          # 지도 표출에 사용되는 라이브러리
import cartopy.feature as cfeature  # 지도 표출에 사용되는 라이브러리
from cartopy.crs import PlateCarree # 지도 표출에 사용되는 라이브러리

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

# 비교할 분포도는 1000hPa과 200hPa입니다.
# 분포도 내의 최대, 최소 값을 이용해 값의 범위를 구합니다.
isobaric_1 = 1000
isobaric_2 = 200
isobaric_1_arr = xr_array['temperature'][np.where(pres_levels == isobaric_1)[0][0]]
isobaric_2_arr = xr_array['temperature'][np.where(pres_levels == isobaric_2)[0][0]]
vmin1, vmax1 = isobaric_1_arr.min() - isobaric_1_arr.min() % 5, isobaric_1_arr.max() - isobaric_1_arr.max() % 5 + 5
vmin2, vmax2 = isobaric_2_arr.min() - isobaric_2_arr.min() % 5, isobaric_2_arr.max() - isobaric_2_arr.max() % 5 + 5

# 이미지를 그릴 크기를 지정합니다.
kw = {'projection': PlateCarree()}
fig, axes = plt.subplots(1, 2, figsize=(15, 5), subplot_kw=kw)
fig.tight_layout(pad=3)

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

# 이미지를 각 영역에 그립니다.
im1 = axes[0].imshow(isobaric_1_arr, extent=[lon_1, lon_2, lat_1, lat_2], cmap='turbo', origin='lower', vmin=vmin1, vmax=vmax1)
im2 = axes[1].imshow(isobaric_2_arr, extent=[lon_1, lon_2, lat_1, lat_2], cmap='turbo', origin='lower', vmin=vmin2, vmax=vmax2)

# 이미지에 해안선과 경계를 표시합니다.
axes[0].coastlines()
axes[1].coastlines()
axes[0].add_feature(cfeature.STATES.with_scale('50m'), edgecolor='k', alpha=0.2, zorder=0)
axes[1].add_feature(cfeature.STATES.with_scale('50m'), edgecolor='k', alpha=0.2, zorder=0)

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

# 색상표에 표시될 글자 크기, 제목, 범위를 설정합니다.
cbar1 = fig.colorbar(im1, cax=cax1)
cbar2 = fig.colorbar(im2, cax=cax2)
cbar1.ax.tick_params(labelsize=8)
cbar1.ax.set_yticks(np.linspace(vmin1, vmax1, 11))
cbar1.ax.set_title('K', fontsize=8)
cbar2.ax.tick_params(labelsize=8)
cbar2.ax.set_yticks(np.linspace(vmin2, vmax2, 11))
cbar2.ax.set_title('K', fontsize=8)

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

마지막으로 연직단면도를 표출합니다.

from metpy.interpolate import cross_section # 연직단면도 계산을 위한 함수
from pyproj import CRS  # 좌표계 형식을 사용하기 위한 라이브러리


# 연직단면도를 위한 두 점의 좌표를 위·경도 좌표로 입력합니다.
start = (36.9, 120.7)
end = (37.5, 133.4)

# metpy를 활용하기 위해 데이터가 좌표계(CRS)를 인식할 수 있도록 CRS를 등록(assign_crs)합니다.
xr_array = xr_array.metpy.assign_crs(CRS('EPSG:4326').to_cf())

# 데이터를 두 점을 연결하는 선으로 자른 연직단면도로 변환합니다.
cross = cross_section(xr_array, start, end).set_coords(('lat', 'lon'))

# 변환된 연직단면도 데이터를 표출합니다.
# 기온 변수의 차원이 lat, lon좌표에서 index(100개)로 변환된 것을 확인할 수 있습니다.
cross

  <xarray.Dataset> Size: 22kB
  Dimensions:      (isobaric: 24, index: 100)
  Coordinates:
  * isobaric     (isobaric) int64 192B 1000 975 950 925 900 ... 150 100 70 50
      metpy_crs    object 8B Projection: latitude_longitude
      lon          (index) float64 800B 120.7 120.8 121.0 ... 133.1 133.3 133.4
      lat          (index) float64 800B 36.9 36.91 36.93 36.94 ... 37.5 37.5 37.5
  * index        (index) int64 800B 0 1 2 3 4 5 6 7 ... 92 93 94 95 96 97 98 99
  Data variables:
      temperature  (isobaric, index) float64 19kB 303.6 303.5 ... 211.0 210.7
  Attributes:
      group:    KIMG
      nwp:      NE36
      data:     P
      name:     T
      map:      S
      sub:      833,927,1185,1123
      tmfc:     2023080906
      hf:       24

# 연직단면도를 그릴 크기를 지정합니다.
fig, ax = plt.subplots(figsize=(16, 9))

# 기온을 표출하기 위해 등치선을 채우는 방식(filled contour, contourf)을 사용합니다.
# 레벨은 최소값과 최대값을 이용해 구합니다.
temperature_contour = ax.contourf(
    cross['lon'], cross['isobaric'], cross['temperature'], cmap='turbo',
    levels=np.arange(
        cross['temperature'].min() - cross['temperature'].min() % 10,
        cross['temperature'].max() - cross['temperature'].max() % 10 + 10, 5
    )
)
# 등치선의 색상표를 표출합니다.
temperature_colorbar = fig.colorbar(temperature_contour, pad=0)
temperature_colorbar.ax.set_title('K')

# Y축을 로그 스케일(log scale)로 변환합니다. 또, 표출할 범위를 지정합니다.
ax.set_yscale('symlog')
ax.set_ylim(cross['isobaric'].max(), cross['isobaric'].min())
ax.set_yticks(np.arange(1000, 50, -100))
ax.set_yticklabels(np.arange(1000, 50, -100))

# 데이터의 좌표계를 정의하고, 지도를 오른쪽 상단에 표시합니다.
data_crs = cross['temperature'].metpy.cartopy_crs
ax_inset = fig.add_axes([0.85, 0.65, 0.25, 0.25], projection=data_crs)

# start, end 두 점의 위·경도 좌표를 지도가 사용하는 좌표계로 변환하여 지도에 표출합니다.
# 또, 앞서 그린 분포도를 표출합니다.
endpoints = data_crs.transform_points(
    ccrs.Geodetic(), *np.vstack([start, end]).transpose()[::-1]
)
ax_inset.imshow(isobaric_2_arr, origin='lower', extent=[lon_1, lon_2, lat_1, lat_2], cmap='turbo', vmin=vmin2, vmax=vmax2)
ax_inset.scatter(endpoints[:, 0], endpoints[:, 1], c='k', zorder=2)
ax_inset.plot(cross['lon'], cross['lat'], c='k', zorder=2)

# 지도에 해안선 및 경계를 표출합니다.
ax_inset.coastlines()
ax_inset.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='k', alpha=0.2, zorder=0)

# 단면도의 제목과 축의 라벨을 지정합니다.
ax_inset.set_title('')
ax.set_title(f'Forecast after {hf} hours from {tmfc_utc} Cross-Section {start} to {end} Temperature(K)')
ax.set_ylabel('Pressure (hPa)')
ax.set_xlabel('Longitude (degrees east)')
temperature_colorbar.set_label('Temperature (dimensionless)')

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