[R] ECMWF Interim 수치 예측 모델로부터 지오포텐셜 고도 (Geopotential Height)를 이용한 주성분 분석 및 가시화

 정보

• 업무명     : ECMWF Interim 수치 예측 모델로부터 지오포텐셜 고도 (Geopotential Height)를 이용한 주성분 분석 및 가시화

• 작성자     : 이상호

• 작성일     : 2020-03-15

• 설   명      :

• 수정이력 :

 

 내용

[개요]

• 안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.

• 현대 기상 예측에서 수치예보는 40 %로서 큰 비중을 차지합니다. 기상청에 따르면 "수치 예보는 물리학의 방정식에 의해 바람과 기온 등의 기온 변화를 컴퓨터로 계산하여 미래의 대기 상태를 예측하는 방법"이라고 합니다.
  

• 즉 컴퓨터를 통해 지구의 상태를 시뮬레이션하여 기상 예측하는 기술입니다. 이러한 시뮬레이션 방법은 여러 종류가 있으나 결과는 "GPV (Grid Point Value)" 형식으로 제공하는 것이 대부분입니다.
  

• 기상 예측은 시간과 격자에 대한 여러 기상 정보를 내포하고 "NetCDF" 형식으로 배포하고 있습니다.

• 이러한 NetCDF 형식은 대기과학에서 자주 사용되기 때문에 프로그래밍 언어 (R, Python, GrADS, NCL)에 따른 자료 처리 및 가시화가 필연입니다. 이와 관련한 정보는 링크를 참조해주시기 바랍니다.

[R] "ncdf4" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (2)

 

[R] "RNetCDF" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (3)

 

[Python] 파이썬 NetCDF 형식인 Aqua/CERES 기상위성 자료를 이용한 가시화

 

[GrADS, ShellScript] 그라즈 NetCDF 형식인 NOAA 월평균 해수면 온도 자료를 이용한 가시화

 

[NCL] NetCDF 형식 및 NCL 변수 구조

 

• 그리고 주성분 분석은 경험적 직교함수 (Empirical Orthogonal Function, EOF)로서 기상 자료에 자유도를 최소화하면서 원래의 자료 지닌 현상을 간단하게 묘사할 수 있습니다. 이와 관련한 정보는 링크를 참조하시기 바랍니다.

[강릉원주대 대기환경과학과] 2015년 2학기 전필 고급기상통계학 소개 및 과제물

 

• 따라서 ECMWF Interim 수치 예측 모델 자료로부터 지오포텐셜 고도를 이용한 주성분 분석 및 가시화를 소개해 드리고자 합니다. 
  

 

[특징]

• 주성분 분석을 이해하기 위해서 자료 처리 및 가시화가 요구되며 이 프로그램은 이러한 목적을 달성하기 위한 소프트웨어

 

[기능]

• NetCDF 파일 읽기

• 다중 자료에 대한 데이터 프레임 처리

• 고유근에 따른 설명력 전처리 및 가시화

• 1-5번째 주성분에 대한 고유벡터 전처리 및 가시화 

• 1-5번째 주성분에 대한 지오포텐셜 기압 전처리 및 가시화

• 1-3번째 주성분 점수 시계열 
  

 

[활용 자료]

• 자료명 : interim_daily_06uvz_200901_15.nc

• 자료 종류 : 

• 확장자 : NetCDF

• 영역 : 전지구

• 기간 : 2009년 01월 15일 - 2018년 12월 15일 (매월 15일)

• 시간 해상도 : 일 1개

• 제공처 : ECMWF

 

 

[자료 처리 방안 및 활용 분석 기법]

• 없음

 

[사용법]

• 소스 코드 참조

 

[사용 OS]

• Windows 10

 

[사용 언어]

• R v3.6.3

• R Studio v1.2.5033

 

 소스 코드

[명세]

• 전역 설정

  • 최대 10 자리 설정

  • 메모리 해제

  • 영어 인코딩 설정

  • 폰트 설정

# Set Option
memory.limit(size = 9999999999999)
options(digits = 10)
Sys.setlocale("LC_TIME", "english")
font = "Palatino Linotype"
cbMatlab =  colorRamps::matlab.like(11)

 

