[R] 새해 경자년에 해돋이/일출/해 뜨는 시간 가시화

 정보

• 업무명     : 소프트웨어 개발

• 작성자     : 이상호

• 작성일     : 2019-12-31

• 설   명      :

• 수정이력 :

 

 내용

[특징]

• 새해 경자년에 해돋이 (일출) 시간을 이해하기 위해서 가시화가 요구되며 이 프로그램은 이러한 목적을 달성하기 위한 소프트웨어

 

[기능]

• 기상 관측소 정보를 이용한 전국 해돋이 (일출) 시각 맵핑

• 각 관측소에 대한 2019년 일출/일몰 시간 시계열

 

[활용 자료]

• 자료명 : 기상 관측소 정보

• 자료 종류 : 관측소 지점 번호, 경도, 위도, 관측소명

 

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

• 없음

 

[사용법]

• 입력 자료를 동일 디렉터리 위치
• 소스 코드를 실행 (Rscript Visualization_and_Preprocessing_Using_Weather_Station_Data_for_Sunrise_and_Sunset_Times.R)
• 가시화 결과 확인

 

[사용 OS]

• Window 10

 

[사용 언어]

• R v 3.6.2

• R Studio v1.2.5033

 

 소스 코드

[명세] 해돋이 (일출) 시각 맵핑 

• 전역 설정

  • 자리수 설정

  • 영문 설정

# Set Option
options(digits = 10)
Sys.setlocale("LC_TIME", "english") 

 

• 라이브러리 읽기

# Library Load
library(data.table)
library(suncalc) 
library(tidyverse)
library(scales)
library(sf)
library(gstat)
library(sp)
library(metR)
library(colorRamps)
library(ggrepel)
library(extrafont)
library(timeDate)

 

• 기상 관측소 정보 읽기

# File Read
dfStation = fread("INPUT/Station_Information.dat", sep = "\t", header = FALSE, skip = 1)
colnames(dfStation) = c("station", "lon", "lat", "stationName") 

dplyr::tbl_df(dfStation) 

 

 

• Data Frame 형태로 초기화

# Set Data Frame
dfData = dfStation %>%
   dplyr::mutate(date = lubridate::ymd("2019-01-01"))

dplyr::tbl_df(dfData) 

 

 

• Data Frame를 통해 L1 전처리

  • getSunlightTimes 함수를 통해 일출 (해돋이) 및 일몰 시간 계산

  • dfDataL1을 기준으로 기상 관측소 정보를 좌측 조인

  • xranYmdHms, xranHms 변수 초기화

