[R] 직달 일사계 (CHP1, MS56, DR02, GWNU) 비교 관측 자료를 이용한 감도정수 보정 및 시계열 가시화

정보

업무명 : 직달 일사계 (CHP1, MS56, DR02, GWNU) 비교 관측 자료를 이용한 감도정수 보정 및 시계열 가시화
작성자 : 이상호
작성일 : 2020-03-01
설 명 :
수정이력 :

내용

[개요]

안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
오늘은 직달 일사계 비교 관측을 통해 감도 정수 보정 및 시계열 가시화에 대한 내용을 다루고자 합니다.

[특징]

직달 일사계 비교 관측 분석을 위해서 감도정수 보정 및 가시화가 요구되며 이 프로그램은 이러한 목적을 달성하기 위한 소프트웨어

[기능]

직달 일사계 감도 정수 보정
직달 일사량 시계열 가시화

[활용 자료]

자료명 : CR3000_SolysGD_Solar_Min2.dat
자료 종류 : 직달 일사계 (CHP1, MS56, DR02, GWNU)에 대한 날짜 및 시간, mV, 일사량 등
확장자 : txt
기간 2016년 09월 29일 - 2016년 10월 01일
시간 해상도 : 1분 평균
제공처 : 복사-위성 연구소

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

직달 일사량과 감도 정수의 변환
- 직달 일사량은 직달 일사계에 측정된 전압 (mV) 및 장비 내 감도 정수를 통해 환산됩니다.
- 이러한 감도 정수는 생산 당시에 장비마다 고유값을 의미합니다.

감도 정수 보정
- 시간에 따라 센서의 반응도가 떨어지기 때문에 감도 정수 보정은 필연입니다.
- 따라서 실측치 (기준 일사량)와 예측치 (비교 일사량 * 상수)를 통해 최소 오차가 되는 감도 정수 (Sensitivity Constant)를 계산합니다.
- 여기서 실측치 및 비교 일사량은 각각 CHP1의 직달 일사량 및 GWNU의 직달 일사량을 의미합니다.
- 자세한 소스 코드는 하기 "감도 정수 보정" 참조하시기 바랍니다.

[사용법]

입력 자료를 동일 디렉터리 위치
소스 코드를 실행
(Rscript Time_Series_Visualization_and_Sensitivity_Constant_Correction_Using_Direct_Pyrheliometer_Data.R)

가시화 결과를 확인

[사용 OS]

Windows 10

[사용 언어]

R v3.6.2
R Studio v1.2.5033

소스 코드

[명세]

전역 설정
- 최대 10 자리 설정
- 메모리 해제

# Set Option
options(digits = 10)
memory.limit(size = 9999999999999)

라이브러리 읽기

# Library Load
library(tidyverse)
library(data.table)
library(lubridate)
library(timeDate)
library(extrafont)
library(colorRamps)
library(RColorBrewer)
library(tidyverse)
library(data.table)
library(Metrics)
library(ggplot2)
library(scales)
library(hash)

함수 정의
- Map 형태인 전역 변수를 설정 (fnSetMap) 및 가져오기 (fnGetMap)를 위한 유용한 함수
- 실측치 (nActual) 그리고 예측치 (nPredicted) * 상수값 (nFactor)가 최소 오차인 상수값을 구하는 함수 (fnGetCalibFactor)

# Function Load
fnSetMap = function(maData, sKey, sValue) {
   
   for (iCount in 1:length(sKey)) {
      # check whether the key exists
      isKey = has.key(sKey[iCount], maData)
      
      if (isKey == TRUE) { 
         cat(sKey[iCount], "is Duplicated", "\n")
      }
      
      maData[sKey[iCount]] = sValue[iCount]
   }
}

fnGetMap =  function(maData, sKey) {
   
   # check whether the key exists
   isKey = has.key(sKey, maData)
   
   if (isKey == FALSE) {
      cat("Key is not exist")
   }
   
   return (maData[[sKey]])
}


