[python] 북극 진동 (Arctic Oscillation) 패턴 계산 & 그리기 2: 자료 수집과 처리

- 개요

이번 포스트에서는 저번에 정한 자료 수집, 자료 처리 목표를 달성합니다.

목표는 자료 수집에서 1개, 자료 처리에서 3개로 총 4가지입니다.

 

- 자료 수집

목표 1. NCEP2 재분석 자료 수집: NOAA 사이트에서 월평균 1000 hPa 지위고도 다운로드 (2, 3)

먼저 아래 링크로 들어갑니다.

https://psl.noaa.gov/data/gridded/data.ncep.reanalysis2.html

다음 그림의 설명대로 하시면 NCEP2 재분석 자료인 hgt.mon.mean.nc파일을 다운로드 받을 수 있습니다.

자료수집

 

 

- 자료 처리

목표 2. 자료 추출: 북위 20도 이상, 1979~2000년 (4, 8)

nc파일을 읽기 위해 xarray 라이브러리를 씁니다.

isel과 sel 그리고 slice를 이용하여 원하는 시간과 영역의 지위고도장을 추출합니다.

import xarray as xr

data = xr.open_dataset("hgt.mon.mean.nc")

syr = 1979
fyr = 2000
nyr = fyr - syr + 1
hgtRaw = data['hgt'].isel(level=0).sel(lat=slice(None, 20), time=slice(f"{syr}-01", f"{fyr}-12")).values

 

isel을 쓰면 특정 coorindate에서 원하는 인덱스의 자료를 추출할 수 있고,

sel과 slice를 혼합하여 이용하면 원하는 coordinate 범위를 추출할 수 있습니다.

 

'hgt'의 coordinate는 [time, level, lat, lon]입니다.

hgtRaw 행의 isel(level=0)은 level의 0번째값이 1000 hPa로 1000 hPa을 추출합니다.

sel에서는 slice 방법을 쓰는데요.

NCEP2 자료는 lat이 90부터 -90까지 내림차순입니다.

lat=slice(None, 20)으로 설정하면 처음부터 위도 20도까지 추출하겠다는 뜻입니다. 이러면 북위 90도부터 20도까지 추출할 수 있습니다.

다음으로 time=slice() 부분은 syr=1979, fyr=2000이므로 1979-01(1979년 1월)부터 2000-12(2000년 12월)까지 자료를 추출합니다.

 

목표 3. 계절 변동을 제거하여 아노말리 구하기 (5, 9)

계절 변동을 제거한다는 것은 모든 기간 각 월에 대해 평균값(기후값)을 구하고 각각의 월평균값에 이 평균값을 뺀다는 것입니다. 보통 실제값(raw value)과 평균값의 차이를 아노말리라고 부릅니다.

예를 들어 1979년 1월의 아노말리는 1979년 1월의 실제값에서 1979~2000년 1월 평균값을 빼서 구합니다.

import numpy as np
"""
평균값과 아노말리값을 저장할 배열의 차원 설정을 위해 lat, lon 배열의 크기 선언
"""
lat = data['lat'].sel(lat=slice(None, 20)).values
lon = data['lon'].values
nlat, nlon = len(lat), len(lon)

"""
hgtClim 변수에는 모든 기간에 대해 1, 2, 3 ... 10, 11, 12월의 평균값을 저장
hgtAno는 실제값(hgtRaw)과 평균값(hgtClim)의 차이
1년에 무조건 12달이 있으므로 i%12를 활용
"""
hgtClim = np.zeros((12, nlat, nlon))
for i in range(nyr*12):
    hgtClim[i%12, :, :] += hgtRaw[i, :, :]
hgtClim /= nyr

hgtAno = np.zeros(hgtRaw.shape)
for i in range(nyr*12):
    hgtAno[i, :, :] = hgtRaw[i,:,:] - hgtClim[i%12, :,:]

 

 

목표 4. 영역별 가중치 적용 (11)

영역별 가중치를 적용하는 이유는 위도에 따라 격자의 간격이 다르기 때문입니다.

적도에서는 위도 1도가 약 110 km지만 북위 37도 즈음되는 한국에서는 위도 1도가 88 km죠.

이렇게 위도가 올라갈수록 위도의 실제 간격은 좁아지므로 고위도에서는 같은 영역에 더 많은 격자점이 들어갑니다.

만약 영역별 가중치를 적용하지 않으면 고위도인 북극 지역은 좁은 지역에 많은 점이 있으므로 EOF분석 원리상 북극의 변동성이 더 크게 잡힙니다.

이를 보정하기 위해 NOAA에서는 sqrt(cos(위도))를 곱해줍니다.

cos 함수는 0도에서 1, 90도에서 0값을 반환하므로 이 값을 곱하면 위도가 높아질수록 같은 영역에 점이 많아지는 효과를 상쇄합니다.

lat_rad = np.radians(lat)  # 위도를 라디안으로 변환
area_weighting = np.sqrt(np.cos(lat_rad))

"""
hgtAno의 차원은 [시간, 위도, 경도]
area_weighting의 차원은 [위도]

area_weighting[:, np.newaxis]로 hgtAno에 곱하면
area_weighting의 차원은 [위도, 경도]가 되면서 hgtAno의 각 시간마다 곱해짐
(실제 프로그래밍 상에서 변수가 이렇게 변하는지는 확인 안해봤으나 설명은 맞음)
"""
hgtAnoWgt = hgtAno * area_weighting[:, np.newaxis]

 

주의사항

여기서 제 코드에서는 아래의 오류가 발생합니다.

RuntimeWarning: invalid value encountered in sqrt

area_weighting = np.sqrt(np.cos(lat_rad))

cos(90도)는 원래 0입니다만

np.sqrt(np.cos(lat_rad)) <<< 이 값을 확인하면 위도가 90도인 첫 번째 부분의 값은 -4.3*10^-8이입니다.

10^-8이면 사실상 0에 가까운 값이긴한데 음수이기 때문에 sqrt를 취할 수 없습니다. (복소수까지 고려한 변수 타입이면 가능할듯)

그래서 위도 90도인 위치에서는 값이 nan이 나옵니다.

아무튼 위도 90도에 위치한 값들은 0이 곱해져서 모든 시간에서 값이 0이 되어야하죠.

모든 시간에 대해 0만 있는 자료는 EOF 분석에 영향을 줄 수 없으므로 nan인 상태로 계산해도 무방합니다.

이게 좀 그렇다면 hgtAno 변수에서 위도의 첫 번째 위치에 0을 넣어주시면 됩니다.

 

 

여기까지 자료 수집과 처리 부분의 4가지 목표를 달성했습니다.

다음 포스트에서는 EOF 분석을 해보겠습니다.