[파이썬] 해양수치모델 검증 - 1. nc 파일읽기, 특정 지점(위치) 찾기, 최근접 격자 찾기

 오늘은 수치예측모델의 정확도를 분석하기 위해 관측정보와 비교하는 과정에 대해 소개해보도록 하겠다. 이를 위해서는 모델자료를 읽고, 관측점에 해당하는 데이터를 추출해야한다. 또한, txt 또는 csv형식의 관측정보를 읽어야 하고, 모델과 관측의 데이터를 동일한 시간간격으로 맞춰줘야 한다.

 

 이번 글에서는 모델자료를 읽는 법과 특정 위치의 데이터만 추출하는 법에 대해 다뤄보도록 하겠다.

 

 

---

| 자료 개요

• 기간 : 2016. 10. 02 ~ 10. 08 (태풍 차바)
• 모델 : HYCOM(해양), RDAPS(기상), RWW3(파랑)
• 관측 : 기상청 서귀포 해양기상부이, 국립해양조사원 서귀포 조위관측소
• 분석항목 : 해양 - 조위, 기상 - 풍속, 파랑 - 파고

 

 분석 샘플로 활용할 자료는 2010년 이후 가장 큰 피해를 입힌 태풍 중 하나인 2016년 10월의 차바 기간의 자료를 활용하였다. 수치모델자료는 3가지 정보를 활용하였고, 해양은 hycom, 기상은 기상청 rdaps, 파랑은 기상청 rww3를 활용했다. 사전에 다운로드 받은 각 자료를 동일한 해상도로 regrid 한뒤, 병합하여 한 파일로 생성하였고, 이 자료를 분석에 활용하였다. 

 

 각 형식 변환 및 시간변환은 이전글에 설명한적이 있으니(해양, 기상 정보 및 분석툴 카테고리), 참고하면 되고, regrid는 esmf의 ncl코드를 활용하였다. regrid는 추후 자세히 소개해보도록 하겠다.

 

 관측정보는 기상청 해양기상부이와 국립해양조사원 조위관측소 자료를 활용하였으며, 태풍 경로에 위치한 서귀포 정점의 정보를 활용하였다.

 

 

---

| 수치모델자료 읽기 netcdf read

 수치모델 자료는 위에서 언급했듯이 hycom, rdaps, rww3를 동일 시간, 동일 격자로 regrid 한 뒤 병합한 netcdf 파일을 사용하였다. netcdf 자료에 대한 기본적인 내용 및 읽는 방법에 대해서는 이전에 작성한 적이 있으니 아래 글을 참고하면 될 것 같다.

 

[파이썬] 해양수치모델 NetCDF 자료처리. 파이썬 Nc 파일 읽기 ncread

 

 

---

1) 라이브러리(NetCDF4)  불러오기

 분석을 위해서 NetCDF를 다룰 수 있는 netcdf4, 그 외 기본적인 계산 및 연산에 필요한 numpy, pandas, datetime, matplotlib, glob, os 등도 불러와야한다. 이 라이브러리는 어떤 분석을 해도 기본적으로 필요한 라이브러리들이니, 항상 코드 상단에 박아놓는 편이다.

import numpy as np
import pandas as pd
import datetime
from matplotlib import dates
import matplotlib.ticker as ticker
import matplotlib.dates as mdates
import re
import matplotlib.pyplot as plt
import os
import glob
import netCDF4 as nc

 

 

 

---

2) 모델자료 및 각 변수 읽기

ncdata = nc.Dataset('파일명') // ncvars = ncdata['변수명']

 위 링크의 글에서 설명한대로 netcdf4 모듈의 Dataset 함수를 활용한다. 모델데이터 변수명을 mdata로 가정하면, mdata = nc.Dataset('파일명') 과 같이 입력한다. 모델 파일의 data를 mdata로 불러들이게 되고 mdata 안의 각 변수를 읽으면 된다. mdata는 print(mdata)로 아래 그림과 같이 정보를 확인할 수 있다.

파이썬 nc파일 변수 확인

 

 여기서 분석할 변수인 조위, 파고, 풍속 변수는 zos, VHM0, UGRD & VGRD(U, V)이다. 각 변수와 시간(time), 위도(lat), 경도(lon)도 함께 읽는다.

path = '/home/opr/jspark/00_PROJECT/2022y/02_1MW/02_product/01_typhoon/20161005/'
ncname = 'typhoon_20161005_v2.nc'
mdata = nc.Dataset(path+ncname)