fnGetCalibFactor = function(nActual, nPredicted, nMin, nMax, nInterval, isPlot = FALSE) {
   
   nFactor = seq(nMin, nMax, by = nInterval)
   
   # RMSE Fitting
   liResultTmp = lapply(1:length(nFactor), function(iCount) Metrics::rmse(nActual, nPredicted * nFactor[iCount]))   
   
   liResult = unlist(liResultTmp)
   
   if (isPlot == TRUE) {
      plot(liResult)   
   }
   
   # Best Factor Index
   iIndex = which(liResult == min(liResult, na.rm = TRUE))

   return (nFactor[[iIndex]])
}

전역 변수 초기화
- 각 직달 일사계마다 감도 정수를 설정

# Global Variable Init
maData = hash()

sKey = c("CHP1_SENS_CONS", "MS56_SENS_CONS", "DR02_SENS_CONS", "GWNU_SENS_CONS")
sValue = c(8.49, 8.531, 10.59, 17.8928)

fnSetMap(maData, sKey, sValue)

파일 읽기

# File Read
dfCr3000 = data.table::fread("CR3000_SolysGD_Solar_Min2.dat", header = FALSE, sep = ",", skip = 4)

dplyr::tbl_df(dfCr3000)

Data Frame 설정
- 날짜 및 시간 정보를 추가

# Set Data Frame
dfData = dfCr3000 %>%
   dplyr::mutate(
      dtDateTime = readr::parse_datetime(V1, format = "%Y-%m-%d %H:%M:%S")
      , nYear = lubridate::year(dtDateTime)
      , nMonth = lubridate::month(dtDateTime)
      , nDay = lubridate::day(dtDateTime)
      , nHour = lubridate::hour(dtDateTime)
      , nMinute = lubridate::minute(dtDateTime)
      , nSec = lubridate::second(dtDateTime)
      , nLastDayInMonth = lubridate::day(timeDate::timeLastDayInMonth(format(dtDateTime, "%Y%m%d"), format = "%Y%m%d", zone = "Asia/Seoul"))
      , xranYmdHms = 
         nYear + ((nMonth - 1) / 12.0) + ((nDay - 1) / (12.0 * nLastDayInMonth))  + (nHour / (12.0 * nLastDayInMonth * 24.0))
      + (nMinute / (12.0 * nLastDayInMonth * 24.0 * 60.0))
      + (nSec / (12.0 * nLastDayInMonth * 24.0 * 60.0 * 60.0))
      , xranYmd = nYear + ((nMonth - 1) / 12.0) + ((nDay - 1) / (12.0 * nLastDayInMonth))
      , xranHms = nHour + (nMinute / 60.0) + (nSec / (60.0 * 60.0))
   )

colnames(dfData) = c("datetime", "count", "nRadCHP1", "nRadMS56", "nRadDR02", "CHP1_temp", "MS56_temp", "DR02_temp", "Vel_CMP22", "Vel_CMP22_temp", "nMvGWNU", "GWNU2_mv", "dtDateTime", "nYear", "nMonth", "nDay", "nHour", "nMinute", "nSec", "nLastDayInMonth", "xranYmdHms", "xranYmd", "xranHms")

dplyr::glimpse(dfData)

Data Frame를 이용한 L1 전처리
- "상기 직달 일사량과 감도 정수의 변환"에서와 같이 전압 (mV) 및 감도 정수를 변환
- 날짜 및 시간 설정 (10월 01일 08-17시)

# L1 Processing Using Data Frame
dfDataL1 = dfData %>%
   dplyr::mutate(
      nMvCHP1 = nRadCHP1 * (fnGetMap(maData, "CHP1_SENS_CONS") / 1000.0)
      , nMvMS56 = nRadMS56 * (fnGetMap(maData, "DR02_SENS_CONS") / 1000.0)
      , nMvDR02 = nRadDR02 * (fnGetMap(maData, "DR02_SENS_CONS") / 1000.0) 
      , nRadGWNU = nMvGWNU * (1000.0 / fnGetMap(maData, "GWNU_SENS_CONS"))
      ) %>%
   dplyr::filter(
      nMonth == 10
      , nDay == 01
      , dplyr::between(xranHms, 08, 17)
   )