• 라이브러리 읽기

# Library Load
library(RNetCDF)
library(tidyverse)
library(extrafont)
library(lubridate)
library(data.table)
library(noncompliance)
library(factoextra)
library(metR)
library(ncdf4)
library(netf4.helpers)
library(colorRamps)
library(RColorBrewer)
library(synoptReg)

 

• NetCDF 파일 설정 및 열기

  • RNetCDF::open.nc를 통해 파일 열기

sFileDirName = Sys.glob("INPUT/interim_daily_06uvz_2009_2018_1.5/*.nc")

iCount = 1

oOpenNc = RNetCDF::open.nc(sFileDirName[iCount])

 

 

• 메타 데이터 이해

  • RNetCDF 패키지에서 print.nc 함수를 사용하여 파일의 메타 데이터를 확인

  • 이러한 메타 데이터에서 dimensions은 위도, 경도 및 시간과 같은 벡터 형식 변수가 포함

  • 추가로 메타 데이터에서 variables마다 세부 정보  포함

# NetCDF Meta Information
RNetCDF::print.nc(oOpenNc)

 

> RNetCDF::print.nc(oOpenNc)

netcdf offset64 {
dimensions:
	longitude = 54 ;
	latitude = 34 ;
	time = UNLIMITED ; // (31 currently)
variables:
	NC_FLOAT longitude(longitude) ;
		NC_CHAR longitude:units = "degrees_east" ;
		NC_CHAR longitude:long_name = "longitude" ;
	NC_FLOAT latitude(latitude) ;
		NC_CHAR latitude:units = "degrees_north" ;
		NC_CHAR latitude:long_name = "latitude" ;
	NC_INT time(time) ;
		NC_CHAR time:units = "hours since 1900-01-01 00:00:00.0" ;
		NC_CHAR time:long_name = "time" ;
		NC_CHAR time:calendar = "gregorian" ;
	NC_SHORT z(longitude, latitude, time) ;
		NC_DOUBLE z:scale_factor = 0.092301838863626 ;
		NC_DOUBLE z:add_offset = 13405.5261147056 ;
		NC_SHORT z:_FillValue = -32767 ;
		NC_SHORT z:missing_value = -32767 ;
		NC_CHAR z:units = "m**2 s**-2" ;
		NC_CHAR z:long_name = "Geopotential" ;
		NC_CHAR z:standard_name = "geopotential" ;
	NC_SHORT u(longitude, latitude, time) ;
		NC_DOUBLE u:scale_factor = 0.000992487215187029 ;
		NC_DOUBLE u:add_offset = 1.37864651762288 ;
		NC_SHORT u:_FillValue = -32767 ;
		NC_SHORT u:missing_value = -32767 ;
		NC_CHAR u:units = "m s**-1" ;
		NC_CHAR u:long_name = "U component of wind" ;
		NC_CHAR u:standard_name = "eastward_wind" ;
	NC_SHORT v(longitude, latitude, time) ;
		NC_DOUBLE v:scale_factor = 0.000970069021761937 ;
		NC_DOUBLE v:add_offset = 5.54898395962974 ;
		NC_SHORT v:_FillValue = -32767 ;
		NC_SHORT v:missing_value = -32767 ;
		NC_CHAR v:units = "m s**-1" ;
		NC_CHAR v:long_name = "V component of wind" ;
		NC_CHAR v:standard_name = "northward_wind" ;

// global attributes:
		NC_CHAR :Conventions = "CF-1.6" ;
		NC_CHAR :history = "2019-10-22 07:39:12 GMT by grib_to_netcdf-2.13.0: grib_to_netcdf /data/data03/scratch/34/80/_mars-atls02-a562cefde8a29a7288fa0b8b7f9413f7-5cmpQI.grib -o /data/data02/scratch/55/a1/_grib2netcdf-atls18-a82bacafb5c306db76464bc7e824bb75-5d9XvN.nc -utime" ;
}

 

• NetCDF 파일 읽기

  • RNetCDF::read.nc를 통해 파일 읽기

# NetCDF File Read 
liReadNc = RNetCDF::read.nc(oOpenNc)

 