mtime = mdata['time'][:];
mlat = mdata['lat'][:];
mlon = mdata['lon'][:];
windu = mdata['UGRD'][:];
windv = mdata['VGRD'][:];
hm = mdata['VHM0'][:];
tide = mdata['zos'][:];

 

 

 mdata['변수명'] 뒤에 [:]를 붙이는 건 해당 변수의 데이터만 읽어온다는 의미이다. 위 코드처럼 변수를 읽어오면, 각 변수 모든 차원의 데이터를 읽어오고, 각 변수별 차원은 아래와 같다.

• mlat(1208) : 위도 (y, 1208개 격자)
• mlon(1805) : 경도 (x, 1805개 격자)
• mtime(169) : 시간 (t, 169개 time step)
• hm(169, 1205, 1805) : 파고 (t, y, x), 1805x1205의 cell이 169개 time step으로 구성
• windu(169, 1205, 1805) : 풍속u (t, y, x), 1805x1205의 cell이 169개 time step으로 구성
• windv(169, 1205, 1805) : 풍속u (t, y, x), 1805x1205의 cell이 169개 time step으로 구성
• time(169, 1205, 1805) : 조위 (t, y, x), 1805x1205의 cell이 169개 time step으로 구성

파이썬 nc데이터 변수 및 차원

 

 위 사진에서처럼 수치모델 각 변수의 차원은 t(시간), y(위도), x(경도)의 3차원 데이터로 구성되어 있다. 보통 우리가 익숙한 x, y, t 차원이다. 1205x1805격자의 데이터가 총 169타임스텝만큼 있다는 의미이다.

 

 관측 데이터와 비교/분석을 하려면 1205x1805개의 데이터 중에서 관측점 위치에 해당하는 1개 격자의 데이터만 있으면 된다. 아래에서 해당 정점의 위치만 찾는 방법에 대해 소개해보도록 하겠다.

 

 

---

 3) 특정 위치의 데이터 추출

- 가장 가까운 격자 찾기 : 관측과 모델의 위경도 최소값

 비교하려는 관측점을 덕적도 해양기상부이로 가정해보도록 하겠다. 덕적도 부이의 위치는 위도 37.2361, 경도 126.0188이다. 따라서, 1205(위도) x 1805(경도)의 격자 중 덕적도 부이와 가장 가까운 위치의 격자를 찾으면 된다.

 

 가장 가까운 격자는 모델 위도값(mlat, 1205개) 중 관측 위도인 37.2361과의 차이가 가장 작은 격자, 모델 경도값(mlon, 1805개) 중 관측 경도인 126.0188과의 차이가 가장 작은 격자 번호를 찾으면 된다. 

 

 코드 수식으로 표현하면, mlat에서 37.2361을 뺀 값 중 최소값을 구하면 되고, mlon에서 126.0188을 뺀 값 중 최소값을 구하면 된다. 뺄셈을 할 때, 해당 값은 거리이므로 절대값을 씌워야하고, 각 위도와 경도의 최소값이 위치한 격자번호를 찾으면 된다. 코드로 구현하면 아래와 같다.

## 덕적도 부이 위치 => ploc
ploc = [37.23611, 126.01875]#duck juck

## obs lat, lon 정의
olat = ploc[0]
olon = ploc[1]

## 위도/경도 차이 최소값 구하기, 최소값이 위치한 격자번호 찾기
bbb = np.array(abs(mlat - olat))
ccc = np.array(abs(mlon - olon))
bbb2 = np.where(bbb==min(bbb))
ccc2 = np.where(ccc==min(ccc))

 

 

 위도/경도 차이의 최소값이 위치한 격자번호는 위도 737, 경도 963이다. 따라서, 덕적도 부이와 가장 가까운 격자는 1205x1805 데이터 중 737, 963의 데이터가 된다. 아래 오른쪽 사진의 mlat, mlon을 보면 737, 963의 데이터는 37.2425, 126.014로 덕적도 부이와 가장 가까운 위치라는 것을 확인할 수 있다.

파이썬 수치모델 가장 가까운 위치 찾기

 

 

 

---

 - 특정 위치의 데이터만 추출 // ncvars = ncdata['변수명'][t, y, x]

 관측점과 가장 가까운 격자는 y(위도) 737번, x(경도) 963번 격자이다. 따라서, [169 x 1205 x 1805]의 3차원 데이터 중 [:, 737, 963] 값만 추출하면 된다. 따라서, 모델 변수를 읽을 때, 아래와 같이 읽을 변수의 위치를 지정해주면 된다.

파이썬 특정 위치의 데이터 추출

 

 

 

---

***** 위경도 값이 행렬 구조인 경우 - 최근접점, 최근접 격자 찾기 *****