dplyr::glimpse(dfDataL1)

전압 시계열을 위한 초기 설정

# Set Value for Visualization 
xAxis = dfDataL1$dtDateTime
yAxis = dfDataL1$nMvGWNU
yAxis2 = dfDataL1$nMvCHP1
yAxis3 = dfDataL1$nMvMS56
yAxis4 = dfDataL1$nMvDR02

sDateYmd = format(xAxis, "%Y-%m-%d") %>% first()
nMaxXran = max(xAxis, na.rm = TRUE)
nMinXran = min(xAxis, na.rm = TRUE)
font = "Palatino Linotype"

plot을 이용한 전압 시계열

# mV Times Series Using plot
png(file = paste0("FIG/mV_", sDateYmd, ".png"), 1200, 600)
par(mar = c(5.1, 7, 5.1, 5.1), cex.main = 2.0, cex.lab = 1.6, cex.axis = 1.6, cex = 1.3, font.axis = 2, font.lab = 2, font.main = 2, font = 2, family = font)

plot(
   xAxis, yAxis
   , xlab = "", ylab = ""
   , yaxs = "i", xaxs = "i"
   , xaxt = "n", yaxt = "n"
   , main = paste0(sDateYmd, " : mV")
   , pch = 21, col = "black", bg = "black", cex = 0.8, type = "l"
   , xlim = c(nMinXran, nMaxXran)
   , ylim = c(0, 7)
)

axis.POSIXct(1, xAxis, format="%H:%M", at = seq(nMinXran, nMaxXran, by = "1 hour"), las = 1 )
mtext("Time  [Hour : Minute]", side = 1, col = "black", line = 3, cex = 2.2)

axis(2, las = 1, at = seq(0, 7, 1))
mtext(expression(paste(bold("Sensitivity  Constant"))), side = 2,  line = 3, cex = 2.2)

points(xAxis, yAxis2, col = "red", pch = 21, bg = "black", cex = 0.8, type = "l")
points(xAxis, yAxis3, col = "blue", pch = 21, bg = "red", cex = 0.8, type = "l")
points(xAxis, yAxis4, col = "green", pch = 21, bg = "blue", cex = 0.8, type = "l")

legend(
   "topright"
   , legend = c("GWNU", "CHP1", "MS56", "DR02")
   , col = c("black", "red", "blue", "green")
   , cex = 1.4, text.font = 2, bty = "n", xpd = TRUE 
   , pt.bg = c("black", "red", "blue", "green")
   , text.col = c("black", "red", "blue", "green")
   , lwd = 2
   )

dev.off()

전압 차이 시계열을 위한 초기 설정

# Set Value for Visualization
nMvCHP1 = dfDataL1$nMvCHP1
nMvGWNU = dfDataL1$nMvGWNU
nMvMS56 = dfDataL1$nMvMS56
nMvDR02 = dfDataL1$nMvDR02

xAxis = dfDataL1$dtDateTime
yAxis = nMvCHP1 - nMvGWNU
yAxis2 = nMvCHP1 - nMvMS56
yAxis3 = nMvCHP1 - nMvDR02

sDateYmd = format(xAxis, "%Y-%m-%d") %>% first()
nMaxXran = max(xAxis, na.rm = TRUE)
nMinXran = min(xAxis, na.rm = TRUE)
font = "Palatino Linotype"

plot을 이용한 전압 차이 시계열

# Difference mV Times Series Using plot
png(file = paste0("FIG/Diff_mV_", sDateYmd, ".png"), 1200, 600)
par(mar = c(5.1, 7, 5.1, 5.1), cex.main = 2.0, cex.lab = 1.6, cex.axis = 1.6, cex = 1.3, font.axis = 2, font.lab = 2, font.main = 2, font = 2, family = font)

plot(
   xAxis, yAxis
   , xlab = "", ylab = ""
   , yaxs = "i", xaxs = "i"
   , xaxt = "n", yaxt = "n"
   , main = paste0(sDateYmd, " : mV Difference")
   , pch = 21, col = "black", bg = "black", cex = 0.8, type = "l"
   , xlim = c(nMinXran, nMaxXran)
   , ylim = c(-4, 4)
)