• 속성 정보 가져오기

  • RNetCDF::dim.inq.nc 및 RNetCDF::att.get.nc를 통해 시간 속성 정보 가져오기
  • RNetCDF::att.get.nc를 통해 해수면 온도에 대한 scale_factor, add_offset, _FillValue 가져오기

# Get Time Unit 
nFileNumber = RNetCDF::dim.inq.nc(oOpenNc, "time")$length
sTimeUnits = RNetCDF::att.get.nc(oOpenNc, "time", "units")

# Get Attribute 
fSCaleFactor = RNetCDF::att.get.nc(oOpenNc, "z", "scale_factor") 
fAddOffset = RNetCDF::att.get.nc(oOpenNc, "z", "add_offset")
fFillValue = RNetCDF::att.get.nc(oOpenNc, "z", "_FillValue")

 

 

• 시간 변수 가져오기

  • lubridate::make_datetime를 통해 숫자형 (year, month, day, hour, minute, second)을 날짜형 (dtDateTime)으로 변환

  • format를 통해 날짜형 (dtDateTime)을 문자형 (sColName)으로 변환

# Get Time Value
nTime = RNetCDF::utcal.nc(sTimeUnits, liReadNc$time)

dfDateTime = data.frame(nTime) %>%
   dplyr::mutate(
      dtDateTime = lubridate::make_datetime(year, month, day, hour, minute, second)
      , sColName = format(dtDateTime, "%Y-%m-%d")
   )

dplyr::tbl_df(dfDateTime)

 

 

• 변수 가져오기

  • 경도 (nLon), 위도 (nLat), 지오포텐셜 고도 (nVal)

  • 지오포텐셜 고도의 경우 Fill Value의 여부에 따라 Scale Factor 및 Add Off Set 그리고 NA로 설정

# Get Value
nLon = liReadNc[["longitude"]]
nLat = liReadNc[["latitude"]]
nTmpVal = liReadNc[["z"]]
nVal = ifelse(nTmpVal != fFillValue, (nTmpVal * fScaleFactor) + fAddOffset, NA)
   
dim(nTime)
dim(nLon)
dim(nLat)
dim(nVal)

 

 

• 데이터 프레임 (Data Frame) 변환

  • 주성분 분석을 위해 행 (위/경도에 따른 지오포텐셜 고도), 열 (시간)에 대해서 변환

 

# Convert 2D Matrix to 1D Matrix
maData = matrix(nVal, nrow = dim(nLon) * dim(nLat), ncol = nFileNumber)

dfData = data.frame(maData)
colnames(dfData) = c(dfDateTime$sColName)

dplyr::tbl_df(dfData)

 

 

• 다중 자료를 위한 데이터 프레임 (Data Frame) 열 추가 및 L1 처리

  • 단일 자료 (iCount == 1)의 경우 noncompliance::expand.grid.DT를 통해 위도, 경도에 따른 데이터 프레임을 설정하고 dplyr::bind_cols를 통해 열 추가 (위도, 경도, 시간에 따른 지오포텐셜 고도)

  • 다중 자료에서는 dplyr::bind_cols를통해 시간에 따른 지오포텐셜 고도를 추가 

  • 이는 벡터화를 사용하기 때문에 성능 저하 문제 해결

# Add Coloum
if (iCount == 1) {
   dfGeoData = noncompliance::expand.grid.DT(
      nLat, nLon
      , col.names = c("nLat", "nLon")
   )
   
   dfDataL1 = dplyr::bind_cols(dfGeoData, dfData)
} else {
   dfDataL1 = dplyr::bind_cols(dfDataL1, dfData)
}

dplyr::tbl_df(dfDataL1)

 

 

• Data Frame L1을 이용한 L2 전처리 

  • 전체 기간 (2009년 1월 15일 - 2018년 12월 15일)에서 필요한 사례에 대해 선택하는 과정

  • 즉 열 이름 (sColName)으로 필터하기 위해 필요한 사례 데이터셋 (dfFilterData)를 이용하여 날짜형 (dtDateTime) 및 문자형 (sColName)으로 변환

  • dplyr::select를 통해 두 데이터셋 (dtDataL1, dfFilterDataL1)에서 위도, 경도, 특정 사례에 대한 열을 선택 

