정보
-
업무명 : 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)
정보 업무명 : "ncdf4" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (2) 작성자 : 이상호 작성일 : 2020-02-13 설 명 : 수정이력 : 내용 [개요]..
shlee1990.tistory.com
[R] "RNetCDF" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (3)
정보 업무명 : "RNetCDF" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (3) 작성자 : 이상호 작성일 : 2020-03-15 설 명 : 수정이력 : 내용 [개..
shlee1990.tistory.com
[Python] 파이썬 NetCDF 형식인 Aqua/CERES 기상위성 자료를 이용한 가시화
정보 업무명 : NetCDF 형식인 Aqua/CERES 기상위성 자료를 이용한 가시화 작성자 : 이상호 작성일 : 2019-10-22 설 명 : 수정이력 : 2020-02-05 : 소스 코드 명세 추가 내용 [개요] 안녕하세요? 웹 개발 및 연구..
shlee1990.tistory.com
[GrADS, ShellScript] 그라즈 NetCDF 형식인 NOAA 월평균 해수면 온도 자료를 이용한 가시화
정보 업무명 : 그라즈 NetCDF 형식인 NOAA 월평균 해수면 온도 자료를 이용한 가시화 작성자 : 이상호 작성일 : 2019-09-07 설 명 : 수정이력 : 내용 [특징] NetCDF 결과를 이해하기 위해 가시화 도구가 필요하며..
shlee1990.tistory.com
[NCL] NetCDF 형식 및 NCL 변수 구조
정보 업무명 : NetCDF 형식 및 NCL 구조 작성자 : 이상호 작성일 : 2020-02-27 설 명 : 수정이력 : 내용 [개요] 안녕하세요? 웹 개발 및 연구 개발을 담당하고 있는 해솔입니다. NCL (NCAR Command Language)은 미..
shlee1990.tistory.com
-
그리고 주성분 분석은 경험적 직교함수 (Empirical Orthogonal Function, EOF)로서 기상 자료에 자유도를 최소화하면서 원래의 자료 지닌 현상을 간단하게 묘사할 수 있습니다. 이와 관련한 정보는 링크를 참조하시기 바랍니다.
[강릉원주대 대기환경과학과] 2015년 2학기 전필 고급기상통계학 소개 및 과제물
정보 업무명 : 고급기상통계학 작성자 : 이상호 작성일 : 2019-12-20 설 명 : 수정이력 : 내용 [특징] 고급기상통계학 수업에 대한 이해를 돕기위해 작성 [기능] 소개 주제별 과제물 자료의 특성 자료의 상관성..
shlee1990.tistory.com
-
따라서 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)
[전체]
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #====================================================================================================================================================== | |
| # Routine : Main R program | |
| # | |
| # Purpose : Visualization and Principal Components Analysis Using Geopotential Height From ECMWF Interim Numerical Prediction Model | |
| # | |
| # Author : MS. Sang-Ho Lee | |
| # | |
| # Revisions: V1.0 March 16, 2020 First release (MS. Sang-Ho Lee) | |
| #====================================================================================================================================================== | |
| # 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(colorRamps) | |
| library(RColorBrewer) | |
| library(synoptReg) | |
| library(raster) | |
| sFileDirName = Sys.glob("INPUT/interim_daily_06uvz_2009_2018_1.5/*.nc") | |
| # iCount = 1 | |
| dfDataL1 = NA | |
| for (iCount in 1:length(sFileDirName)) { | |
| oOpenNc = RNetCDF::open.nc(sFileDirName[iCount]) | |
| # NetCDF Meta Information | |
| # RNetCDF::print.nc(oOpenNc) | |
| # NetCDF File Read | |
| liReadNc = RNetCDF::read.nc(oOpenNc) | |
| # 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") | |
| # 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) | |
| # 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) | |
| # 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) | |
| # 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) | |
| } | |
| cat(sFileDirName[iCount], " ", dim(dfData), " ", dim(dfDataL1), "\n") | |
| } | |
| dplyr::tbl_df(dfDataL1) | |
| # 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) | |
| # L3 Processing Using Data Frame L2 | |
| dfDataL3 = dfDataL2 %>% | |
| dplyr::filter( | |
| dplyr::between(nLon, 100, 140) | |
| , dplyr::between(nLat, 20, 50) | |
| ) | |
| dplyr::tbl_df(dfDataL3) | |
| # L4 Processing Using Data Frame L3 | |
| dfDataL4 = dfDataL3[ , -c(1:2)] | |
| dplyr::tbl_df(dfDataL4) | |
| # L5 Processing Using Data Frame L4 | |
| dfDataL5 = t(dfDataL4) %>% | |
| as.data.frame() | |
| dplyr::tbl_df(dfDataL5) | |
| #============================================================================= | |
| # Principal Component Analysis | |
| #============================================================================= | |
| liPcaDec = synoptReg::pca_decision(dfDataL5) | |
| liPcaDec | |
| summary(liPcaDec) | |
| # 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) | |
| # Eigen Vector | |
| nEigVec = liPcaDec$PCA$loadings %>% | |
| as.data.frame.matrix() | |
| dfDataL6 = data.frame( | |
| nLon = dfDataL3$nLon | |
| , nLat = dfDataL3$nLat | |
| , nEigVec | |
| ) | |
| dplyr::tbl_df(dfDataL6) | |
| # 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) | |
| # 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) | |
| } | |
| #============================================================================= | |
| # Geopotential Height 150 hPa | |
| #============================================================================= | |
| 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) | |
| 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) | |
| 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 = "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 = 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) | |
| # 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) | |
| #============================================================================= | |
| # Principal Component Analysis Scores | |
| #============================================================================= | |
| 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) | |
| # 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
'프로그래밍 언어 > R' 카테고리의 다른 글
| [R] R을 이용한 통계 분석 및 데이터 시각화 : 데이터형 (0) | 2020.03.23 |
|---|---|
| [R] R를 이용한 고급 기상 통계학 실습 : R 기초, 자료의 특성, 자료의 상관성, 통계적 유의수준, 조화 분석, 주성분 분석, 기계 학습 (0) | 2020.03.17 |
| [R] "rsoi" 패키지를 통해 남방 진동 지수 (Southern Oscillation Index) 자료 다운로드 및 가시화 (0) | 2020.03.08 |
| [R] "rsoi" 패키지를 통해 엘니뇨-남방진동 (EI Nino-Southern Oscillation) 자료 다운로드 및 가시화 (0) | 2020.03.08 |
| [R] 육지 및 해양에 대한 마스킹 (Masking) 처리 및 가시화 (0) | 2020.03.08 |