[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