• 직선거리 :    ((x2 - x1)^2 + (y2 - y1)^2)^(1/2)

 수치모델 중 위도와 경도가 1차원이 아닌 2차원 행렬 매트릭스 구조인 경우가 있다. 주로 기상모델에서 이런 경우가 많은데, 모델 도메인(영역)을 직교좌표격자가 아닌 구면격자로 나눈 데이터들이 그러하다. 아래 그림처럼 각 격자의 행 또는 열이 같더라도, 해당 격자의 위도 또는 경도의 값이 상이한 경우이다.

수치모델 격자 형식이 상이한 위도와 경도

수치모델 격자 형식이 상이한 위도와 경도 데이터

 

 

 이러한 경우에는 모델과 관측점 위/경도의 최소값이 아닌 각 격자의 위치와 관측점의 위치와의 직선거리를 구하고, 직선거리의 최소값을 구해서 최근접점을 찾아야 한다.

 

 다들 알다시피, 직선거리는 (x2 - x1)^2 + (y2 - y1)^2 의 제곱근으로 구한다. x1과 x2에는 관측, 모델의 경도, y1와 y2에는 관측, 모델의 위도를 넣어주면 된다. 따라서, (mlon - 37.2361)^2 + (mlat - 126.0188)^2의 제곱근을 구하면 된다. 

 

 

• 행렬 차원 및 배열(구조) 변환 : numpy.reshape(data,(차원))
• 3x5 행렬로 변환시, numpy.reshape(data,(3,5))
• 다차원 행렬을 1차원으로 변환시, numpy.reshape(data, -1)

 코드로 구현하기 위해서는 mlon을 먼저, 1차원 행렬로 변환해야 한다. 위 데이터에서 예로 든 기상모델의 위경도는 419 x 491의 격자로 구성되어 있는데, 이를 1차원 행렬인 205729 x 1의 행렬로 바꿔줘야한다. 이는 numpy의 reshape를 활용한다. 그럼 아래와 같이, 1차원 행렬로 변형할 수 있고, 경도도 동일하게 변형해준다. numpy의 reshape함수는 변환할 데이터와 차원을 입력하면 된다. 차원은 (n, n)과 같이 입력하면 되지만, 다차원 행렬을 1차원 행렬로 변형할 때는 -1을 입력하면 된다.

파이썬 위도와 경도 행렬변환 - reshape

 

 

1차원으로 바꾼 위도와 경도를 활용해서 관측 위치와의 직선거리를 구하면 된다. 모델 위도를 mlat3, 경도를 mlon3, 관측 위도를 olat, 관측 경도를 olon이라고 하면, 수식은 아래와 같다.

line_distance = ((mlat3-olat)**2 + (mlon3-olon)**2)**(1/2)

 

 

 이제, 직선거리를 구한 line_distance의 최소값을 구하면 된다. 아래와 같이 최소값과 최소값이 위치한 index를 구하면 된다.

ind2 = np.where(line_distance==min(line_distance))

파이썬 위도와 경도, 최근접점 찾기

 

 그럼 관측위치와 가장 가까운 격자는 205729개 중 99427번째 격자라는 것을 알 수 있다. 하지만 이 99427번이라는 격자는 419 x 491 행렬을 1차원인 205729x1로 바꾼 행렬의 index이다. 따라서, 이 99427번이 다시 2차원 행렬로 바꿨을 때의 index를 알아야 한다.

 

 방법은 간단하다. np.reshape를 활용하거나 단순 계산을 해도 된다. 아래 행렬 구조를 바꾸는 reshape의 개념을 나타내보았다. 205729x1 행렬을 원래대로(419 x 491) reshape 하게되면 첫 491개는 1행, 492~982까지는 2행.... 과 같이 재배치하게 된다. 따라서, 99427번째 데이터를 491로 나눈 몫이 행번호, 나머지가 열번호가 된다. 따라서, 최근접점의 인덱스는 (202, 245)가 된다.

파이썬 위도와 경도, 행렬 변환

 

 

 만약, 계산이 틀렸을 수도 있으니, 데이터를 확인해보면 된다. 아래 이미지 모델 위도와 경도의 데이터다. 해당 격의 위도는 37.2115, 경도는 125.98이다. 덕적도 부이 위치인 위도 37.2361, 경도 126.0188과 가장 가까운 격자라는 것을 확인할 수 있다.

파이썬 위도와 경도, 최근접점 찾기

 

---

 여기까지 관측과 모델값을 비교하기 위해서 모델데이터의 특정 위치 정보를 추출하는 방법에 대해 알아보았다. 다음에는 실제 비교/분석을 위해서 time step/scale을 맞추고, 비교해보는 과정에 대해 소개해보도록 하겠다.