# L2 Processing Using Data Frame L1 : Filter Exception Case
dfFilterData = data.table::fread("INPUT/total_sb_date_step6_gn.csv", header = FALSE, sep = ",")

dfFilterDataL1 = dfFilterData %>%
   dplyr::mutate(
      dtDateTime = lubridate::make_datetime(V1, V2, V3)
      , sColName = format(dtDateTime, "%Y-%m-%d")
   )

dplyr::tbl_df(dfFilterDataL1)

dfDataL2 = dfDataL1 %>%
   dplyr::select(one_of(c("nLon", "nLat", dfFilterDataL1$sColName)))

dplyr::tbl_df(dfDataL2)

 

 

 

• Data Frame L2을 이용한 L3 전처리

  • dplyr::filter를 통해 동아시아에 대한 위도 (20-50), 경도 (100-140)를 사용 

# L3 Processing Using Data Frame L2
dfDataL3 = dfDataL2 %>%
   dplyr::filter(
      dplyr::between(nLon, 100, 140)
      , dplyr::between(nLat, 20, 50)
      )

dplyr::tbl_df(dfDataL3)

 

 

• Data Frame L3을 이용한 L4 전처리 

  • 주성분 분석을  위해 위도, 경도를 제외한 데이터 셋 구성

# L4 Processing Using Data Frame L3
dfDataL4 = dfDataL3[ , -c(1:2)]

dplyr::tbl_df(dfDataL4)

 

 

• Data Frame L4을 이용한 L5 전처리 

  • 주성분 분석을  위해 전치 행렬로 변환

 

# L5 Processing Using Data Frame L4
dfDataL5 = t(dfDataL4) %>%
   as.data.frame()

dplyr::tbl_df(dfDataL5)

 

 

• 주성분 분석

  • 수행 결과 총 7종 (sdev, loadings, center, scale, n.obs, scores, call)

liPcaDec = synoptReg::pca_decision(dfDataL5)

summary(liPcaDec)

 

 

• 고유근에 따른 설명력

  • 수행 결과 총 7종 (sdev, loadings, center, scale, n.obs, scores, call)

liPcaDec = synoptReg::pca_decision(dfDataL5)

summary(liPcaDec)

 

 

 

• 가시화

  • 데이터 세트에서 정보는 여전히 숨겨져 있기 때문에 가시화 필요

  • 가시화는 일반적인 정적 탐색 데이터 분석에서 웹 브라우저의 동적 대화식 데이터 시각화에 이르기까지 다양함

  • 특히 R의 기본 plot으로 여러 미학적 측면을 제어 할 수 있으나 Hadley Wickham (2016)이 개발한 ggplot2는 새로운 방법으로 시각화하기 때문에 이를 사용 

 

• ggplot2를 이용한 고유근에 따른 설명력

# Eigen Value
nEigVal  = liPcaDec$PCA$sdev^2  

# Visualization Using ggplot2
liEof = princomp(dfDataL5, cor = TRUE)
nEigValPer = (nEigVal / sum(nEigVal, na.rm = TRUE)) * 100.0
nEigValCumTop10 = cumsum(nEigValPer)[1:10]

factoextra::fviz_screeplot(
   liEof
   , addlabels = TRUE
   , x = "Principal  Components"
   , y = "Percentage of Variance [%]"
   , ncp = 10
   ) + 
   geom_point(aes(y = nEigValCumTop10), col = "red") + 
   geom_line(aes(y = nEigValCumTop10), col = "red") +
   theme_bw() +
   theme(
      plot.title = element_text(face = "bold", size = 18, color = "black")
      , axis.title.x = element_text(face = "bold", size = 18, colour = "black")
      , axis.title.y = element_text(face = "bold", size =18, colour = "black", angle=90)
      , axis.text.x  = element_text(face = "bold", size = 18, colour = "black")
      , axis.text.y  = element_text(face = "bold", size = 18, colour = "black")
      , legend.title = element_text(face = "bold", size = 14, colour = "white")
      , legend.position = c(0, 1)
      , legend.justification = c(0, 0.96)
      , legend.key = element_blank()
      , legend.text = element_text(size = 14, face = "bold", colour = "white")
      , legend.background = element_blank()
      , text = element_text(family = font)
      , plot.margin = unit(c(0, 8, 0, 0), "mm")
   ) +
   ggsave(filename = paste0("FIG/PrinComp1.png"), width = 12, height = 8, dpi = 600)

 

 