axis.POSIXct(1, xAxis, format="%H:%M", at = seq(nMinXran, nMaxXran, by = "1 hour"), las = 1 )
mtext("Time  [Hour : Minute]", side = 1, col = "black", line = 3, cex = 2.2)

axis(2, las = 1, at = seq(-4, 4, 1))
mtext(expression(paste(bold("Voltage  [mV]"))), side = 2,  line = 4, cex = 2.2)

abline(h = 0, col = "green", font = 2, lw = 2)
          
points(xAxis, yAxis2, col = "red", pch = 21, bg = "red", cex = 0.8, type = "l")
points(xAxis, yAxis3, col = "blue", pch = 21, bg = "blue", cex = 0.8, type = "l")

legend(
   "topright"
   , legend = c(
      paste0("CHP1 - GWNU (RMSE = ", round(Metrics::rmse(nMvCHP1, nMvGWNU), 2), ")")
      , paste0("CHP1 - MS56 (RMSE = ", round(Metrics::rmse(nMvCHP1, nMvMS56), 2), ")")
      , paste0("CHP1 - DR02 (RMSE = ", round(Metrics::rmse(nMvCHP1, nMvDR02), 2), ")")
      )
   , col = c("black", "red", "blue")
   , cex = 1.4, text.font = 2, bty = "n", xpd = TRUE
   , pt.bg = c("black", "red", "blue")
   , text.col = c("black", "red", "blue")
   , lwd = 2
)

dev.off()

일사량 시계열을 위한 초기 설정

# Set Value for Visualization
xAxis = dfDataL1$dtDateTime
yAxis = dfDataL1$nRadGWNU
yAxis2 = dfDataL1$nRadCHP1
yAxis3 = dfDataL1$nRadMS56
yAxis4 = dfDataL1$nRadDR02

sDateYmd = format(xAxis, "%Y-%m-%d") %>% first()
nMaxXran = max(xAxis, na.rm = TRUE)
nMinXran = min(xAxis, na.rm = TRUE)
font = "Palatino Linotype"

plot을 이용한 일사량 시계열

# Difference Direct Shortwave Radiation Times Series Using plot
png(file = paste0("FIG/Radion_", sDateYmd, ".png"), 1200, 600)
par(mar = c(5.1, 7, 5.1, 5.1), cex.main = 2.0, cex.lab = 1.6, cex.axis = 1.6, cex = 1.3, font.axis = 2, font.lab = 2, font.main = 2, font = 2, family = font)

plot(
   xAxis, yAxis
   , xlab = "", ylab = ""
   , yaxs = "i", xaxs = "i"
   , xaxt = "n", yaxt = "n"
   , main = paste0(sDateYmd, " : Direct Shortwave Radiation")
   , pch = 21, col = "black", bg = "black", cex = 0.8, type = "l"
   , xlim = c(nMinXran, nMaxXran)
   , ylim = c(0, 400)
)

axis.POSIXct(1, xAxis, format="%H:%M", at = seq(nMinXran, nMaxXran, by = "1 hour"), las = 1 )
mtext("Time  [Hour : Minute]", side = 1, col = "black", line = 3, cex = 2.2)

axis(2, las = 1, at = seq(0, 400, 50))
mtext(expression(paste(bold("Direct  Shortwave  Radiation  [Wm"^bold("-2")*"]"))), side = 2,  line = 4, cex = 2.2)

points(xAxis, yAxis2, col = "red", pch = 21, bg = "red", cex = 0.8, type = "l")
points(xAxis, yAxis3, col = "blue", pch = 21, bg = "blue", cex = 0.8, type = "l")
points(xAxis, yAxis4, col = "green", pch = 21, bg = "green", cex = 0.8, type = "l")

legend(
   "topright"
   , legend = c("GWNU", "CHP1", "MS56", "DR02")
   , col = c("black", "red", "blue", "green")
   , cex = 1.4, text.font = 2, bty = "n", xpd = TRUE 
   , pt.bg = c("black", "red", "blue", "green")
   , text.col = c("black", "red", "blue", "green")
   , lwd = 2
)
dev.off()

