[Python/Cartopy]모델자료 활용한 일기도 묘화

주제 :  모델자료 활용한 일기도 묘화 작업

데이터: 국지예보모델(LDAPS)
              기상자료개방포털 ▶ 데이터  ▶수치모델  ▶ 국지예보모델(LDAPS)
             data.kma.go.kr/data/rmt/rmtList.do?code=340&pgmNo=65

버  젼 :  conda 4.9.0 / Python 3.8.3 / numpy 1.19.2 / pandas 1.1.3
              - cartopy 설치 :  conda install -c conda-forge cartopy
              - pygrib 설치 : conda install pygrib 
              ※ conda 및 pip를 병행하여 사용하는 경우, 충돌 가능성이 있음으로 개발자는 하나로 통일하여 사용하기를 권장

---

[Python/Basemap]모델자료를 활용한 일기도 묘화

이전 포스팅은 Basemap을 활용하였다면, 이번에는 Cartopy를 사용하여 일기도 묘화를 하고자 함
더 이상 Basemap은 업데이트가 진행되지 않으며, 파이선 버젼 업됨에 따라 호환성 이슈 발생
이에 Basemap보다는 Cartopy를 활용하는 사례가 많아짐

이전 포스팅에서 설명한 내용을 확인하려면 "더보기"를 클릭

# 데이터 다운로드 및 설명

기상자료개방포털 > 데이터 > 수치모델 > 지역예보모델

일반적으로 접하는 excel, csv 파일에 비해, 수치모델자료는 크기가 큼
: 수치모델자료는 기상자료개방포털에서 다운로드할 수 있으며, 그 방법은 자신이 원하는 데이터를 조회하여 담은 뒤(①~③), 마이페이지에서 자료를 신청하면, E-mail이나 SNS로 FTP계정의 일시적인 ID/PW를 부여 받음
  이 FTP 서버를 활용하여 데이터 다운로드 가능(④)함

수치예보모델 종류 (출처: 기상청 수치모델관리실)

국지예보모델 이해(출처: 국가기상슈퍼컴퓨터센터)

## 지역예보모델(RDAPS)

지역예보모델(UM 12kmL70)의 공간 해상도는 12km이며, 연직으로 약 80km까지 70층으로 구성되며, 3시간 간격으로 전지구예보모델로부터 경계장을 제공받아 1일 4회(00, 06, 12, 18UTC) 87시간 예측을 수행함

지역예보모델은 4차원 변분자료 동화 기법을 이용하여 6시간 간격의 분석-예측 순환 체계로 운영되고 있으며, 지역예보모델의 산출자료는 등압면, 모델면, 단일면 자료 3종류가 제공되고, 데이터 형식은 WMO에서 제시한 GRIB2 형식으로 제공함

## 국지예보모델(LDAPS)

국지예보모델(UM 1.5kmL70)의 공간 해상도는 1.5km이며, 연직으로 약 40km까지 70층으로 구성되며, 3시간 간격으로 전지구모델로부터 경계장을 제공받아 1일 8회(00, 06, 12, 18UTC : 36시간 예측, 03, 09, 15, 18UTC : 3시간 예측) 수행함

국지예보모델은 3차원 변분자료 동화 기법을 이용하여 각각의 자체 분석 - 예측 순환 체계로 운영하고 있으며, 국지예보모델의 산출자료는 등압면, 모델면, 단일면 자료 3종류가 제공되며, 데이터 형식은 WMO에서 제시한 GRIB2 형식으로 제공함

수치예보모델의 종류 및 비교자료(출처: 국가기상슈퍼컴퓨터센터)

# GRIB1/GRIB2 

GRIB는 수치모델에서 산출되는 격자자료를 WMO의 표준규격에 따라 Binary 형식으로 구성된 데이터로,
데이터의 정보와 그 자체를 동시에 압축 저장이 가능하여 저장과 시간 측면에서 효과적임
(기상) 대기/해양에서 많이 사용함(netCDF(. nc) 파일도 많이 사용함)

## GRIB1/GRIB2 데이터 읽기

grib파일을 가장 쉽게 읽는 방법은 pygrib 라이브러리 활용 (xarray 또한 가능)

import pygrib

FILE = './Data/l015v070erlounish000.2020102900.gb2'
grbs = pygrib.open(FILE)


for g in grbs :                                   # 해당 데이터의 어떤 변수들이 있는 확인가능
    print(g)

 

"print(g)" 출력결과 일부

모델자료 중 일기도 묘화를 위해 Pressure reduced to MSL 활용
단위는 Pa이나, (지상) 일기도에서는 hPa 사용함으로 (÷) 100 하여 단위 환산 적용 (1 hPa = 100 Pa)

grb = grbs.select(name= "Pressure reduced to MSL")[0]
data = grb.values/100                                 # 입력자료 단위는 Pa임으로 일기도묘화를 위해 hPa로 변환
lat, lon = grb.latlons()

1. 일기도 묘화

Cartopy 활용한 일기도 묘화작업

## 라이브러리 불러오기
import pygrib
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import scipy.ndimage as ndimage


## 필요한 데이터 읽기
FILE = './Data/l015v070erlounish000.2020102900.gb2'
grbs = pygrib.open(FILE)
grb = grbs.select(name= "Pressure reduced to MSL")[0]

data = grb.values/100 
data = ndimage.gaussian_filter(data, sigma=5, order=0)
lat, lon = grb.latlons()


ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()                                  # 지도의 배경그림을 설정
ax.set_extent([122,131.1,43,33])
ax.coastlines('50m')

cs =ax.contour(lon, lat, data,levels = np. arange(1000.,1250.,2.), colors ="black")
ax.clabel(cs, inline=True, fontsize=8)
ax.gridlines(crs =ccrs.PlateCarree(), draw_labels= True)

plt.show()

---

Reference

• Mapping with cartopy in python  semba-blog.netlify.app/07/04/2020/mapping-with-cartopy-in-python/