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

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

데이터 : 국지예보모델(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
              - basemap 설치 : conda install basemap
              - pygrib 설치 : conda install pygrib 
              - conda 업데이트 주의 (conda update -n base conda) 
                 * (2020.11.01 기준) 업데이트 시, basemap 지원 문제로 충돌 발생
              ※ conda 및 pip를 병행하여 사용하는 경우, 충돌 가능성이 있음으로 개발자는 하나로 통일하여 사용하기를 권장

---

1. 데이터 다운로드 및 설명

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

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

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

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

1.1. 지역예보모델(RDAPS)

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

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

1.2. 국지예보모델(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 형식으로 제공함

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

2. GRIB1/GRIB2 

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

2.1. 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()

3. 일기도 묘화

pcolormesh 활용

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

m = Basemap(projection='mill', lat_ts=10,
               llcrnrlon=lon.min(),
               urcrnrlon=lon.max(),
               llcrnrlat=lat.min(),
               urcrnrlat=lat.max(),
               resolution=None)

x, y = m(lon,lat)
cs = m.pcolormesh(x,y,data,shading='flat',cmap=plt.cm.jet)

m.drawparallels(np.arange(-90.,80.,5.),labels=[1,0,0,0])
m.drawmeridians(np.arange(-180.,180.,5.),labels=[0,0,0,1])
plt.colorbar(cs,orientation='vertical')
plt.show()

3.1. contour 활용한 등압선 그리기

 

contour 활용하여 등압선 표시

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

plt.figure(figsize=(10,10))
m = Basemap(projection='merc', lat_ts=10,
               llcrnrlon=lon.min(),
               urcrnrlon=lon.max(),
               llcrnrlat=lat.min(),
               urcrnrlat=lat.max(),
               resolution='h')

m.drawcoastlines()
m.drawcountries()
m.drawmapboundary()

parallels = np.arange(30.,45.,2.)
m.drawparallels(parallels,labels=[True,False,False,False])
meridians = np.arange(120.,150.,5.)
m.drawmeridians(meridians,labels=[True,False,False,True])

x, y = m(lon,lat)
cs = m.contour(x,y,data,shading='flat',levels = np. arange(1000,1250,2), colors= "blue")
plt.colorbar(cs,orientation='vertical')
plt.show()

 

3.2. 등압선은 Smooth 하게 그리기

import scipy.ndimage as ndimage

data = ndimage.gaussian_filter(data, sigma=5, order=0)

3.3. (최종) 모델자료를 활용하여 일기도 묘화 완성

모델자료를 활용하여 일기도 묘화 완성

# 라이브러리
import pygrib
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as ndimage
from mpl_toolkits.basemap import Basemap

# 수치모델자료 Read
FILE = './Data/l015v070erlounish000.2020102900.gb2'
grbs = pygrib.open(FILE)
grb = grbs.select(name= "Pressure reduced to MSL")[0]

## 수치모델자료 변경
dat = grb.values/100                                    # 단위변환(Pa → hPa)
data = ndimage.gaussian_filter(dat, sigma=5, order=0)   # 가우시안 필터로 smooth하게 데이터 처리
lat, lon = grb.latlons()


# Basemap로 한반도 그리기
plt.figure(figsize=(10,10))
m = Basemap(projection='merc', lat_ts=10,
               llcrnrlon=lon.min(),
               urcrnrlon=lon.max(),
               llcrnrlat=lat.min(),
               urcrnrlat=lat.max(),
               resolution='h')

m.drawcoastlines()
m.drawcountries()
m.drawmapboundary()

## 위경도표시
parallels = np.arange(30.,45.,2.)
m.drawparallels(parallels,labels=[True,False,False,False])        #True=1, False =0 으로 표현가능
meridians = np.arange(120.,150.,5.)
m.drawmeridians(meridians,labels=[True,False,False,True])

## 등압선 그리기
x, y = m(lon,lat)
cs = m.contour(x,y,data,shading='flat',
               levels = np. arange(1000,1250,2), colors= "blue")  # 한반도데이터로 압력변화가 작게 나타나 
                                                                  # 2간격으로 표시
plt.clabel(cs, inline=True, fontsize=8)                           # 등압선 압력표시
plt.show()

---

Reference

• 국가기상슈퍼컴퓨터센터 www.kma.go.kr/aboutkma/intro/supercom/model/model_category.jsp
• 기상산업 활성화를 위한 수치예보모델자료 활용 가이드(기상청 수치모델관리실) 발간등록번호 11-1360395-000252-01
• Density and Contour Plots jakevdp.github.io/PythonDataScienceHandbook/04.04-density-and-contour-plots.html
• Basemap tutorial basemaptutorial.readthedocs.io/en/latest/basemap.html#passing-the-bounding-box