일사량 차이 시계열을 위한 초기 설정

# Set Value for Visualization
nRadCHP1 = dfDataL1$nRadCHP1
nRadGWNU = dfDataL1$nRadGWNU
nRadMS56 = dfDataL1$nRadMS56
nRadDR02 = dfDataL1$nRadDR02

xAxis = dfDataL1$dtDateTime
yAxis = nRadCHP1 - nRadGWNU
yAxis2 = nRadCHP1 - nRadMS56
yAxis3 = nRadCHP1 - nRadDR02

sDateYmd = format(xAxis, "%Y-%m-%d") %>% first()
nMaxXran = max(xAxis, na.rm = TRUE)
nMinXran = min(xAxis, na.rm = TRUE)
font = "Palatino Linotype"

plot을 이용한 일사량 차이 시계열

# Difference Direct Shortwave Radiation Times Series Using plot
png(file = paste0("FIG/Diff_Rad_", sDateYmd, ".png"), 1200, 600)
par(mar = c(5.1, 7, 5.1, 5.1), cex.main = 2.0, cex.lab = 1.6, cex.axis = 1.6, cex = 1.3, font.axis = 2, font.lab = 2, font.main = 2, font = 2, family = font)

plot(
   xAxis, yAxis
   , xlab = "", ylab = ""
   , yaxs = "i", xaxs = "i"
   , xaxt = "n", yaxt = "n"
   , main = paste0(sDateYmd, " : Direct Shortwave Radiation Difference")
   , pch = 21, col = "black", bg = "black", cex = 0.8, type = "l"
   , xlim = c(nMinXran, nMaxXran)
   , ylim = c(-20, 20)
)

axis.POSIXct(1, xAxis, format="%H:%M", at = seq(nMinXran, nMaxXran, by = "1 hour"), las = 1 )
mtext("Time  [Hour : Minute]", side = 1, col = "black", line = 3, cex = 2.2)

axis(2, las = 1, at = seq(-20, 20, 5))
mtext(expression(paste(bold("Direct  Shortwave Radiation  [Wm"^bold("-2")*"]"))), side = 2,  line = 4, cex = 2.2)

points(xAxis, yAxis2, col = "red", pch = 21, bg = "red", cex = 0.8, type = "l")
points(xAxis, yAxis3, col = "blue", pch = 21, bg = "blue", cex = 0.8, type = "l")

abline(h = 0, col = "green", font = 2, lw = 2)

legend(
   "topright"
   , legend = c(
      paste0("CHP1 - GWNU (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadGWNU), 2), ")")
      , paste0("CHP1 - MS56 (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadMS56), 2), ")")
      , paste0("CHP1 - DR02 (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadDR02), 2), ")")
   )
   , col = c("black", "red", "blue")
   , cex = 1.4, text.font = 2, bty = "n", xpd = TRUE
   , pt.bg = c("black", "red", "blue")
   , text.col = c("black", "red", "blue")
   , lwd = 2
)

dev.off()

감도 정수 보정
- 이론적 배경은 상기 감도 정수 보정에서 참조
- GWNU 직달 일사계의 감도 정수는 17.89 > 18.48로 보정

# Sensitivity Constant Correction
nGwnuSensCons = fnGetMap(maData, "GWNU_SENS_CONS")

nActual = dfDataL1$nRadCHP1
nPredicted = dfDataL1$nRadGWNU

nCalibFactor = fnGetCalibFactor(nActual, nPredicted, -100, 100, 0.001, FALSE)
nGwnuSensConsNew = nGwnuSensCons / nCalibFactor

dfDataL2 = dfDataL1 %>%
   dplyr::mutate(
      nRadNewGWNU = nMvGWNU * (1000.0 / nGwnuSensConsNew)
   )

cat(
   sprintf("Old Sensitivity Constant : %.2f", nGwnuSensCons)
   , sprintf("\nNew Sensitivity Constant : %.2f", nGwnuSensConsNew)
)

