이전에 동아시아 기상도를 동적으로 그린 바 있다. 이번에는 아시아-유럽 기상도와 북반구 기상도를 작성해보자. 여기서는 500hPa 데이터를 예로 들어 설명한다.

데이터 출처 및 다운로드

ECMWF(유럽 기상 예보 센터)의 재분석 데이터인 ERA5를 사용하여 시각화를 진행한다. 온도(t(K)), 지형고도(z(m)) 데이터를 활용한다.

ERA5 재분석 데이터 다운로드 링크: ERA5 시간별 압력층 데이터(1940년부터 현재까지)(copernicus.eu)

아시아-유럽 기상도

아시아-유럽 지역의 등고선 및 등온선 분포를 탐구하기 위해 중국 지역을 확장하여 시각화한다. 2019년 11월 24일 00시 데이터를 예시로 사용한다.

다음 사이트에서 제공하는 shp 데이터로 행정경계를 그린다.

CN-border: 중국 국경 및 성 경계 데이터 - GMT 중국 매뉴얼(gmt-china.org)

코드 예시:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
from scipy.ndimage import gaussian_filter

# 라이브러리 불러오기
ds = xr.open_dataset(r"D:\ERA5\Z-T-U-V-2019.11.24-00h.nc")

# 데이터 선택 및 범위 설정
sel_time = '2019-11-24T00:00:00.000000000'
sel_lev = 500
geopot = ds.z.loc[sel_time, sel_lev, 80:-10, 0:360]
temp = ds.t.loc[sel_time, sel_lev, 80:-10, 0:360]
u_wind = ds.u.loc[sel_time, sel_lev, 80:-10, 0:360]
v_wind = ds.v.loc[sel_time, sel_lev, 80:-10, 0:360]

# 단위 변환
temp -= 273.15  # 섭씨 단위
lon, lat = temp.longitude, temp.latitude
geo_data = mpcalc.geopotential_to_height(geopot * units('m²/s²')) / 10

# 시각화 설정
fig = plt.figure(figsize=[12.5, 8], dpi=500)
projection = ccrs.LambertConformal(central_latitude=90, central_longitude=104)
ax = plt.axes(projection=projection)

# 배경 요소 추가
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='white', edgecolor='gray')
ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='#e7ffff')

# 제목 및 정보 추가
plt.title('UTC ' + sel_time[:-13] + ' ' + str(sel_lev) + 'hPa' + ' 온도($^o$C) 고도(dagpm)', loc='left', c='blue', fontsize=15)
plt.text(0.0, -0.03, 'ECMWF Reanalysis v5(ERA5)', transform=ax.transAxes, fontsize=12, color='blue', ha='left')

# 경계선 데이터 로드
with open(r"china-geospatial-data-GB2312\CN-border-La.gmt") as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

# 경계선 그리기
for line in borders:
    ax.plot(line[0::2], line[1::2], '-', color='#9b9f91', lw=1, transform=ccrs.Geodetic())

# 지도 범위 설정
ax.set_extent([15, 150, 5, 60], crs=ccrs.PlateCarree())

# 격자선 추가
gl = ax.gridlines(draw_labels=False, linewidth=0.5, color='k', alpha=0.6, linestyle='--')

# 좌표 라벨 설정
x_ticks = np.arange(0, 180+1, 10)
y_ticks = np.arange(0, 70+1, 10)
gl.xlocator = mticker.FixedLocator(x_ticks)
gl.ylocator = mticker.FixedLocator(y_ticks)

# 라벨 표시
for xtick in np.arange(0, 180+1, 30):
    ax.text(xtick, 42, f'{xtick}°E', horizontalalignment='center', transform=ccrs.PlateCarree(), color='k', fontsize=8)
for ytick in np.arange(0, 75+1, 15):
    ax.text(60, ytick, f'{ytick}°N', verticalalignment='center', transform=ccrs.PlateCarree(), color='k', fontsize=8)

# 필터 적용
smooth_geo = gaussian_filter(geo_data, sigma=2)
smooth_temp = gaussian_filter(temp, sigma=3)