# L1 Processing Using Data Frame
dfDataL1 = getSunlightTimes(
   data = dfData
   , keep = c("sunrise", "sunriseEnd", "sunset", "sunsetStart")
   ,  tz = "Asia/Seoul"
) %>%
   dplyr::left_join(dfStation, by=c("lat" = "lat", "lon" = "lon")) %>%
   dplyr::mutate(
      dtDate = sunrise
      , nYear = lubridate::year(dtDate)
      , nMonth = lubridate::month(dtDate)
      , nDay = lubridate::day(dtDate)
      , nHour = lubridate::hour(dtDate)
      , nMinute = lubridate::minute(dtDate)
      , nSec = lubridate::second(dtDate)
      , nLastDayInMonth = lubridate::day(timeDate::timeLastDayInMonth(format(dtDate, "%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))
      , xranHms = 
         nHour + (nMinute / 60.0) + (nSec / (60.0 * 60.0))
      )

dplyr::tbl_df(dfDataL1)
dplyr::glimpse(dfDataL1)

 

 

 

• 맵핑을 위한 설정

# Set Value for Visualization
cbMatlab =  colorRamps::matlab.like(11)
font = "Palatino Linotype"
mapKor = read_sf("INPUT/gadm36_KOR_shp/gadm36_KOR_1.shp")
mapPrk = read_sf("INPUT/gadm36_PRK_shp/gadm36_PRK_1.shp")
mapJpn = read_sf("INPUT/gadm36_JPN_shp/gadm36_JPN_1.shp")

spData = dfDataL1
coordinates(spData) = ~lon + lat

# plot(spData)

 

 

• 기상 관측소 정보를 통해 전처리

  • 내삽을 위한 위/경도 설정

  • 신규 격자 (spNewData)에 대해 "역거리 가중치 (IDW)" 공간 내삽 수행

# Processing Using Station Data
dfStationData = data.frame(
   nLat = spData$lat
   , nLon = spData$lon
   , sStationName  = as.character(spData$stationName)
)

# Min/Max Longitude of the Interpolation Area
xRange = as.numeric(c(124, 132))  

# Min/Max Latitude of the Interpolation Area
yRange = as.numeric(c(33, 39)) 

# Expand Points to Grid
spNewData = expand.grid(
   x = seq(from = xRange[1], to = xRange[2], by = 0.1)
   , y = seq(from = yRange[1], to = yRange[2], by = 0.1)
)

coordinates(spNewData) = ~ x + y
gridded(spNewData) = TRUE

# Apply IDW Model for the Data
spDataL1 = gstat::idw(formula = xranHms ~ 1, locations = spData, newdata = spNewData)  

 

• L1을 통해 L2 전처리

# L2 Processing Using L1 Data Frame
dfDataL2 = spDataL1 %>%
   as.data.frame() %>%
   dplyr::rename(
      nLon = x
      , nLat = y
      , nVal = var1.pred
   ) %>%
   dplyr::mutate(isLand = metR::MaskLand(nLon, nLat, mask = "world")) %>%
   dplyr::filter(isLand == TRUE)

dplyr::tbl_df(dfDataL2)

 

 

• L2 자료를 이용한 맵핑

  • 해돋이 (일출) 시각에 대한 등고선

  • 기상 관측소 정보 기입

  • 남한/북한/일본 Shp

# Map Visualization Using ggplot2
ggplot() + 
   coord_fixed(ratio = 1.1) +
   theme_bw() +
   geom_tile(data = dfDataL2, aes(x = nLon, y = nLat, fill = nVal)) +
   geom_point(data = dfStationData, aes(x = nLon, y = nLat), colour = "black", size= 5, alpha = 0.3, show.legend = FALSE) +
   geom_text_repel(data = dfStationData, aes(x = nLon, y = nLat, label = sStationName, colour=sStationName), point.padding = 0.25, box.padding = 0.25,  nudge_y = 0.1, size = 4, colour = "black", family = font) +
   geom_contour2(data = dfDataL2, aes(x = nLon, y = nLat, z = nVal), color = "red", alpha = 0.75, breaks = seq(7.5, 8.0, 0.05), show.legend = TRUE) +
   geom_text_contour(data = dfDataL2, aes(x = nLon, y = nLat, z = nVal), stroke = 0.2, breaks = seq(7.5, 8.0, 0.05), show.legend = TRUE) +
   scale_fill_gradientn(colours = cbMatlab, limits=c(7.5, 8), breaks = seq(7.5, 8.0, 0.25), na.value = cbMatlab[length(cbMatlab)], labels = c("7.50 (07:30)", "7.75 (07:45)" ,"8.00 (08:00)")) + 
   geom_sf(data = mapKor, color = "black", fill = NA) +
   geom_sf(data = mapPrk, color = "black", fill = NA) +
   geom_sf(data = mapJpn, color = "black", fill = NA) +
   metR::scale_x_longitude(expand = c(0, 0), breaks = seq(124, 132, 2), limits = c(124, 132)) +
   metR::scale_y_latitude(expand = c(0, 0), breaks = seq(32, 40, 1), limits = c(33, 39)) +
   labs(x = ""
        , y = ""
        , 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/SunRise.png'), width = 8, height = 10, dpi = 600)

 

그림. 2019년 1월 1일 (경자년)에 대한 일출 (해돋이) 시각 맵핑.

 

[명세] 해돋이 (일출) 시각 시계열

• Data Frame 형태로 초기화

# Set Data Frame
dtDate = seq(lubridate::ymd("2019-01-01"), lubridate::ymd("2020-01-01"), by = "1 day")

for (iCount in 1:length(dtDate)) {
   
   dfTmp = dfStation %>%
      mutate(date = dtDate[iCount])
   
   if (iCount == 1) {
      dfTimeData = tibble::add_column(dfTmp)
   } else {
      dfTimeData = bind_rows(dfTimeData, dfTmp)
   }
}

dplyr::tbl_df(dfTimeData)

 

 

• Data Frame를 통해 L1 전처리

  • getSunlightTimes 함수를 통해 일출 (해돋이) 및 일몰 시간 계산

  • dfDataL1을 기준으로 기상 관측소 정보를 좌측 조인

dfTimeDataL1 = getSunlightTimes(
   data = dfTimeData
   , keep = c("sunrise", "sunriseEnd", "sunset", "sunsetStart")
   , tz = "Asia/Seoul"
   ) %>%
   dplyr::left_join(dfStation, by=c("lat" = "lat", "lon" = "lon")) %>%
   dplyr::mutate(
      dtDate = as.POSIXct(date) - lubridate::hours(9)
      , dtSunrise = sunrise - dtDate
      , dtSunriseEnd = sunriseEnd - dtDate
      , dtSunset =  sunset - dtDate
      , dtSunsetStart = sunsetStart - dtDate
      )

dplyr::glimpse(dfTimeDataL1)

 

 

• 시계열을 위한 설정

  • 폰트 설정

# Set Value for Visualization
font = "Palatino Linotype"

 

• L1 자료를 이용한 시계열

  • 2019년 1월 1일부터 2020년 1월 1일까지 하루 간격으로 해돋이 (일출) 및 일몰 시각에 대한 선 그래프

  • 동/서/남해에 대한 대표 도시 (강릉, 인천, 부산) 선정

# Time series Visualization Using ggplot2
for (iCount in 1:nrow(dfStation)) {
   
   sStationName = as.character(dfStation[iCount, 4])
   sSaveFilePathName = paste0("FIG/Sunrise_and_Sunset_", sStationName, ".png")
   
   dfTimeDataL2 = dfTimeDataL1 %>%
      dplyr::filter(stationName == sStationName) %>%
      dplyr::arrange(date)
   
   dplyr::tbl_df(dfTimeDataL2)
   
   ggplot(data = dfTimeDataL2, aes(x = dtDate, ymin = dtSunrise, ymax = dtSunset)) +
      theme_bw() +
      geom_ribbon(fill = "#FDE725FF", alpha = 0.8) +
      scale_x_datetime(breaks = seq(min(dfTimeDataL2$dtDate), max(dfTimeDataL2$dtDate), "month"), expand = c(0, 0), labels = date_format("%b-%d\n%Y", tz = "Asia/Seoul"), minor_breaks = NULL) +
      scale_y_continuous(limits = c(4, 22), breaks = seq(4, 22, 4), expand = c(0, 0), minor_breaks = NULL) +
      labs(
         x = "Date [Month-Day Year]",
         y = "Sunrise and Sunset [Hour]",
         title = sprintf("Sunrise and Sunset for %s\n%s ", sStationName
                         , paste0(range(dfTimeData$date), sep = " ", collapse = "to "))
         ) +
      theme(
         panel.background = element_rect(fill = "#180F3EFF")
         , panel.grid = element_line(colour = "grey", linetype = "dashed")
         , plot.title = element_text(hjust = 0.5, 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")
         , text = element_text(family = font)
         , plot.margin = unit(c(5, 10, 5, 5), "mm")
         ) +
      ggsave(filename = sSaveFilePathName, width = 12, height = 8, dpi = 600)
}

 

• 강릉

  • 일출 시작/종료 : 07시 40분 55초 - 07시 44분 01초

  • 일몰 시간/종료 : 17시 13분 33초 - 17시 16분 39초 

 

그림. 2019년 1월 1일 (경자년)에 강릉의 일출 (해돋이) 및 일몰 시각 시계열.

 

• 인천

  • 일출 시작/종료 : 07시 49분 14초 - 07시 52분 18초

  • 일몰 시간/종료 : 17시 23분 24초 - 17시 26분 29초 

 

그림. 2019년 1월 1일 (경자년)에 인천의 일출 (해돋이) 및 일몰 시각 시계열.

 

• 부산

  • 일출 시작/종료 : 07시 33분 14초 - 07시 36분 11초

  • 일몰 시간/종료 : 17시 20분 15초 - 17시 23분 12초 

 

그림. 2019년 1월 1일 (경자년)에 부산의 일출 (해돋이) 및 일몰 시각 시계열.

 

[전체] 

 

 참고 문헌

[논문]

• 없음

[보고서]

• 없음

[URL]

• 없음

 

 문의사항

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

• sangho.lee.1990@gmail.com

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

• saimang0804@gmail.com