• ggplot2를 이용한 주성분에 따른 고유벡터 가시화

  • 주성분 분석 결과 (liPcaDec)에서 loadings를 이용한 고유벡터 데이터 셋 (nEigVec) 구성

  • 즉 위도, 경도, 1-540번째 주성분에 따른 고유벡터로 구성

# Eigen Vector
nEigVec = liPcaDec$PCA$loadings %>%
   as.data.frame.matrix()

dfDataL6 = data.frame(
   nLon = dfDataL3$nLon
   , nLat = dfDataL3$nLat
   , nEigVec
   )

dplyr::tbl_df(dfDataL6)

 

 

• ggplot2를 이용한 1번째 주성분에 대한 고유벡터

# Visualization Using ggplot2
ggplot() +
   theme_bw() +
   geom_raster(data = dfDataL6, aes(x = nLon, y = nLat, fill =  Comp.1), interpolate = TRUE) +
   metR::geom_contour2(data = dfDataL6, aes(x = nLon, y = nLat, z = Comp.1), color = "black", alpha = 0.3) +
   metR::geom_text_contour(data = dfDataL6, aes(x = nLon, y = nLat, z = Comp.1), stroke = 0.2, check_overlap = TRUE, rotate = TRUE, na.rm = TRUE) +
   scale_fill_gradient2(limits=c(-0.1, 0.1), breaks = seq(-0.1, 0.1, 0.05)) +
   geom_path(data = ggplot2::map_data("world2"), aes(long, lat, group = group), color = "black") +
   metR::scale_x_longitude(expand = c(0, 0), breaks = seq(100, 140, 10), limits = c(100, 140)) +
   metR::scale_y_latitude(expand = c(0, 0), breaks = seq(20, 50, 10), limits = c(20, 50)) +
   labs(
      x = ""
      , y = ""
      , fill = "Eigen Vector"
      , colour = ""
      , title  = "ERA Interim Geopotential Height 150 hPa"
      , subtitle = "Period : January 15, 2009 - December 15, 2018"
      , caption = "Source : ECMWF"
   ) +
   theme(
      plot.title = element_text(face = "bold", size = 18, color = "black")
      , axis.title.x = element_text(face = "bold", size = 18, colour = "black")
      , axis.title.y = element_text(face = "bold", size =18, colour = "black", angle=90)
      , axis.text.x  = element_text(face = "bold", size = 18, colour = "black")
      , axis.text.y  = element_text(face = "bold", size = 18, colour = "black")
      , legend.title = element_text(face = "bold", size = 14, colour = "black")
      , legend.position = c(0.01, 0.99)
      , legend.justification = c(0, 0.96)
      , legend.key = element_blank()
      , legend.text = element_text(size = 14, face = "bold", colour = "black")
      , legend.background = element_blank()
      , text = element_text(family = font)
      , plot.margin = unit(c(0, 8, 0, 0), "mm")
   ) +
   ggsave(filename = paste0("FIG/PrinComp2_Comp.1.png"), width = 12, height = 8, dpi = 600)

 

 

• ggplot2 및 반복문을 이용한 1-5번째 주성분에 대한 고유벡터