# 등고선 및 등온선 그리기
contour = ax.contour(lon, lat, smooth_geo, levels=np.arange(500, 601, 4), colors="blue", linewidths=1, transform=ccrs.PlateCarree())
contour_temp = ax.contour(lon, lat, smooth_temp, levels=np.arange(-48, 5, 4), colors="red", linewidths=1, transform=ccrs.PlateCarree())

# 값 표시
ax.clabel(contour, fmt="%2.0f", fontsize=10, colors="blue")
ax.clabel(contour_temp, fmt="%2.0f", fontsize=10, colors="red")

# 특정 등고선 추가
ax.contour(lon, lat, smooth_geo, 25, levels=[588], linewidths=1.5, colors='blue', transform=ccrs.PlateCarree())

# 최대/최소값 표시
def mark_extremes(ax, lon, lat, data, ext_type, size, symbol, color, transform):
    if ext_type == 'max':
        filtered = maximum_filter(data, size, mode='nearest')
    elif ext_type == 'min':
        filtered = minimum_filter(data, size, mode='nearest')
    else:
        raise ValueError('Invalid type')
    m_xy, m_xx = np.where(filtered == data)
    for i in range(len(m_xy)):
        ax.text(lon[m_xy[i], m_xx[i]], lat[m_xy[i], m_xx[i]], symbol, color=color, size=12, transform=transform)

mark_extremes(ax, lon, lat, smooth_geo, 'max', 65, 'H', 'b', transform=ccrs.PlateCarree())
mark_extremes(ax, lon, lat, smooth_geo, 'min', 65, 'L', 'r', transform=ccrs.PlateCarree())
mark_extremes(ax, lon, lat, smooth_temp, 'max', 65, 'W', 'r', transform=ccrs.PlateCarree())
mark_extremes(ax, lon, lat, smooth_temp, 'min', 65, 'C', 'b', transform=ccrs.PlateCarree())

# 풍향선 그리기
u_wind, v_wind = np.array(u_wind), np.array(v_wind)
ax.barbs(lon[::32], lat[::32], u_wind[::32, ::32], v_wind[::32, ::32], barbcolor='k', lw=0.5, length=5, 
         barb_increments=dict(half=2, full=4, flag=20), transform=ccrs.PlateCarree())

# 파일 저장
plt.subplots_adjust(left=0.01, right=0.99, top=0.995, bottom=0.005)
fig.savefig(f'아시아-유럽 기상도-{sel_time[:-16]}-500hpa.png')

결과:

북반구 기상도

북반구의 대규모 고기압/저기압 분포를 시각화하기 위해 극지방 투영법을 사용한다. 2024년 10월 25일 UTC 00시(베이징 시간 오전 8시) 500hPa 데이터를 예로 들어 설명한다.

코드 예시:

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
from scipy.ndimage import gaussian_filter, maximum_filter, minimum_filter
from cartopy.util import add_cyclic_point

# 최대/최소값 표시 함수
def mark_extremes(ax, lon, lat, data, ext_type, size, symbol, color, transform):
    if ext_type == 'max':
        filtered = maximum_filter(data, size, mode='nearest')
    elif ext_type == 'min':
        filtered = minimum_filter(data, size, mode='nearest')
    else:
        raise ValueError('Invalid type')
    m_xy, m_xx = np.where(filtered == data)
    for i in range(len(m_xy)):
        ax.text(lon[m_xy[i], m_xx[i]], lat[m_xy[i], m_xx[i]], symbol, color=color, size=16, transform=transform)

# 데이터 로드
ds = xr.open_dataset(r"D:\ERA5\Z-T-U-V-2024.10.25.26-00.12h.nc")
sel_time = '2024-10-25T00:00:00.000000000'
sel_lev = 500
geopot = ds.z.loc[sel_time, sel_lev, 90:-20, :]
temp = ds.t.loc[sel_time, sel_lev, 90:-20, :]

# 단위 변환
temp -= 273.15
lon, lat = temp.longitude, temp.latitude
geo_data = mpcalc.geopotential_to_height(geopot * units('m²/s²')) / 10

# 시각화 설정
fig = plt.figure(figsize=[8, 8], dpi=500)
projection = ccrs.NorthPolarStereo(central_longitude=110)
ax = plt.axes(projection=projection)

# 배경 요소 추가
ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor='#bdd2ff')
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='#ffffff')