dplyr::glimpse(dfDataL2)

plot을 이용한 일사량 시계열 (감도 정수 보정)

# Set Value for Visualization
xAxis = dfDataL2$dtDateTime
yAxis = dfDataL2$nRadNewGWNU
yAxis2 = dfDataL2$nRadCHP1
yAxis3 = dfDataL2$nRadMS56
yAxis4 = dfDataL2$nRadDR02

sDateYmd = format(xAxis, "%Y-%m-%d") %>% first()
nMaxXran = max(xAxis, na.rm = TRUE)
nMinXran = min(xAxis, na.rm = TRUE)
font = "Palatino Linotype"

# New Direct Shortwave Radiation Times Series Using plot
png(file = paste0("FIG/NEW_Radion_", sDateYmd, ".png"), 1200, 600)
par(mar = c(5.1, 7, 5.1, 5.1), cex.main = 2.0, cex.lab = 1.6, cex.axis = 1.6, cex = 1.3, font.axis = 2, font.lab = 2, font.main = 2, font = 2, family = font)

plot(
   xAxis, yAxis
   , xlab = "", ylab = ""
   , yaxs = "i", xaxs = "i"
   , xaxt = "n", yaxt = "n"
   , main = paste0(sDateYmd, " : Direct Shortwave Radiation (Sensitivity Constant Correction)")
   , pch = 21, col = "black", bg = "black", cex = 0.8, type = "l"
   , xlim = c(nMinXran, nMaxXran)
   , ylim = c(0, 400)
)

axis.POSIXct(1, xAxis, format="%H:%M", at = seq(nMinXran, nMaxXran, by = "1 hour"), las = 1 )
mtext("Time  [Hour : Minute]", side = 1, col = "black", line = 3, cex = 2.2)

axis(2, las = 1, at = seq(0, 400, 50))
mtext(expression(paste(bold("Direct  Shortwave  Radiation  [Wm"^bold("-2")*"]"))), side = 2,  line = 4, cex = 2.2)

points(xAxis, yAxis2, col = "red", pch = 21, bg = "red", cex = 0.8, type = "l")
points(xAxis, yAxis3, col = "blue", pch = 21, bg = "blue", cex = 0.8, type = "l")
points(xAxis, yAxis4, col = "green", pch = 21, bg = "green", cex = 0.8, type = "l")

legend(
   "topright"
   , legend = c("GWNU", "CHP1", "MS56", "DR02")
   , col = c("black", "red", "blue", "green")
   , cex = 1.4, text.font = 2, bty = "n", xpd = TRUE 
   , pt.bg = c("black", "red", "blue", "green")
   , text.col = c("black", "red", "blue", "green")
   , lwd = 2
)
dev.off()

plot을 이용한 일사량 차이 시계열 (감도 정수 보정)

# Set Value for Visualization
nRadCHP1 = dfDataL2$nRadCHP1
nRadGWNU = dfDataL2$nRadNewGWNU
nRadMS56 = dfDataL2$nRadMS56
nRadDR02 = dfDataL2$nRadDR02

xAxis = dfDataL2$dtDateTime
yAxis = nRadCHP1 - nRadGWNU
yAxis2 = nRadCHP1 - nRadMS56
yAxis3 = nRadCHP1 - nRadDR02

sDateYmd = format(xAxis, "%Y-%m-%d") %>% first()
nMaxXran = max(xAxis, na.rm = TRUE)
nMinXran = min(xAxis, na.rm = TRUE)
font = "Palatino Linotype"

# Difference New Direct Shortwave Radiation Times Series Using plot
png(file = paste0("FIG/Diff_NEW_Rad_", sDateYmd, ".png"), 1200, 600)
par(mar = c(5.1, 7, 5.1, 5.1), cex.main = 2.0, cex.lab = 1.6, cex.axis = 1.6, cex = 1.3, font.axis = 2, font.lab = 2, font.main = 2, font = 2, family = font)