# Visualization Using For Loop and ggplot2
for (iCount in 1:5) {
 
   sVal = paste0("Comp.", iCount)
   sFigDirName = paste0("FIG/PrinComp2_Comp.", iCount, ".png")
   
   ggplot() +
      theme_bw() +
      geom_raster(data = dfDataL6, aes(x = nLon, y = nLat, fill = !! as.name(sVal)), interpolate = TRUE) +
      metR::geom_contour2(data = dfDataL6, aes(x = nLon, y = nLat, z = !! as.name(sVal)), color = "black", alpha = 0.3) +
      metR::geom_text_contour(data = dfDataL6, aes(x = nLon, y = nLat, z = !! as.name(sVal)), stroke = 0.2, check_overlap = TRUE, rotate = TRUE, na.rm = TRUE) +
      scale_fill_gradient2(limits=c(-0.1, 0.1), breaks = seq(-0.1, 0.1, 0.05)) +
      geom_path(data = ggplot2::map_data("world2"), aes(long, lat, group = group), color = "black") +
      metR::scale_x_longitude(expand = c(0, 0), breaks = seq(100, 140, 10), limits = c(100, 140)) +
      metR::scale_y_latitude(expand = c(0, 0), breaks = seq(20, 50, 10), limits = c(20, 50)) +
      labs(
         x = ""
         , y = ""
         , fill = "Eigen Vector"
         , colour = ""
         , title  = "ERA Interim Geopotential Height 150 hPa"
         , subtitle = "Period : January 15, 2009 - December 15, 2018"
         , caption = "Source : ECMWF"
      ) +
      theme(
         plot.title = element_text(face = "bold", size = 18, color = "black")
         , axis.title.x = element_text(face = "bold", size = 18, colour = "black")
         , axis.title.y = element_text(face = "bold", size =18, colour = "black", angle=90)
         , axis.text.x  = element_text(face = "bold", size = 18, colour = "black")
         , axis.text.y  = element_text(face = "bold", size = 18, colour = "black")
         , legend.title = element_text(face = "bold", size = 14, colour = "black")
         , legend.position = c(0.01, 0.99)
         , legend.justification = c(0, 0.96)
         , legend.key = element_blank()
         , legend.text = element_text(size = 14, face = "bold", colour = "black")
         , legend.background = element_blank()
         , text = element_text(family = font)
         , plot.margin = unit(c(0, 8, 0, 0), "mm")
      ) +
      ggsave(filename = sFigDirName, width = 12, height = 8, dpi = 600)
}

 

 

• Data Frame L5를 이용한 주성분에 따른 지오포텐셜 기압 전처리 (1)

  • synoptReg::synoptclas를 통해 종관 분류를 계산

    • 이 과정에서 분류를 위한 임계값 (ncomp)을 -5 < ncomp < 5 설정

  • synoptReg::raster_clas를 통해 raster 형식으로 변환되고 이를 데이터 프레임으로 재 변환

  • 즉 위도, 경도, 1-10번째 주성분에 대한 지오포텐셜 고도

liPcaCla = synoptReg::synoptclas(dfDataL5, ncomp = 5)

dfTmpDataL7 = synoptReg::raster_clas(
   longitude = unique(dfDataL3$nLon)
   , latitude = unique(dfDataL3$nLat)
   , grouped_data = liPcaCla$grouped_data
   ) %>%
   raster::as.data.frame()

dplyr::tbl_df(dfTmpDataL7)

dfDataL7 = data.frame(
   nLon = dfDataL3$nLon
   , nLat = dfDataL3$nLat
   , dfTmpDataL7
   ) 

dplyr::tbl_df(dfDataL7)

 

 

 

• Data Frame L7를 이용한 주성분에 따른 지오포텐셜 기압 전처리 (2)

  • dplyr::select를 통해 위도, 경도, 1-5번째 주성분에 대한 지오포텐셜 고도 선택

  • 위도, 경도, 1-5번째 주성분에 대한 키 (sKey) 및 값 (nVal)으로 정제하기 위해 dplyr::gather 이용

  • dplyr::mutate를 통해 지오포텐셜 고도를 기압으로 변환 

dfDataL8 = dfDataL7 %>%
   dplyr::select(nLon, nLat, CT1:CT5) %>%
   tidyr::gather(key = "sKey", value = "nVal", -nLon, -nLat) %>% 
   dplyr::mutate(nVal = nVal / 100) 

dplyr::tbl_df(dfDataL8)

 

 

• ggplot2를 이용한 1번째 주성분에 대한 지오포텐셜 기압 가시화

dfDataL9 = dfDataL8 %>%
   dplyr::filter(sKey == "CT1")

dplyr::tbl_df(dfDataL9)

