반응형
정보
-
업무명 : 파이썬 Grib 및 Grb2 형식인 수치예측모델 (ECMWF) 자료를 이용한 가시화
-
작성자 : 이상호
-
작성일 : 2019-10-22
-
설 명 :
-
수정이력 :
-
2020-02-08 : 소스 코드 명세 추가
-
내용
[개요]
-
안녕하세요? 웹 개발 및 연구 개발을 담당하고 있는 해솔입니다.
-
현대 기상 예측에서 수치예보는 40 %로서 큰 비중을 차지합니다. 기상청에 따르면 "수치 예보는 물리학의 방정식에 의해 바람과 기온 등의 기온 변화를 컴퓨터로 계산하여 미래의 대기 상태를 예측하는 방법"이라고 합니다.
-
즉 컴퓨터를 통해 지구의 상태를 시뮬레이션하여 기상 예측하는 기술입니다. 이러한 시뮬레이션 방법은 여러 종류가 있으나 결과는 "GPV (Grid Point Value)" 형식으로 제공하는 것이 대부분입니다.
-
기상 예측은 시간과 격자에 대한 여러 기상 정보를 내포하고 "Grib" 또는 "Grb2" 형식으로 배포하고 있습니다.
-
따라서 대기과학에서 사용되는 Grib 및 Grb2 설명뿐만 아니라 Python을 이용하여 자료 전처리 및 가시화를 소개해 드리고자 합니다.
[특징]
-
Grib 및 Grb2 형태인 수치예측모델 자료를 이해하기 위해 가시화 도구가 필요하며 이 프로그램은 이러한 목적을 달성하기 위해 고안된 소프트웨어
[기능]
-
수치예측모델 자료를 이용한 가시화
-
새로운 규칙 격자에 대한 공간 내삽
[활용 자료]
-
모델명 : ECMWF
-
자료종류 : 2m 대기온도
-
영역 : 전지구
-
해상도 : 25 km
-
확장자 : .grib 및 .grb2
-
기간 : 2014년 01월
[자료 처리 방안 및 활용 분석 기법]
-
원도우 (Windos 10)에서 최소한의 노력으로 "Grib" 및 "Grb2"을 가시화하기 위해서 "NetCDF"로 변환하여 사용
-
그에 따라 wgrib, wgrib2, cdo와 같은 배치 파일에 대해 "시스템 환경변수 (PATH)" 설정 필연
-
해당 "배치 파일에 대한 설치" 포스팅 참조
[Python] 파이썬 원도우 (Window 10)에서 수치예측모델 (Grib, Grb2) 자료 처리를 위한 wgrib, wgrib2, cdo 설치 방법
정보 업무명 : 파이썬 원도우 (Window 10)에서 수치예측모델 (Grib, Grb2) 자료 처리를 위한 wgrib, wgrib2, cdo 설치 방법 작성자 : 이상호 작성일 : 2020-02-09 설 명 : 수정이력 : 내용 [개요] 안녕하세요? 웹..
shlee1990.tistory.com
-
"Grib"을 "NetCDF"로 변환
os.system('cdo -f nc copy ERA_interim-sfc-201606.grib grib.nc')
-
"Grb2"을 "NetCDF"로 변환
os.system('wgrib2 20140101_0600_003.grb2 -netcdf test2.nc')
[사용법]
-
입력 자료를 동일 디렉터리 위치
-
소스 코드를 실행 (Jupyter notebook 실행)
-
가시화 결과를 확인
[사용 OS]
-
Window 10
[사용 언어]
-
Python v2.7
[사용 배치 프로그램]
-
wgrib v1.7.3.1
-
wgrib2 v0.1.9.9.9
-
cdo v1.6.4
소스 코드
[명세]
-
라이브러리 읽기
# Library
import pandas as pd
import numpy as np
import sys
# import iris
import os
import matplotlib.pyplot as plt
from dplython import *
# (DplyFrame, X, diamonds, select, sift, sample_n,
# sample_frac, head, arrange, mutate, group_by, summarize, DelayFunction)
from scipy.stats import linregress
from matplotlib import pyplot as plt
from IPython.display import Image
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import Normalize
import matplotlib
import matplotlib.cm as cm
import seaborn as sns
from scipy.stats import linregress
from matplotlib import rcParams
from netCDF4 import Dataset
import struct
import binascii
from mpl_toolkits.basemap import addcyclic
from netCDF4 import num2date, date2num, date2index
import datetime
from pyhdf.SD import SD, SDC
import h5py
-
NetCDF 파일 읽기
-
위도, 경도, 시간, 변수 (u10, v10, t2m)
-
ncfile_path = 'ECMWF/era-i_sfc_2014mon.nc'
ncfile = Dataset(ncfile_path)
print(ncfile)
print(ncfile.variables.keys())
<type 'netCDF4._netCDF4.Dataset'>
root group (NETCDF3_64BIT_OFFSET data model, file format NETCDF3):
Conventions: CF-1.6
history: 2016-01-25 07:44:11 GMT by grib_to_netcdf-1.14.3: grib_to_netcdf /data/data01/scratch/_mars-atls12-95e2cf679cd58ee9b4db4dd119a05a8d-O7rdSz.grib -o /data/data01/scratch/_grib2netcdf-atls07-95e2cf679cd58ee9b4db4dd119a05a8d-Zw1tQZ.nc -utime
dimensions(sizes): longitude(480), latitude(241), time(12)
variables(dimensions): float32 longitude(longitude), float32 latitude(latitude), int32 time(time), int16 u10(time,latitude,longitude), int16 v10(time,latitude,longitude), int16 t2m(time,latitude,longitude)
groups:
[u'longitude', u'latitude', u'time', u'u10', u'v10', u't2m']
-
NetCDF 헤더 보기
-
원도우 (Window)에서 ncdump를 사용할 경우 Anaconda Cloud에서 "conda install -c anaconda netcdf4" 설치 필요
-
ncdump 외에 "ncBrowse", "Panoply" 등 다양한 도구가 있으니 참고 바램
-
shlee1990.tistory.com
!ncdump -h ECMWF/era-i_sfc_2014mon.nc
netcdf ECMWF/era-i_sfc_2014mon {
dimensions:
longitude = 480 ;
latitude = 241 ;
time = UNLIMITED ; // (12 currently)
variables:
float longitude(longitude) ;
longitude:units = "degrees_east" ;
longitude:long_name = "longitude" ;
float latitude(latitude) ;
latitude:units = "degrees_north" ;
latitude:long_name = "latitude" ;
int time(time) ;
time:units = "hours since 1900-01-01 00:00:0.0" ;
time:long_name = "time" ;
time:calendar = "gregorian" ;
short u10(time, latitude, longitude) ;
u10:scale_factor = 0.00040673147983914 ;
u10:add_offset = -0.840284103807546 ;
u10:_FillValue = -32767s ;
u10:missing_value = -32767s ;
u10:units = "m s**-1" ;
u10:long_name = "10 metre U wind component" ;
short v10(time, latitude, longitude) ;
v10:scale_factor = 0.000413141076998056 ;
v10:add_offset = 0.9845084143553 ;
v10:_FillValue = -32767s ;
v10:missing_value = -32767s ;
v10:units = "m s**-1" ;
v10:long_name = "10 metre V wind component" ;
short t2m(time, latitude, longitude) ;
t2m:scale_factor = 0.00168685886309188 ;
t2m:add_offset = 258.497920608654 ;
t2m:_FillValue = -32767s ;
t2m:missing_value = -32767s ;
t2m:units = "K" ;
t2m:long_name = "2 metre temperature" ;
// global attributes:
:Conventions = "CF-1.6" ;
:history = "2016-01-25 07:44:11 GMT by grib_to_netcdf-1.14.3: grib_to_netcdf /data/data01/scratch/_mars-atls12-95e2cf679cd58ee9b4db4dd119a05a8d-O7rdSz.grib -o /data/data01/scratch/_grib2netcdf-atls07-95e2cf679cd58ee9b4db4dd119a05a8d-Zw1tQZ.nc -utime" ;
}
-
NetCDF 파일에서 세부 정보 보기
-
위도, 경도에 대해 단일 배열을 지님
-
변수 (u10, v10, t2m)에 대해 2차원 배열을 지님
-
lon = ncfile.variables[u'longitude'][:]
lat = ncfile.variables[u'latitude'][:]
# time = ncfile.variables[u'time'][:]
u10 = ncfile.variables[u'u10'][:]
v10 = ncfile.variables[u'v10'][:]
t2m = ncfile.variables[u't2m'][:]
# print(lat, lon)
print(lat.shape, lon.shape)
print(max(lat), min(lat))
print(max(lon), min(lon))
((241L,), (480L,))
(90.0, -90.0)
(359.25, 0.0)
-
시간에 대해 1차원 변환
-
dplython를 이용하기 위해서 1차원 필요
-
시간은 단일 12개 배열로서 위/경도 격자 (241 * 480개)에 대해 1차원 생성
-
time = ncfile.variables[u'time']
print(times)
dtime = num2date(time[:], time.units)
times = pd.to_datetime(dtime)
time_fn = times.strftime("%B %Y")
times1D = np.repeat(times, 241*480)
print(times1D.shape)
DatetimeIndex(['2014-01-01', '2014-02-01', '2014-03-01', '2014-04-01',
'2014-05-01', '2014-06-01', '2014-07-01', '2014-08-01',
'2014-09-01', '2014-10-01', '2014-11-01', '2014-12-01'],
dtype='datetime64[ns]', freq=None)
(1388160L,)
-
위/경도, 변수에 대해 1차원 변환
-
dplython를 이용하기 위해서 1차원 필요
-
2차원 변수 (t2m)를 1차원으로 변환
-
lon2D, lat2D = np.meshgrid(lon, lat)
lon1D = np.reshape(lon2D, (1,np.product(lon2D.shape)))[0]
lat1D = np.reshape(lat2D, (1,np.product(lat2D.shape)))[0]
# val = [1,2,3,4]
# np.tile(val,2) = [1,2,3,4,1,2,3,4]
# np.repeat(val,2) = [1,1,2,2,3,3,4,4]
lats1D = np.tile(lat1D, 12)
lons1D = np.tile(lon1D, 12)
t2m1D = np.reshape( t2m, ( 1, np.product(t2m.shape) ) )[0]
print(lons1D.shape)
print(t2m1D.shape)
(1388160L,)
(1388160L,)
-
Data Frame 설정
-
1차원 시간, 위도, 경도, 변수 (t2m)
-
data = pd.DataFrame( np.column_stack( [times1D.year, times1D.month, times1D.day,
times1D.hour, times1D.minute, times1D.second, times1D.microsecond,
lons1D, lats1D, t2m1D] ),
columns=['year', 'month', 'day', 'hour', 'minute', 'second', 'microsecond',
'lon', 'lat', 't2m'], )
data.head()
# data.tail()
year month day hour minute second microsecond lon lat t2m
0 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 0.00 90.0 252.401613
1 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 0.75 90.0 252.401613
2 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 1.50 90.0 252.401613
3 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 2.25 90.0 252.401613
4 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 3.00 90.0 252.401613
-
Data Frame를 통해 L1 전처리
-
최근 데이터 분석에서 사용되는 R의 "dplyr" 라이브러리와 유사한 "dplython"를 사용
-
시간 (2014년 1월)에 해당되는 자료 선택
-
data_L1 = ( DplyFrame(data) >>
sift( X.year == 2014) >>
sift( X.month == 1)
)
data_L1.head()
year month day hour minute second microsecond lon lat t2m
0 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 0.00 90.0 252.401613
1 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 0.75 90.0 252.401613
2 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 1.50 90.0 252.401613
3 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 2.25 90.0 252.401613
4 2014.0 1.0 1.0 0.0 0.0 0.0 0.0 3.00 90.0 252.401613
- 가시화를 위한 설정
- 초기 설정
- 크기 : figure
- 스타일 : style
- 폰트 : rc, rcParams
- 지도 해상도 : Basemap
- 변수 설정
- 위도 (lon), 경도 (lat), 2m 대기 온도 (t2m)
- 그림 설정
- 산점도 : pcolormesh
- 컬러바 : colorbar
- 해안선 : drawcostlines, drawcountries, drawmapboundary
- 수평/수직 그리드 : drawparallels, drawmeridians
- 그림 제목 : title
%matplotlib inline
# define plot size in inches (width, height) & resolution(DPI)
plt.figure(figsize=(12, 10))
# style
plt.style.use('seaborn-darkgrid')
# define font size
plt.rc("font", size=22)
rcParams['font.family'] = 'New Century Schoolbook'
# plt.rcParams["font.weight"] = "bold"
# plt.rcParams['axes.labelsize'] = 26
# plt.rcParams['xtick.labelsize'] = 26
# plt.rcParams['ytick.labelsize'] = 26
# plt.rcParams["axes.labelweight"] = "bold"
# plt.rcParams["axes.titleweight"] = "bold"
m = Basemap(projection='cyl', lon_0 = 180, lat_0 = 0)
# m = Basemap(projection='mill', lon_0 = 180, lat_0 = 0)
# X = data_L1["lon"].values
# Y = data_L1["lat"].values
# VAL = data_L1["t2m"].values
# m.scatter(X, Y, c=VAL, s=15, marker="s", zorder=1, vmin=220, vmax=320, cmap=plt.cm.get_cmap('jet'), alpha=1.0)
X, Y = m(*np.meshgrid(lon, lat))
VAL = t2m[0][:][:]
m.pcolormesh(X, Y, VAL, shading='flat', vmin=220, vmax=320, cmap=plt.cm.get_cmap('jet'), alpha=1.0)
# m.colorbar(location='right',label='2 Meter Temperature [K]')
m.colorbar(location='bottom',label='2 Meter Temperature [K]', pad=0.5)
m.drawcoastlines()
m.drawcountries()
m.drawmapboundary(fill_color='white')
m.drawparallels(np.arange(-90,91,30), labels=[1,0,0,0], dashes=[2,2])
m.drawmeridians(np.arange(0,361,60), labels=[0,0,0,1], dashes=[2,2])
plt.title('ECMWF / ERA : %04d-%02d \n' %( data_L1.year[0].astype('int'), data_L1.month[0].astype('int') ) )
# plt.savefig("FIG/Scane_analysis_Himawari-8_AHI_RSR.png")
plt.show()
-
일반적으로 ECMWF 자료를 이용하여 고/저 해상도로 격자화하기 위해서 규칙 및 불규칙 격자에 대해 공간 내삽 필연
-
라이브러리 읽기
from scipy.interpolate import griddata
import numpy as np
-
신규 격자 생성
-
기존 격자 (241 * 480개)에 대해 신규 격자 (150 * 150개) 생성
-
# 기존 격자
lon2D, lat2D = np.meshgrid(lon, lat)
val2D = t2m[0][:][:]
val1D = np.reshape(val2D, (1,np.product(val2D.shape)))[0]
lon1D = np.reshape(lon2D, (1,np.product(lon2D.shape)))[0]
lat1D = np.reshape(lat2D, (1,np.product(lat2D.shape)))[0]
print(val1D.shape)
print(lon1D.shape)
print(lat1D.shape)
# 새로운 격자 생성
nlon = np.linspace(0,359.25,150)
nlat = np.linspace(90,-90,150)
nlon2D, nlat2D = np.meshgrid(nlon, nlat)
print(nlon2D.shape)
(115680L,)
(115680L,)
(115680L,)
(150L, 150L)
-
Data Frame를 통해 L1 전처리
-
최근 공간 내삽에서 사용되는 "griddata" 라이브러리를 통해 "nearest"와 "linear" 및 "cubic"를 사용
-
공간 내삽할 경우 기존 격자 (1차원 위/경도, 변수) 및 신규 격자 (2차원 위/경도) 필요
-
신규 격자에 대한 2차원을 1차원으로 변환
-
from scipy.interpolate import griddata
# nval2D = griddata( (lon1D,lat1D), val1D, (nlon2D, nlat2D), method='nearest' )
nval2D = griddata( (lon1D,lat1D), val1D, (nlon2D, nlat2D), method='linear' )
# nval2D = griddata( (lon1D,lat1D), val1D, (nlon2D, nlat2D), method='cubic' )
print(nval2D.shape)
nval1D = np.reshape(nval2D, (1,np.product(nval2D.shape)))[0]
nlat1D = np.reshape(nlat2D, (1,np.product(nlat2D.shape)))[0]
nlon1D = np.reshape(nlon2D, (1,np.product(nlon2D.shape)))[0]
(150L, 150L)
-
공간 내삽에 대한 가시화
%matplotlib inline
# define plot size in inches (width, height) & resolution(DPI)
plt.figure(figsize=(12, 10))
# style
plt.style.use('seaborn-darkgrid')
# define font size
plt.rc("font", size=22)
rcParams['font.family'] = 'New Century Schoolbook'
# plt.rcParams["font.weight"] = "bold"
# plt.rcParams['axes.labelsize'] = 26
# plt.rcParams['xtick.labelsize'] = 26
# plt.rcParams['ytick.labelsize'] = 26
# plt.rcParams["axes.labelweight"] = "bold"
# plt.rcParams["axes.titleweight"] = "bold"
m = Basemap(lon_0 = 180, lat_0 = 0, projection='cyl', resolution='c')
# X, Y = m(*np.meshgrid(nlon, nlat))
# VAL = nval2D
# m.pcolormesh(X, Y, VAL, shading='flat', vmin=220, vmax=320, cmap=plt.cm.get_cmap('jet'), alpha=1.0)
X, Y = m(nlon1D, nlat1D)
VAL = nval1D
m.scatter(X, Y, c=VAL, s=15, marker="s", zorder=1, vmin=220, vmax=320, cmap=plt.cm.get_cmap('jet'), alpha=1.0)
# m.colorbar(location='right',label='2 Meter Temperature [K]')
m.colorbar(location='bottom',label='2 Meter Temperature [K]', pad=0.5)
m.drawcoastlines()
m.drawcountries()
m.drawparallels(np.arange(-150.,120.,30.),labels=[1,0,0,0],dashes=[2,2])
m.drawmeridians(np.arange(-180.,180.,60.),labels=[0,0,0,1],dashes=[2,2])
plt.title('ECMWF / ERA : %04d-%02d \n' %( data_L1.year[0].astype('int'), data_L1.month[0].astype('int') ) )
# plt.savefig("FIG/Scane_analysis_Himawari-8_AHI_RSR.png")
plt.show()
[전체]
참고 문헌
[논문]
- 없음
[보고서]
- 없음
[URL]
- 없음
문의사항
[기상학/프로그래밍 언어]
- sangho.lee.1990@gmail.com
[해양학/천문학/빅데이터]
- saimang0804@gmail.com
반응형
'프로그래밍 언어 > Python' 카테고리의 다른 글
| [Python] 파이썬 NetCDF 형식인 Aqua/CERES 기상위성 자료를 이용한 가시화 (1) | 2019.10.23 |
|---|---|
| [Python] 파이썬 HDF 형식인 천리안위성 1A호 (COMS/MI) 기상위성 자료를 이용한 가시화 (4) | 2019.10.23 |
| [Python] 파이썬 바이너리 형식인 히마와리 8호 (Himawari-8/AHI) 기상 위성 자료를 이용한 가시화 (0) | 2019.10.22 |
| [Python] 파이썬 라디오 미터 및 라디오 존데 관측 자료를 이용한 고도별 온/습도 등고선 가시화 (0) | 2019.09.07 |
| [Python] 파이썬 라디오미터 밝기온도 자료를 이용한 시계열 그래프 (0) | 2019.09.07 |
최근댓글