plot(
   xAxis, yAxis
   , xlab = "", ylab = ""
   , yaxs = "i", xaxs = "i"
   , xaxt = "n", yaxt = "n"
   , main = paste0(sDateYmd, " : Difference Direct Shortwave Radiation (Sensitivity Constant Correction)")
   , pch = 21, col = "black", bg = "black", cex = 0.8, type = "l"
   , xlim = c(nMinXran, nMaxXran)
   , ylim = c(-20, 20)
)

axis.POSIXct(1, xAxis, format="%H:%M", at = seq(nMinXran, nMaxXran, by = "1 hour"), las = 1 )
mtext("Time  [Hour : Minute]", side = 1, col = "black", line = 3, cex = 2.2)

axis(2, las = 1, at = seq(-20, 20, 5))
mtext(expression(paste(bold("Direct  Shortwave Radiation  [Wm"^bold("-2")*"]"))), side = 2,  line = 4, cex = 2.2)

points(xAxis, yAxis2, col = "red", pch = 21, bg = "red", cex = 0.8, type = "l")
points(xAxis, yAxis3, col = "blue", pch = 21, bg = "blue", cex = 0.8, type = "l")

abline(h = 0, col = "green", font = 2, lw = 2)

legend(
   "topright"
   , legend = c(
      paste0("CHP1 - GWNU (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadGWNU), 2), ")")
      , paste0("CHP1 - MS56 (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadMS56), 2), ")")
      , paste0("CHP1 - DR02 (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadDR02), 2), ")")
   )
   , col = c("black", "red", "blue")
   , cex = 1.4, text.font = 2, bty = "n", xpd = TRUE
   , pt.bg = c("black", "red", "blue")
   , text.col = c("black", "red", "blue")
   , lwd = 2
)

dev.off()

GWNU 직달 일사계에 대해 감도 정수 보정한 결과 5.29 > 4.93으로 개선되었습니다.
그럼에도 불구하고 2016년 10월 1일 (흐린 날)에 GWNU 직달 일사계는 타 일사계 (MS56, DR02)보다 다소 오차가 발생하였습니다.
이는 특정 사례에 대한 결과이기 때문에 장기간 비교 관측을 통한 통계 결과가 필요할 것으로 사료됩니다.

ggplot2를 이용한 가시화

최근 데이터 분석에서 자주 사용하는 ggplot2를 통해 가시화 해보았습니다.
일반적으로 연구에서는 기본 plot을 주로 이용하나 학술 논문 및 발표에서는 ggplot2를 사용하고 합니다.
이는 효과적인 전달하는 데 도와줍니다.

ggplot2 가시화를 위한 초기 설정
- 가시화에 필요한 컬럼 추출 및 차이 시계열을 위한 데이터 정제 (dfDataL3)
- tidyr::gather를 통해 key 및 value로 데이터 정제 (dfDataL4)
- dplyr::filter를 통해 시계열 일사량 (dfDataL5)과 차이 일사량 (dfDataL6)에 대해 데이터 정제

# Set Value for Visualization
nRadCHP1 = dfDataL2$nRadCHP1
nRadGWNU = dfDataL2$nRadNewGWNU
nRadMS56 = dfDataL2$nRadMS56
nRadDR02 = dfDataL2$nRadDR02

dfDataL3 = dfDataL2 %>%
   dplyr::select(dtDateTime, nRadCHP1, nRadMS56, nRadDR02, nRadNewGWNU) %>%
   dplyr::mutate(
      nDiffGWNU = nRadCHP1 - nRadNewGWNU
      , nDiffMS56 = nRadCHP1 - nRadMS56
      , nDiffDR02 = nRadCHP1 - nRadDR02
   )

dplyr::tbl_df(dfDataL3)

dfDataL4 = dfDataL3 %>%
   tidyr::gather(key = "key", value = "value", -dtDateTime)
   
dplyr::tbl_df(dfDataL4)

dfDataL5 = dfDataL4 %>%
   dplyr::filter(! str_detect(key, pattern = "nDiff*"))
 
dplyr::tbl_df(dfDataL5)

dfDataL6 = dfDataL4 %>%
   dplyr::filter(str_detect(key, pattern = "nDiff*"))