# Visualization Using ggplot2
ggplot() +
   theme_bw() +
   geom_raster(data = dfDataL9, aes(x = nLon, y = nLat, fill =  nVal), interpolate = TRUE) +
   metR::geom_contour2(data = dfDataL9, aes(x = nLon, y = nLat, z = nVal), color = "black", alpha = 0.3) +
   metR::geom_text_contour(data = dfDataL9, aes(x = nLon, y = nLat, z = nVal), stroke = 0.2, check_overlap = TRUE, rotate = TRUE, na.rm = TRUE) +
   scale_fill_gradientn(colours = cbMatlab, limits=c(130, 150), breaks = seq(130, 150, 10), na.value = cbMatlab[length(cbMatlab)]) +
   geom_path(data = ggplot2::map_data("world2"), aes(long, lat, group = group), color = "black") +
   metR::scale_x_longitude(expand = c(0, 0), breaks = seq(100, 140, 10), limits = c(100, 140)) +
   metR::scale_y_latitude(expand = c(0, 0), breaks = seq(20, 50, 10), limits = c(20, 50)) +
   labs(
      x = ""
      , y = ""
      , fill = "Perssure [hPa]"
      , colour = ""
      , title  = "ERA Interim Geopotential Height 150 hPa"
      , subtitle = "Period : January 15, 2009 - December 15, 2018"
      , caption = "Source : ECMWF"
   ) +
   theme(
      plot.title = element_text(face = "bold", size = 18, color = "black")
      , axis.title.x = element_text(face = "bold", size = 18, colour = "black")
      , axis.title.y = element_text(face = "bold", size =18, colour = "black", angle=90)
      , axis.text.x  = element_text(face = "bold", size = 18, colour = "black")
      , axis.text.y  = element_text(face = "bold", size = 18, colour = "black")
      , legend.title = element_text(face = "bold", size = 14, colour = "black")
      , legend.position = c(0.01, 0.99)
      , legend.justification = c(0, 0.96)
      , legend.key = element_blank()
      , legend.text = element_text(size = 14, face = "bold", colour = "black")
      , legend.background = element_blank()
      , text = element_text(family = font)
      , plot.margin = unit(c(0, 8, 0, 0), "mm")
   ) +
   ggsave(filename = paste0("FIG/PrinComp3_Comp.1.png"), width = 12, height = 8, dpi = 600)

 

 

 

• ggplot2 및 facet_wrap을 이용한 1-5번째 주성분에 대한 지오포텐셜 기압 가시화

# Visualization Using ggplot2
ggplot() +
   theme_bw() +
   geom_raster(data = dfDataL8, aes(x = nLon, y = nLat, fill =  nVal), interpolate = TRUE) +
   metR::geom_contour2(data = dfDataL8, aes(x = nLon, y = nLat, z = nVal), color = "black", alpha = 0.3) +
   metR::geom_text_contour(data = dfDataL8, aes(x = nLon, y = nLat, z = nVal), stroke = 0.2, check_overlap = TRUE, rotate = TRUE, na.rm = TRUE) +
   scale_fill_gradientn(colours = cbMatlab, limits=c(130, 150), breaks = seq(130, 150, 10), na.value = cbMatlab[length(cbMatlab)]) +
   geom_path(data = ggplot2::map_data("world2"), aes(long, lat, group = group), color = "black") +
   metR::scale_x_longitude(expand = c(0, 0), breaks = seq(100, 140, 10), limits = c(100, 140)) +
   metR::scale_y_latitude(expand = c(0, 0), breaks = seq(20, 50, 10), limits = c(20, 50)) +
   labs(
      x = ""
      , y = ""
      , fill = "Pressure [hPa]"
      , colour = ""
      , title  = "ERA Interim Geopotential Height 150 hPa"
      , subtitle = "Period : January 15, 2009 - December 15, 2018"
      , caption = "Source : ECMWF"
   ) +
   theme(
      plot.title = element_text(face = "bold", size = 18, color = "black")
      , axis.title.x = element_text(face = "bold", size = 18, colour = "black")
      , axis.title.y = element_text(face = "bold", size =18, colour = "black", angle=90)
      , axis.text.x  = element_text(face = "bold", size = 18, colour = "black")
      , axis.text.y  = element_text(face = "bold", size = 18, colour = "black")
      , legend.title = element_text(face = "bold", size = 18, colour = "black")
      , legend.position = c(0.6, 0.15)
      , legend.key = element_blank()
      , legend.text = element_text(size = 18, face = "bold", colour = "black")
      , legend.background = element_blank()
      , text = element_text(family = font)
      , plot.margin = unit(c(0, 8, 0, 0), "mm")
   ) +
   facet_wrap(. ~ sKey, ncol = 2) +
   ggsave(filename = paste0("FIG/PrinComp3_Comp.3.png"), width = 12, height = 12, dpi = 600)

 

 

