최근 xarray라는 파이썬 패키지를 활용하여 다차원 시공간 데이터를 처리하는 사례들을 볼 수 있어서 살펴보았다.
XARRAY
왜 xarray를 만들게 되었고, 필요한지 등에 대해서는 홈페이지나 깃헙에 나와 있다. DB 또는 빅데이터/ NoSQL과 결합하여 좀 더 원활하게 Read/Write할 수 있게 된다면 훨씬 더 범용적으로 많이 쓰일 것 같다는 생각이 든다.
xarray의 데이터 구조
위의 docs.xarray.dev의 data-structures.html에서 설명되고 있는데, 그림으로 잘 설명해주고 있다.
netCDF와 Grib
xarray 에서 다루는 대표적인 다차원 배열 구조의 데이터 포맷이 netCDF와 Grib이다.
NetCDF(Network Common Data Form)
시공간 대기(기상) 및 해양 관측 정보를 저장할 수 있는 대표적인 데이터 포맷으로 자세한 설명은 아래의 블로그에 잘 나와 있다. ArcGIS는 물론 QGIS에서도 데이터를 로딩하고 시각화할 수 있다. 일반적으로 알려진 확장자는 .nc이고 xarray에서 처리할 수 있다.
Grib(GRIdded Binary)
netCDF와 함께 대기(기상) 분야에서 많이 쓰이는 포맷으로 WMO(국제기상기구)의 표준 포맷이다. 현재 최신 버전은 Grib2로 확장자는 .grib/ .grib2 또는 .gb2가 쓰이는 것 같다. grib에 대한 자세한 설명 역시 위에서 참조한 “이상우의 IDL 블로그” 에 잘 나와 있다. 이 분야의 전문가 분인 것 같다.
이제 다차원 시공간 데이터 포맷도 알아보고, 이를 다룰 수 있게 해주는 패키지도 알았으니 실제 데이터를 살펴보자.
ECMWF (European Centre for Medium-Range Weather Forecasts)
우리말로 하면 유럽중기예보센터라고 할 수 있겠다. 찾아보면 다른 사이트도 있겠지만 많이 검색되기도 하고, 있어 보이기도 하니 한 번 찾아보자.
과거 이력 데이터와 예측모델을 통해 산출된 예측 데이터도 제공되고 있는데, 우선은 과거 데이터부터 찾아 보았다. 과거 데이터들도 다양한 소스와 방식으로 제공되는데 이 중에서 가장 대표적인 것은 ERA5라는 데이터이다.
ERA5의 5는 5세대를 말하는 것 같고, RA는 재분석(Reanalysis)되었다는 것에서 따온 것 같다. 수일 전 데이터도 있지만 재분석을 통해 보다 정밀하게 보정된 데이터를 제공하는 것 같다.
일부는 유럽 지역에 국한된 데이터도 있지만, 기본적인 데이터들은 전 지구를 대상으로 하고 있고 데이터에 따라 조금씩 다르기는 하지만 1940년 또는 1950년 이후의 기상, 해양, 토양, 강수, 기온, 바람 등과 관련한 다양한 데이터를 제공하고 있다. 👍
데이터를 다운로드하기 위해서는 회원가입 및 로그인이 필요한데, 과도한 정보를 요구하지는 않는다.
이 중에서, “ERA5-Land monthly averaged data from 1950 to present”를 한 번 살펴보면, 1950년부터 월별로 기온 등의 다양한 정보를 제공하고 있다.
“Download data” 탭을 클릭해서, 필요로 하는 측정값 항목, 연도, 월, 지역 범위(좌표) 등을 선택하고, “Show API request” 버튼을 클릭하면, Python에서 해당 조건의 데이터를 내려받을 수 있는 코드를 제공한다.
그러나, 파이썬에서 코드에 제시된 cdsapi 패키지를 설치하고, 위의 코드를 실행하면 에러가 나는데, API Key 정보가 없기 때문이다. 위에서 언급한 코페르니쿠스 웹사이트 우측 상단의 자신의 이름 부분을 클릭하면 나오는 사용자 정보 화면에서 API key에 대한 정보를 확인할 수 있다.
메모장 등의 텍스트 편집기에서 아래와 같은 구성으로 작성하고 C:\Users\사용자계정 폴더에 “.cdsapirc”이라는 파일명으로 저장한 다음 코드를 실행하면 정상적으로 다운로드가 실행된다. 파일은 grib 포맷이다.
url: https://cds.climate.copernicus.eu/api/v2
key: UID(6자리 숫자):API Key (예시 : 4XXXX:4XXXfXXX-XXXf-4XXX-9XXX-7XXXebXXfdXX)API key 환경설정이 어려운 분들은 아래 페이지를 참조하시기 바랍니다.
ERA5-Land monthly averaged data from 1950 to present 데이터 시각화
이제 본격적인 코딩이다. (몇 줄 안된다)
# Download data 탭에서 조건 설정 후 Show API request 버튼 누르면 제공되는 코드
import cdsapi
c = cdsapi.Client()
c.retrieve(
'reanalysis-era5-land-monthly-means',
{
'product_type': 'monthly_averaged_reanalysis',
'variable': [
'10m_u_component_of_wind', '10m_v_component_of_wind', '2m_temperature',
'total_precipitation',
],
'year': [
'1950', '1951', '1952',
'1953', '1954', '1955',
'1956', '1957', '1958',
'1959', '1960', '1961',
'1962', '1963', '1964',
'1965', '1966', '1967',
'1968', '1969', '1970',
'1971', '1972', '1973',
'1974', '1975', '1976',
'1977', '1978', '1979',
'1980', '1981', '1982',
'1983', '1984', '1985',
'1986', '1987', '1988',
'1989', '1990', '1991',
'1992', '1993', '1994',
'1995', '1996', '1997',
'1998', '1999', '2000',
'2001', '2002', '2003',
'2004', '2005', '2006',
'2007', '2008', '2009',
'2010', '2011', '2012',
'2013', '2014', '2015',
'2016', '2017', '2018',
'2019', '2020', '2021',
'2022',
],
'month': '07',
'time': '00:00',
'format': 'grib',
},
'data/era5-land-monthly-means.grib')연월별 전세계 데이터이므로 시간이 좀 걸린다.
# xarray 패키지 임포트
# 아래를 보면 의존성이 있는 아래의 패키지를 먼저 설치해야 한다.
# !pip install numpy packaging pandas
# !pip install xarray
# https://docs.xarray.dev/en/stable/getting-started-guide/installing.html
import xarray as xr 그런데, 내 환경에서는 netCDF와 달리 grib는 바로 열리지 않고 에러가 발생했다. 계속된 구글링과 시행착오 끝에 해결한 방법은 grib 파일을 처리해주는 패키지를 추가로 설치해준 것이다.
# grib 처리를 위한 패키지 설치
#!pip install cfgrib ecmwflibs ecCodes이제 파일을 열어보자. 약 1.65GB 파일을 여는데 약 9.3초가 걸렸다.
ds = xr.open_dataset('data/era5-land-monthly-means.grib', engine='cfgrib', filter_by_keys={'stepType': 'avgid'})
# 데이터를 조회해보면
ds
전체 속성은 50개 정도인데 다운로드 받을 때 7월/0시/ 3개 변수만 선택했기 때문 에1950년 7월1일 0시~2022년 7월 1일 0시의 전 지구 범위의 u10, v10, t2m 변수에 대한 값들을 가지고 있다. 전체 월/시/변수를 모두 선택했다면 파일 용량이 어마어마할 것이다. 반대로, 다운로드 받을 때 입력한 좌표 범위를 만족하는 데이터만 받을 수도 있다.
지구 전체의 수십 년간의 월별 데이터여서 약 1.65GB이기 때문에 우선은 조건 선택(필터링) 기능을 이용해서 한반도 지역 좌표 범위에 해당하는 데이터만 추출했다.
# 한반도 지역 데이터만 추출
ds_kor = ds.where((ds.longitude > 124) & (ds.longitude < 131) & (ds.latitude < 39) & (ds.latitude > 32), drop=True) 나와 같이 pandas가 익숙한 분들을 위해 xarray에서 pandas dataframe으로 변환해주는 기능을 제공하고 있다. 이를 이용해서 데이터 내용을 확인해보자.
df = ds_kor.to_dataframe()
df.head(15)
데이터프레임으로 변환한 상태에서는 쉽게 작업할 수 있겠고, csv 등으로 저장해서 쓸 수도 있다. 우선은 xarray 상태에서 더 진행해보자. 아래는 연도별로 데이터를 분리하는 코드이다.
import pandas as pd
for i in range(1950, 2023):
i2 = pd.to_datetime(str(i) + '-07-01T00:00:00.000000000')
nm = '_kor' + str(i) #연도별 데이터명
# 연도별 데이터명 객체에 연도별 데이터만 추출해서 적용
globals()['ds{}'.format(str(nm))] = ds_kor.where(ds_kor.time == i2, drop=True) 연도별 데이터를 시각화해보자. 먼저 시각화에 필요한 패키지 및 관련한 환경을 준비한다.
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
mpl.rc('font', family='Malgun Gothic') #한글 폰트 적용시
plt.rcParams["figure.figsize"] = (6,6) #차트 사이즈아래는 maptplotlib에서 배경지도를 표시할 수 있는 패키지이다.
# pip install basemap
from mpl_toolkits.basemap import Basemap, cm베이스맵 패키지를 이용해서 베이스맵을 선언한다. 좌표범위는 아까 한반도 추출시의 좌표 범위와 같다.
m= Basemap(llcrnrlon=min(ds_kor.longitude),llcrnrlat=min(ds_kor.latitude),urcrnrlon=max(ds_kor.longitude),urcrnrlat=max(ds_kor.latitude), resolution='i',projection='cyl',lon_0=(max(ds_kor.longitude)-min(ds_kor.longitude)) / 2,lat_0=(max(ds_kor.latitude)-min(ds_kor.latitude)) / 2)1950년 7월 1일의 t2m(지상 2m의 온도)를 시각화해보자. m의 drawcoastlines를 이용해서 해안선도 그려준다.
ds_kor1950.t2m.plot()
m.drawcoastlines()
plt.text(124, 39.5,"1950년 07월 평균 기온", fontdict={'weight': 'bold', 'size': 14})
plt.savefig('image/ds_kor_1950.png', bbox_inches='tight', pad_inches=0.1 )
그런데, t2m 항목의 온도가 절대온도 기준이다. pandas처럼 xarray에서 쉽게 새로운 항목을 만들면서 기존 항목의 값을 연산해서 넣어줄 수 있을까?
ds_kor['t2m2'] = ds_kor['t2m'] - 273.15ds_kor을 조회해보니 t2m2가 생겼고, 원하는 대로 값이 들어갔다!!!
연도별로 분할하는 코드부터 다시 실행하고 나서 이제는 반복문을 이용해서 한 번에 연도별 이미지를 만들자.
l = [ds_kor1950, ds_kor1951, ds_kor1952, ds_kor1953, ds_kor1954, ds_kor1955, ds_kor1956, ds_kor1957, ds_kor1958, ds_kor1959, ds_kor1960, ds_kor1961, ds_kor1962, ds_kor1963, ds_kor1964, ds_kor1965, ds_kor1966, ds_kor1967, ds_kor1968, ds_kor1969, ds_kor1970, ds_kor1971, ds_kor1972, ds_kor1973, ds_kor1974, ds_kor1975, ds_kor1976, ds_kor1977, ds_kor1978, ds_kor1979, ds_kor1980, ds_kor1981, ds_kor1982, ds_kor1983, ds_kor1984, ds_kor1985, ds_kor1986, ds_kor1987, ds_kor1988, ds_kor1989, ds_kor1990, ds_kor1991, ds_kor1992, ds_kor1993, ds_kor1994, ds_kor1995, ds_kor1996, ds_kor1997, ds_kor1998, ds_kor1999, ds_kor2000, ds_kor2001, ds_kor2002, ds_kor2003, ds_kor2004, ds_kor2005, ds_kor2006, ds_kor2007, ds_kor2008, ds_kor2009, ds_kor2010, ds_kor2011, ds_kor2012, ds_kor2013, ds_kor2014, ds_kor2015, ds_kor2016, ds_kor2017, ds_kor2018, ds_kor2019, ds_kor2020, ds_kor2021, ds_kor2022]
k = 1950
for j in l:
j.t2m.plot()
m.drawcoastlines()
plt.text(124, 39.5, str(k) + "년 07월 평균 기온", fontdict={'weight': 'bold', 'size': 14})
plt.savefig('image/' + str(k) + '.png', bbox_inches='tight', pad_inches=0.1) #tight
plt.clf() # 기존 차트에 누적해서 생기는 걸 막아줌
k += 1
이제 이 이미지들을 animated GIF로 묶어보자. 먼저, 필요한 패키지들을 부르고, image 폴더안에 생성된 연도별 지도 이미지들을 리스트에 추가한다.
from pathlib import Path
import imageio #!pip install imageio
image_path = Path('image') #특정 폴더/경로 지정
images = list(image_path.glob('*.png')) #해당 경로의 특정 포맷에 대한 (이미지) 파일명 리스트 생성
image_list = []
for file_name in images: #Animated GIF 생성용 리스트에 추가
image_list.append(imageio.imread(file_name))이제 이미지 리스트를 gif로 만들면 된다.
# Animated GIF 생성
#fps가 1이면 1초에 한 프레임씩. 2면 1초에 2프레임씩
imageio.mimwrite('image/kor_tempmap.gif', image_list, fps=1, loop=False ) 최근에는 2018년의 기온이 높았던 것으로 보인다. 물론, 시도 행정구역에 따라 다를 수는 있겠다.
grib 데이터를 netCDF 등 다른 포맷으로 저장할 수도 있다.
# netCDF 등의 다른 포맷으로 저장 가능
ds_kor.to_netcdf("data/ds_kor.nc", engine="netcdf4")내용은 별 거 없지만 글이 생각보다 길어졌다. 앞으로도 해볼 수 있는 게 많을 것 같다.
- 벡터로 마스킹해서 특정 행정구역 데이터만 추출해서 평균값을 내볼 수 있을까?
- ERA5 등의 과거 데이터말고, 예보 데이터도 다룰 수 있을까? 어떤 차이가 있을까?
- netCDF 포맷의 데이터도 다뤄보자
- 미국이나 우리나라 데이터도 살펴보자.
우선 오늘은 여기까지….