dplyr::tbl_df(dfDataL6)

ggplot2을 이용한 일사량 시계열 (감도 정수 보정)

# New Direct Shortwave Radiation Times Series Using ggplot2
dfDataL5 %>%  
   ggplot(aes(x = dtDateTime, y = value, colour = key)) +
   geom_line() +
   theme_bw() +
   scale_x_datetime(breaks = date_breaks("1 hours"), labels = date_format("%H:%M", tz="UTC"), expand=c(0,0)) +
   scale_y_continuous(breaks=seq(0, 400, by = 50), expand=c(0,0), limits=c(0, 400)) +
   scale_color_discrete(
      breaks = c("nRadCHP1", "nRadDR02", "nRadMS56", "nRadNewGWNU")
      , labels = c("CHP1", "DR02", "MS56", "GWNU")
      ) +
   labs(
      x = "Time  [Hour : Minute]"
      , y = expression(paste(bold("Direct  Shortwave  Radiation  [Wm"^bold("-2")*"]")))
      , fill = ""
      , colour = ""
      , title = ""
   ) +
   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.position = c(1, 1.03)
      , legend.justification = c(1, 1)
      , 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/ggplot2_Figure1.png"), width = 12, height = 8, dpi = 600)

ggplot2을 이용한 차이 일사량 시계열 (감도 정수 보정)

# Difference New Direct Shortwave Radiation Times Series Using ggplot2
dfDataL6 %>% 
   ggplot(aes(x = dtDateTime, y = value, colour = key)) +
   geom_line() +
   theme_bw() +
   scale_x_datetime(breaks = date_breaks("1 hours"), labels = date_format("%H:%M", tz="UTC"), expand=c(0,0)) +
   scale_y_continuous(breaks=seq(-20, 20, by = 5), expand=c(0,0), limits=c(-20, 20)) +
   scale_color_discrete(
      breaks = c("nDiffGWNU", "nDiffMS56", "nDiffDR02")
      , labels = c(
         paste0("CHP1 - GWNU (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadGWNU), 2), ")")
         , paste0("CHP1 - MS56 (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadMS56), 2), ")")
         , paste0("CHP1 - DR02 (RMSE = ", round(Metrics::rmse(nRadCHP1, nRadDR02), 2), ")")
         )) +
   labs(
      x = "Time  [Hour : Minute]"
      , y = expression(paste(bold("Direct  Shortwave  Radiation  [Wm"^bold("-2")*"]")))
      , fill = ""
      , colour = ""
      , title = ""
   ) +
   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.position = c(1, 1.03)
      , legend.justification = c(1, 1)
      , 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/ggplot2_Figure2.png"), width = 12, height = 8, dpi = 600)

[전체]

참고 문헌

[논문]

없음

[보고서]

없음

[URL]

없음

문의사항

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

sangho.lee.1990@gmail.com

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

saimang0804@gmail.com

저작자표시 비영리 변경금지

'프로그래밍 언어 > R' 카테고리의 다른 글

[R] 라디오존데 (Radiosonde) 관측 자료를 이용한 고도별 온도 (Temperature) 및 상대 습도 (Relative Humidity) 등고선 가시화 (0)	2020.03.03
[R] 데이터의 분포를 그래프로 확인하는 "ggridges"패키지 (0)	2020.03.02
[R] 하와이 마우나로아 (Manuna Loa)에서 이산화탄소 (CO2) 농도 자료를 이용한 통계 처리 및 시계열 가시화 (0)	2020.03.01
[R] nonogram (네모네모로직) 해결 알고리즘 및 가시화 (0)	2020.02.23
[R] 파일 관련 기본 명령어 (디렉터리/파일 선택, 생성, 삭제, 복사) 소개 (0)	2020.02.13

정보

내용

[개요]

[특징]

[기능]

[활용 자료]

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

[사용법]

[사용 OS]

[사용 언어]

소스 코드

[명세]

[전체]

참고 문헌

[논문]

[보고서]

[URL]

문의사항

'프로그래밍 언어 > R' 카테고리의 다른 글

티스토리툴바