Mayavi 라이브러리를 이용해서 전구 기압 분포 3D 시각화

오늘은 전구 기압 분포를 3차원 구에 시각화해보겠습니다

---

이번 포스팅을 하게된 계기는 Spherical harmonics를 공부하던 중 전구 기압 분포를 3차원 구에 시각화한 뒤

 

애니메이션으로 회전시키면서 보면 좋을 것 같아서 진행해봤습니다

 

먼저 구를 만들기 위해서는 구면 좌표계를 사용해야합니다

 

이는 plt를 활용하여 그릴때도 똑같이 적용됩니다

 

구면 좌표계는 위키 백과의 정의로

 

"3차원 공간 상의 점들을 나타내는 좌표계의 하나로, 보통 (r, Φ, Θ) 로 나타낸다" 라고 설명되어 있습니다

 

각 성분들은 r은 원점으로부터의 거리, Φ 는 방위각, Θ는 고도각으로 3차원 공간인 구에서 위치를 표현할 수 있습니다

 

구면좌표계는 경우에 따라 좌표값이 한 점을 여러 좌표가 가리키는 경우가 있기 때문에 변수들의 범위는 아래처럼 제한됩니다

 

r ≥ 0 

0 ≤ Θ ≤ π

0 ≤ Φ ≤ 2π

 

여기서 π는 원주율이 아닌 각도의 라디안 값을 나타내며 180˚를 의미합니다

 

일반적으로 저희가 보는 3차원 x, y, z 좌표계는 구면 좌표계로 변환시 아래와 같은 과정을 거칩니다

 

x = r sinΘ cosΦ

y = r sinΘ sinΦ

z = r cosΘ

 

이는 아래 그림을 보면 이유를 아실 수 있습니다

출처 : https://rcasio80.blogspot.com/2016/10/3d.html

---

코드입니다

#-----------------------------------------------------------------------
#Library Import
#
from mayavi import mlab
import numpy as np
from scipy.special import sph_harm
import pygrib as pg
from scipy.ndimage import gaussian_filter
from mayavi.tools.animator import animate

#
#-----------------------------------------------------------------------
#file preprocessing
grbs = pg.open('/Users/lsh/coding/python/fnl_20230717_12_00.grib2')

grb = grbs.select(name = 'Pressure reduced to MSL')[0]

pres = grb.values/100
#
#-----------------------------------------------------------------------
# Create a sphere
r = 1.5
pi = np.pi
cos = np.cos
sin = np.sin

phi, theta = np.mgrid[0:pi:181j, 0:2*pi:360j]
#
#-----------------------------------------------------------------------
# 카테시안 좌표계에서 구면 좌표계로 변환
x = r * sin(phi) * cos(theta)
y = r * sin(phi) * sin(theta)
z = r * cos(phi)
#
#-----------------------------------------------------------------------
#
mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(900, 600))
mlab.clf()
m = mlab.mesh(x, y, z, scalars=pres, colormap = 'summer')
mlab.colorbar(orientation = 'vertical', nb_labels = 10, label_fmt = '%.1f')

@mlab.animate
def anim():
  f = mlab.gcf()
  while 1:
    f.scene.camera.azimuth(10)
    f.scene.render()
    yield
    
a = anim()

mlab.show()
#-----------------------------------------------------------------------

 

코드는 간단합니다

 

여기서 위에 설명과 연결지어서 봐야합니다 

 

헷갈리실 수도 있는 부분이 좌표계를 변경하는 부분입니다

 

위에 설명에는

 

x = r sinΘ cosΦ

y = r sinΘ sinΦ

이렇게 나와 있지만 코드에는 반대로 되어있습니다

 

이는 자료를 기록하고 만드는데 규칙과 관계가 있습니다

 

보통 기압, 기온과 같은 스칼라 형태의 기상자료는 (time, lat, lon) 이렇게 3가지 차원을 가집니다

 

제가 가지고 온 grib2 파일도 같은 방식이고 netCDF 역시 마찬가지입니다

 

하지만 우리가 변환할 좌표계는 (x, y)순서지만 자료가 기록된 파일은 (위도, 경도) 순으로 작성이 되어 있습니다

 

그렇기 때문에 np.mgird를 사용해서 그리드 격자를 생성할 때도 360j, 181j가 아닌 181j, 360j인 것이죠

 

직접 차원을 NCL로 확인해보겠습니다

 

제가 가져온 Pressure reduced to MSL의 변수 정보입니다

 

고도의 영향을 제거하는 감소된 기압으로 관측소에 실제로 존재하는 온도 프로파일을 사용하여 관측소 바로 아래

 

지점에 해수면이 존재할 수 있는지 추정된 압력값입니다

 

변수명 뒤에 (lat_0, lon_0)이렇게 쓰여져있는 모습입니다

 

 

실제로 봐도 181칸과 360칸의 그리드로 구성된 모습이죠

 

그렇기 때문에 좌표계 변환시 방위각과 고도각을 반대로 써줘야 제대로된 구가 만들어집니다

 

원래대로 쓰시면 어떤 결과가 나오는지는 직접 해보시면 좋을 것 같습니다(비틀어진 종이 모양이 나옵니다!)

---

 

여기서 중요한 점은 mayavi 설치 시 pyqt5를 반드시 설치해야하는데 이는 코랩 환경에서는 실행이 안되니

 

반드시 개인 작업 환경에서 하시고 맥북이 있으신 분들은 그대로 작업을 진행하셔도 됩니다!

 

윈도우 상에서는 혹시나 wsl 환경을 쓰시는 분들은 반드시 xming을 설치하셔서

 

x-window가 켜질 수 있도록 해주셔야합니다 (이것은 저도 아직 해결을 못했어요ㅠㅠ 해결하는대로 포스팅해보겠습니다!)

 

저렇게 코드를 실행하면 

이렇게 작은 창이 뜨게되고 Delay 값을 작게 해줄수록 지구가 빠르게 돌아갑니다

 

f.scene.camera.azimuth(10) 이 부분에서 저는 초당 회전하는 방위각을 10도로 해두었지만 

 

이것을 줄이면 더 천천히 돌아가게 됩니다

 

저는 Delay를 1로 해주었구여 돌아가는 건 아래 영상으로 확인해보시면 되겠습니다

 

 

전구 기압 분포는 정확히 958.3723hPa부터 남극 부분에 1093.42hPa 까지 존재하고 있었습니다

 

남극 부근에는 아주 강한 고기압이 자리를 잡고 있었습니다

 

남극에 이렇게 강한 고기압이 나타난 이유를 추정해보면

 

해당 자료는 6월 15일 자료로 북반구의 여름이 오고 있으면서 남반구는 점점 겨울이 오는 시기입니다

 

이 시기에 남극의 온도는 더욱 떨어지고 온도가 낮아질수록 공기괴는 수축하여 공기가 누르는 무게가 

 

더욱 커지기 때문에 지표에서 이렇게 강한 고기압이 나타났다고 생각합니다

---

 

이렇게 3d 시각화도 도전해보았습니다

 

다음 해볼것은 같은 mayavi를 이용해서 3차원 공간의 바람 흐름을 그려보는 것입니다

 

셤 기간이라 빠르게 하지는 못할 것 같고 틈틈이 시간 내서 조금씩 해보겠습니다!

 

긴 글 읽어주셔서 감사합니다!