오늘은 기상청에서 제공하는 수치모델자료를 시각화 해보고자 한다. 활용하는 자료는 기상자료개방포털에서 제공하는 경량화 수치모델 자료를 활용한다.

자료수집

자료 시각화

활용 패키지

활용패키지는 목록은 아래와 같다.

library(ncdf4)#nc file 로드
library(dplyr)# 핸들링 
library(tidyverse)# 핸들링 
library(jsonlite)# leaflet 객체로 활용
library(raster)
library(leaflet)
library(leaflet.extras2)

자료 추출

NC파일을 불러오고 데이터를 먼저 추출해보자.

r=ncdf4::nc_open('text/g120_v070_erea_unis_20kind.2022031400.nc')
names(r$var)# 변수 목록 확인

lat=ncvar_get(r,'latitude')
lon=ncvar_get(r,'longitude')
U=ncvar_get(r,names(r$var)[6])# 10m Wind U Component
V=ncvar_get(r,names(r$var)[7])# 10m Wind V Component
temp=ncvar_get(r,names(r$var)[8])# 5m Temperature
ncdf4::nc_close(r)

raw_df=data.frame(lat=c(lat),lon=c(lon),
        U=c(U[,,1]),V=c(V[,,1]),temp=c(temp[,,1]))

위치정보를 포함한 격자로 변환

경량화 수치모델 파일은 행491 열 419로 격자형태의 자료이나, 위치정보는 구면좌표계의 자료이므로 바로 레스터화할 수 없다. 레스터화 하기 위해서는

from_crs="+proj=longlat +datum=WGS84 +no_defs"
r <- raster(nrow=491*1,ncol=419*1,crs=from_crs,
            xmn=min(lon),xmx=max(lon),ymn=min(lat),ymx=max(lat))

coordinates(raw_df)=~lon+lat
crs(raw_df)=from_crs
ru <- rasterize(raw_df, r, 'U', fun=mean)
rv <- rasterize(raw_df, r, 'V', fun=mean)
rtemp <- rasterize(raw_df, r, 'temp', fun=mean)

시각화

leaflet() %>%
  addProviderTiles('Esri.WorldImagery', group = 'Esri.WorldImagery')%>%
  setView(lng = 127.7319, lat = 36.096, zoom = 6)%>%
  addRasterImage(rtemp,  opacity = 0.3)#opacity : 투명도 조절

풍향풍속도 시각화

풍향풍속도는 아래 두 링크를 참조하였다. 먼저 아래와 같이 header와 data를 json파일로 만들어야 한다.

header에서 parameterCategory는 2 고정값이며, parameterNumber는 U는 2 V는 3 la1, la2는 각각 위도(latitude)의 최소,최대값, lo1, lo2는 경도(Longitude)의 최소, 최대값, nx는 가로 격자 수, ny는 세로 격자 수, dx는 가로 격자 간격, dy는 세로 격자 간격, refTime은 관측시간을 입력하면 된다.

[Solved] leaflet.extras2 Velocity not producing anything

jzadra Asks: leaflet.extras2 Velocity not producing anything I am attempting to produce a leaflet velocity layer (see GitHub - onaci/leaflet-velocity: Visualise velocity data on a leaflet layer) using R (a port to R exists as part of leaflet.extras2 package: GitHub -...

SolveForum

Visualizing wind using Leaflet

Description of collection, formation and visualization of wind data using Leaflet plugin leaflet-velocity

WolfblogMats Melander

header 생성

자료추출에서 만든 raw_df와 rasterize한 U와 V를 활용해보자.

lo1=extent(raw_df)[1]
lo2=extent(raw_df)[2]
la1=extent(raw_df)[3]
la2=extent(raw_df)[4]

nx=r@nrows
ny=r@ncols
dx = (lo2 - lo1) / nx
dy = (la1 - la2 ) / ny
parameterCategory=2
parameterNumber=c(2,3)
refTime="2009-01-01 14:00:00"
header=data.frame(la1,la2,lo1,lo2,nx,ny,dx,dy,parameterCategory,
                  parameterNumber,refTime)

data 생성

raw_df=data.frame(coordinates(r),
        U=ru@data@values,V=rv@data@values)
colnames(raw_df)=c('lon','lat',"U","V")

#데이터프레임 격자화
raw_df2=raw_df%>%complete(
  lat, lon, fill = list(U = 0,V=0))

# 결측된 격자 채우기
raw_df2=raw_df2%>% fill(everything(), .direction = "down")

#자료 정렬
raw_df2=raw_df2%>% arrange(desc(lon), desc(lat))

data <- vector(mode = "list", length = 2)
data[[1]] <- raw_df2$U
data[[2]] <- raw_df2$V

wave_content  <- data.frame(I(data))
wave_content$header=header
wavejson=toJSON(wave_content)

시각화

풍향/풍속 시각화는 아래와 같이  실행하면 된다.

leaflet() %>%
  addProviderTiles('Esri.WorldImagery', group = 'Esri.WorldImagery')%>%
  setView(lng = 127.7319, lat = 36.096, zoom = 3)%>%
  addVelocity(content = wavejson, group = "velo", layerId = "veloid") %>%
  addLayersControl(baseGroups = "base", overlayGroups = "velo")

조금 더 해보기

raster패키지의 getData를 통해 고도나 지적도를 받을 수 있다. 지적도 자료는 공공데이터포털에서 받는게 조금 더 정밀한 것 같지만 편의상 해당 자료를 활용하겠다.

library(raster)
grid <- merge(raster::getData('SRTM', lon=126,lat=37),
              raster::getData('SRTM', lon=126,lat=35))
grid3=aggregate(grid, fact=50)
kor=raster::getData('GADM',country='KOR',level=0)

위에서 만든 기온 레스터를 우리나라의 크기에 맞게 잘라서 시각화해보자.

rtemp_crop=mask(crop(rtemp,kor),kor)
values(rtemp_crop)=values(rtemp_crop)-273.15
grid.pal=colorNumeric(palette = "RdYlGn",values(rtemp_crop),
                      na.color = 'transparent', reverse = TRUE)

leaflet() %>%
  addProviderTiles('Esri.WorldImagery', group = 'Esri.WorldImagery')%>%
  setView(lng = 127.7319, lat = 36.096, zoom = 6)%>%
  addRasterImage(rtemp_crop, colors = grid.pal, opacity = 0.8)%>%
addLegend(pal = grid.pal, values = values(rtemp_crop),title = "기온(도)")%>%
  addVelocity(content = wavejson, group = "velo", layerId = "veloid")