[MetPy] Vertical cross section, 연직단면 그리기

- 개요

기상 자료의 차원은 시간, 연직층, 위도, 경도로 많아야 4개입니다.

특정 시간, 특정 층에서 기상 변수의 수평 분포(위도, 경도)는 matplotlib의 contour, contourf만 잘 쓰면 쉽게 그릴 수 있습니다. 하지만 특정 경로에서 대기 변수의 연직 분포를 알고 싶다면 난이도가 많이 올라갑니다.

정공법으로 연직단면을 그리려면 먼저 경로를 정하고, 연직층마다 이 경로의 좌표(위도, 경도)에 맞게 interpolation을 하여 각 층에서 경로의 좌표에 따른 변수값을 계산해야 합니다.

 

이런 코드를 짠다니 생각만 해도 머리가 아픕니다.

다행히도 MetPy는 특정 경로의 연직단면을 얻을 수 있는 함수를 제공하고 있습니다.

 

- 자료 읽기

이번에는 NCEP2 재분석 자료의 30년 평균 월평균 온도를 사용했습니다.

기본적으로 metpy는 xarray 자료 구조를 사용하기 때문에 xarray 라이브러리로 파일을 읽습니다.

import xarray as xr

fn = 'air.mon.ltm.1991-2020.nc'
data = xr.open_dataset(fn)

 

자료의 구조를 보시면 'air'가 온도 변수이며, 차원은 time(시간), level(연직), lat(위도), lon(경도)입니다.

NCEP2 재분석 자료의 구조

 

 

- 연직단면을 자르기 위해 새로운 xarray 생성

연직단면을 만들어주는 함수를 이용하기 위해서는

1. coords(위도, 경도와 같은 좌표)

2. crs(투영하는 좌표계)

이렇게 두 정보가 필요합니다.

 

재분석 자료를 저장한 data 변수에서는 저 두 값을 설정하여도 연직단면 함수가 동작하지 않았습니다.

원래 불가능한 것인지 제가 그냥 못하는 것인지 잘 모르겠네요.

그래서 연직단면 함수를 이용하기 위해 새로운 xarray를 생성했습니다.

ds = xr.Dataset(
    {'air': (("time", "level", "lat", "lon"), data.air.data)}, 
    #{'변수 이름': ((변수 차원), 변수값)}
    coords = {
        'time': data.time.data, # 시간
        'level': data.level.data, # 연직층
        'lat': data.lat.data, # 위도
        'lon': data.lon.data # 경도
    }
)

ds = ds.metpy.assign_crs(
    grid_mapping_name="latitude_longitude"
)

 

data.air.data는 data의 air 변수에서 값을 가져오겠다는 뜻입니다. pandas의 values와 비슷한 기능을 합니다.

 

다음으로 NCEP2 재분석 자료의 투영법을 생각해봅시다.

위도, 경도는 1차원으로 등간격이며 cartopy 라이브러리 crs의  PlateCarree 투영과 같다고 보시면 됩니다.

이는 assign_crs에서 grid_mapping_name을 latitude_longitude로 선언한 것과 같습니다.

 

만약 다른 투영법으로 된 자료라면 투영법과 관련된 다른 변수들도 선언해줄 필요가 있습니다.

투영법과 관련된 정보는 아래의 링크에 자세히 적혀있습니다.

http://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/cf-conventions.html#appendix-grid-mappings

 

- 연직단면 자르기

새로운 xarray를 만들었으니 연직단면을 잘라봅시다.

start와 end 변수는 위도, 경도 값으로 경로의 시작과 끝점을 의미합니다.

이 경로의 연직단면을 구해봅시다.

from metpy.interpolate import cross_section

start = (35., 125.) # 시작점
end = (40., 130.) # 끝점

cross = cross_section(ds, start, end).set_coords(('lat', 'lon'))

 

다음으로 cross 변수의 정보를 보면 index는 0부터 99까지로 start부터 end까지의 경로를 100개의 점으로 나누었다는 뜻입니다.

lon과 lat의 값을 보면 시작점부터 끝점까지로 정의되어있고 100개의 점에 대한 위도, 경도 값입니다.

 

cross 변수의 정보

 

- 연직단면 그리기

cross['level][0:5]는 1000 hPa부터 600 hPa까지를 의미하고 cross['air'][0, 0:5, :]는 1월의 온도값입니다.

y축을 연직층, x축을 위도로 설정했습니다.

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot()
shading_vert = ax.contourf(cross['lat'], cross['level'][0:5], cross['air'][0,0:5,:])
ax.invert_yaxis() # 연직층이 기압이라 위로 갈수록 감소해야 하므로 y축을 뒤집어야 함

cbar = fig.colorbar(shading_vert)

연직단면

마지막으로 연직단면이 잘 잘라졌는지 확인해봅시다.

아래는 850 hPa 온도의 수평 분포를 그리고, 연직단면을 자른 경로를 그리는 코드입니다.

levels = shading_vert.levels # 연직단면을 그릴 때 쓴 구간을 그대로 가져옴
fig = plt.figure()
ax = fig.add_subplot()

shading_hori = ax.contourf(data['lon'], data['lat'], data['air'][0,2,:,:], levels=levels)
ax.plot(cross['lon'], cross['lat'], color='black')
ax.set_xlim(120, 135)
ax.set_ylim(30, 45)
cbar = fig.colorbar(shading_hori)

 

검은선을 보시면 위도 35도, 경도 125도부터 위도 40도, 경도 130도까지 선이 잘 그어져 있습니다.

colorbar도 완전히 동일하게 설정했으므로 연직단면의 850 hPa의 색과 일치합니다.