• 주성분 결과 (liPcaDec)를 이용한 주성분 점수 시계열

  • tibble::rownames_to_column를 통해 날짜 행을 열로 추가

  • 위도, 경도, 1-3번째 주성분 점수에 대한 키 (sKey) 및 값 (nVal)으로 정제하기 위해 dplyr::select 및 dplyr::gather 이용

  • dplyr::filter를 통해 2018년 01월 01일 이후 사례 선정 

dfTmpDataL10 = liPcaDec$PCA$scores %>%
   as.data.frame()

dfDataL10 = dfTmpDataL10 %>%
   tibble::rownames_to_column() %>%
   dplyr::select(Comp.1:Comp.3, rowname) %>%
   tidyr::gather(key = "sKey", value = "nVal", -rowname) %>%
   dplyr::mutate(
      dtDate = readr::parse_date(rowname)
      , nYear = lubridate::year(dtDate)
   ) %>%
   dplyr::filter(dtDate >= lubridate::as_date("2018-01-01"))

dplyr::tbl_df(dfDataL10)

 

 

• ggplot2를 이용한 주성분 점수 시계열

# Visualization Using ggplot2
sMinDate = lubridate::as_date("2018-01-01")
sMaxDate = lubridate::as_date("2019-01-01")

ggplot(data = dfDataL10, aes(x = dtDate, y = nVal, color = sKey)) +
   theme_bw() +
   geom_point() +
   geom_line() +
   scale_x_date(expand = c(0, 0), date_breaks = "month", date_labels = "%b-%d\n%Y", limits = as.Date(c(sMinDate, sMaxDate))) +
   scale_y_continuous(expand = c(0, 0), breaks = seq(-50, 50, 20), limits=c(-50, 50)) + 
   scale_color_discrete( 
      breaks = c("Comp.1", "Comp.2", "Comp.3") 
      , labels = c("PC1 (40.4 %)", "PC2 (18.3 %)", "PC3 (14.2 %)") 
      ) + 
   labs( 
      x = "Date [Month-Day Year]" 
      , y = "Principal Component Analysis Scores"
      , subtitle = "Period : February 19, 2018 - December 31, 2018"
      , caption = "Source : ECMWF"
      , fill = "" 
      , colour = "" 
   
      ) + 
   theme( 
      plot.title = element_text(face = "bold", size = 18, color = "black") 
      , axis.title.x = element_text(face = "bold", size = 18, colour = "black") 
      , axis.title.y = element_text(face = "bold", size=18, colour = "black", angle = 90) 
      , axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size=18, colour = "black") 
      , axis.text.y = element_text(face = "bold", size=18, colour = "black") 
      , legend.position = c(0.08, 0.95)
      # , legend.justification = c(0, 0.96)
      , legend.key = element_blank() 
      , legend.text = element_text(size = 14, face = "bold") 
      , legend.title = element_text(face = "bold", size = 14, colour = "black") 
      , legend.background=element_blank() 
      , text=element_text(family = font) 
      , plot.margin = unit(c(0, 8, 0, 0), "mm") 
      ) + 
   ggsave(filename = paste0("FIG/PrinComp4.png"), width = 12, height = 8, dpi = 600)

 

 

[전체]

 

 참고 문헌

[논문]

• 없음

[보고서]

• 없음

[URL]

• 없음

 

 문의사항

[기상학/프로그래밍 언어]

• sangho.lee.1990@gmail.com

[해양학/천문학/빅데이터]

• saimang0804@gmail.com