[MetPy] 단열선도 그리기

- 개요

전 대기과학 분야의 그림 중에 단열선도가 직접 그리기 가장 어렵다고 생각합니다.

파이썬을 쓰기 전에도 물론 직접 그리진 않고, NCL에서 제공하는 단열선도로 그림을 그렸죠.

파이썬을 쓰는 지금은 MetPy 라이브러리로 단열선도를 그리면 됩니다. 

 

이번 포스트에서는 MetPy 라이브러리의 parcel_profile_with_lcl_as_dataset 함수를 이용해서 단열선도를 그려보겠습니다.

아래의 예제 코드를 참고했으며, 저는 온도, 이슬점과 관련된 변수만 그릴 것입니다.

단열선도 오른 편에 풍향, 풍속도 나타낼 수 있습니다.

https://unidata.github.io/MetPy/latest/examples/plots/Sounding_LCL_Dataset.html#sphx-glr-examples-plots-sounding-lcl-dataset-py

 

- 함수 설명 및 데이터 처리

1. parcel_profile_with_lcl_as_dataset 함수 설명

metpy.calc.parcel_profile_with_lcl_as_dataset(pressure, temperature, dewpoint)

이 함수에서는 pressure(기압), temperature(온도), dewpoint(이슬점)의 세 가지 입력 변수가 필요합니다.

세 변수 모두 기압의 층 수만큼의 길이를 가진 array 형태입니다.

여기서 기압은 하층부터 상층으로, 즉 내림차순으로 정렬되어있어야 합니다.

출력 변수로는 그림을 그릴 때 필요한 단열선도의 온도, 이슬점 프로파일이 나옵니다.

 

2. 입력 변수의 클래스

MetPy에서 쓰는 변수의 클래스는 단순히 integer, float 같은 것이 아니라

"pint.Quantity"라는 것을 사용합니다.

이 클래스에서는 value(값)과 unit(유닛)이 정의되어 있고, MetPy 라이브러리의 모듈은 이러한 단위를 고려해서 알아서 여러 값들을 계산해줍니다.

 

3. 데이터 처리

전 위성 데이터의 온도와 mixing ratio의 프로파일을 가지고 있으므로 이슬점을 따로 계산했습니다.

그리고 변수들을 pint.Quantity로 바꿔줍니다.

'metpy.units.단위'를 곱하면 pint.Quanity 클래스로 바뀌니 단위를 정확히 적어주세요.

from metpy.calc import saturation_vapor_pressure
from metpy.calc import vapor_pressure
from metpy.calc import dewpoint
from metpy.units import units

# mixing ratio로 vapor pressure를 구함
vp_prof = [vapor_pressure(p[i] * units.hPa, w_prof[i] * units('g/kg')).to('hPa').magnitude for i in range(len(p))]

# vapor pressure로 이슬점을 구함
td_prof = [dewpoint(vp_prof[i] * units.hPa).magnitude for i in range(len(p))]

"""
p는 기압, T는 온도, td_prof은 위에서 구한 이슬점

p와 T는 numpy.ndarray 클래스
Td_prof는 list 클래스

metpy.units을 곱하는 식으로 변수를 선언하면 qint.Quantity 클래스로 변환 가능
단 input의 클래스를 list로 바꿔주어야 함
"""
p = p.tolist() * units.hPa
T = t_prof.tolist() * units.degC
Td = td_prof * units.degC

 

- 단열선도 그리기

앞서 구한 p, T, Td 변수로 단열선도를 그려보겠습니다.

MetPy는 matplotlib 라이브러리를 활용해 단열선도를 그립니다.

import metpy.calc as mpcalc
import matplotlib.pyplot as plt

ds = mpcalc.parcel_profile_with_lcl_as_dataset(p, T, Td)

fig = plt.figure(figsize=(10, 6))
skew = SkewT(fig, rotation=45)

# 온도 선 그리기
skew.plot(ds.isobaric, ds.ambient_temperature, 'r')
# 이슬점 선 그리기
skew.plot(ds.isobaric, ds.ambient_dew_point, 'g')
# 지상의 공기를 강제로 상승 시켰을 때의 온도
skew.plot(ds.isobaric, ds.parcel_temperature.metpy.convert_units('degC'), 'black')

# 단열선도의 배경을 그리는 부분
pressure = np.arange(1000, 499, -50) * units('hPa')
mixing_ratio = np.array([0.1, 0.2, 0.4, 0.6, 1, 1.5, 2, 3, 4,
                        6, 8, 10, 13, 16, 20, 25, 30, 36, 42]).reshape(-1, 1) * units('g/kg')
# 건조 단열선
skew.plot_dry_adiabats(t0=np.arange(233, 533, 10) * units.K, alpha=0.25,
                       colors='orangered', linewidths=1)
# 습윤 단열선      
skew.plot_moist_adiabats(t0=np.arange(233, 400, 5) * units.K, alpha=0.25,
                         colors='tab:green', linewidths=1)
# mixing ratio 선        
skew.plot_mixing_lines(pressure=pressure, mixing_ratio=mixing_ratio, linestyles='dotted',
                       colors='dodgerblue', linewidths=1)
                       
# y축 1050 hPa부터 200 hPa, 위성 자료의 최하층 기압이 1050 hPa이라 이렇게 설정
skew.ax.set_ylim(1050,200)
# x축 영하 10도부터 영상 40도
skew.ax.set_xlim(-10, 40)

plt.title('red: ambient_temperature, green:dew_point_temperature')

plt.show()

 

단열선도 결과

 

 

- 검은선의 해석

검은선은 1050 hPa에 위치한 공기를 강제로 상승시켰을 때의 온도입니다.

빨간선은 온도, 초록선은 이슬점입니다.

가장 하층인 1050 hPa을 보시면 검은선의 온도는 빨간선의 위치와 같습니다.

이 공기를 상승시키면 처음에는 아직 포화가 안되었으므로 건조단열선의 기울기만큼 온도가 변합니다.

그래서 1050 hPa부터 1000 hPa의 검은선은 건조단열선과 평행합니다.

 

다음으로 1050 hPa에서의 초록선의 위치를 봅시다.

이슬점은 공기가 가진 수증기라고 생각하면 되는데 공기가 상승해도 수증기량은 거의 변하지 않습니다.

그러므로 1050 hPa 선과 초록 선이 만나는 지점에서 mixing ratio선(파란 점선)과 평행하게 선을 그어봅시다.(운이 좋게도 여기선 초록선과 일치함)

그러면 1000 hPa에서 방금 말한 mixing ratio선과 평행한 선이 검은선과 만납니다.

 

이 지점에서는 공기가 포화된 것으므로 1000 hPa는 상승할 때 습윤단열선과 평행하게 온도가 변하게 됩니다.

 

공기를 강제로 상승시켜 포화되는 지점을 LCL(Lifting Condensation Level)이라고 합니다)