# 제목 및 정보 추가
plt.title('UTC ' + sel_time[:-13] + ' ' + str(sel_lev) + 'hPa' + ' 온도($^o$C) 고도(dagpm)', loc='left', c='blue', fontsize=15)
plt.text(0.0, -0.03, 'ECMWF Reanalysis v5(ERA5)', transform=ax.transAxes, fontsize=12, color='blue', ha='left')

# 경계선 데이터 로드
with open(r"china-geospatial-data-GB2312\CN-border-La.gmt") as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

# 경계선 그리기
for line in borders:
    ax.plot(line[0::2], line[1::2], '-', color='#9b9f91', lw=1, transform=ccrs.Geodetic())

# 지도 범위 설정
ax.set_extent([0, 360, 0, 90], crs=ccrs.PlateCarree())

# 격자선 추가
gl = ax.gridlines(draw_labels=False, linewidth=0.5, color='k', alpha=0.6, linestyle='--')

# 좌표 라벨 설정
x_ticks = np.arange(-180, 160+1, 20)
y_ticks = np.arange(0, 30+1, 10)
gl.xlocator = mticker.FixedLocator(x_ticks)
gl.ylocator = mticker.FixedLocator(y_ticks)

# 라벨 표시
for xtick in x_ticks:
    ax.text(xtick, 15-3, f'{xtick}°E', horizontalalignment='center', transform=ccrs.PlateCarree(), color='k', fontsize=10)
for ytick in y_ticks:
    ax.text(40-6, ytick, f'{ytick}°N', verticalalignment='center', transform=ccrs.PlateCarree(), color='k', fontsize=10)

# 필터 적용
smooth_geo = gaussian_filter(geo_data, sigma=6)
smooth_temp = gaussian_filter(temp, sigma=6)

# 최대/최소값 표시
mark_extremes(ax, lon, lat, smooth_geo, 'max', 120, 'H', 'b', transform=ccrs.PlateCarree())
mark_extremes(ax, lon, lat, smooth_geo, 'min', 120, 'L', 'r', transform=ccrs.PlateCarree())
mark_extremes(ax, lon, lat, smooth_temp, 'max', 120, 'W', 'r', transform=ccrs.PlateCarree())
mark_extremes(ax, lon, lat, smooth_temp, 'min', 120, 'C', 'b', transform=ccrs.PlateCarree())

# 데이터 주기성 처리
smooth_geo, cycle_lon = add_cyclic_point(smooth_geo, coord=lon)
smooth_temp, cycle_lon = add_cyclic_point(smooth_temp, coord=lon)
cycle_lon, cycle_lat = np.meshgrid(cycle_lon, lat)

# 등고선 및 등온선 그리기
contour = ax.contour(cycle_lon, cycle_lat, smooth_geo, levels=np.arange(500, 601, 4), colors="blue", linewidths=1, transform=ccrs.PlateCarree())
contour_temp = ax.contour(cycle_lon, cycle_lat, smooth_temp, levels=np.arange(-60, 1, 4), colors="red", linewidths=1, transform=ccrs.PlateCarree())

# 값 표시
ax.clabel(contour, fmt="%2.0f", fontsize=10, colors="blue")
ax.clabel(contour_temp, fmt="%2.0f", fontsize=10, colors="red")

# 특정 등고선 추가
ax.contour(cycle_lon, cycle_lat, smooth_geo, 25, levels=[588], linewidths=1.5, colors='blue', transform=ccrs.PlateCarree())
ax.clabel(ax.contour(cycle_lon, cycle_lat, smooth_geo, 25, levels=[588], colors='blue'), fmt="%2.0f", fontsize=10, colors="blue")

# 파일 저장
plt.subplots_adjust(left=0.0, right=1, top=0.9, bottom=0.1)
fig.savefig(f'북반구 기상도-{sel_time[:-16]}-500hpa.png')

결과:

겨울에 가까워짐에 따라 여러 대규모 고기압/저기압 구조를 확인할 수 있다. 예를 들어 동아시아 해안의 고기압, 알라스카 저기압, 북미 동부 고기압, 서유럽 해안 저기압, 서유럽 동부 고기압 등이 나타난다.