[R] 베이지안 선형회귀 모델을 통한 ASOS 종관 기상 관측지점 자료를 이용하여 AWS 관측 지점의 기온 추정

정보

업무명 : 베이지안 선형회귀 모델을 통한 ASOS 종관 기상 관측지점 자료를 이용하여 AWS 관측 지점의 기온 추정
작성자 : 박진만
작성일 : 2021-01-16
설 명 :
수정이력 :

내용

[개요]

안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
공간 자료를 활용한 기상 자료의 추정의 경우 공간적 효과를 고려해야 더욱 정확한 추정이 가능한 것으로 알려져 있습니다.
이는 기상학적 데이터가 지형지물과 시간 등의 영향을 절대적으로 받고있기 때문입니다.
따라서 이러한 공간적 추정을 일반적 선형모형으로 가정하게 되는 경우 발생하는 문제가 모형 추정오류의 원인이 될 수도 있습니다.
따라서 본 글에서는 일반적인 선형 모형 또는 내삽에 의한 기상 자료의 공간적 추정이 아닌 자기회귀 베이지안 모델을 이용하여 기상 변수를 추정하는 방법을 제시 해 보고자 합니다.
실제 코드는 R을 통해 구현되었으며, 종관 기상 관측 지점의 기온을 이용하여, 위도/경도/고도를 독립변수로 하여, 실제 AWS 기상 관측 지점의 기온을 추정하고 이를 비교 해 보고자 합니다.

[특징]

R을 통해 베이지안 모델을 설계하고 공간적 추정을 실시 하였.

[기능]

단일변량 베이지안 추정을 통한 기온 추정

[활용 자료]

자료명 : 종관 기상 관측자료 / AWS 기상 관측자료
자료 종류 : 시간별 기상관측 자료 및 지점별 위치자료
확장자 : csv
영역 : 대한민국 전역
기간 : 2019년 1년간의 자료를 사용하였음

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

자료처리 : tidyverse

주요 분석 패키지 : Spatial 및 splm
시각화 : ggplot2

[사용법]

소스 코드 참조

[사용 OS]

window 10

[사용 언어]

R 4.0.3

소스 코드

[명세]

라이브러리 불러오기
- 실행에 필요한 패키지를 불러온다.

library(spatstat)
library(spatial)
library(splancs)
library(gstat)
library(geoR)
library(fields)
library(spBayes)
library(RandomFields)
library(sgeostat)
library(vardiag)
library(spdep)
library(DCluster)
library(spgwr)
library(ade4)
library(mgcv)
library(R2WinBUGS)
library(R2BayesX)
library(dplyr)
library(tidyr)
library(tidyverse)
library(Metrics)
library(ggplot2)
library(Rfast)
library(dplyr)
library(slam)
library(tidyr)
library(readr)
library(data.table)

기상 관측자료 읽기
- 종관 기상관측자료 및 AWS 기상관측 자료를 불러온다.
- 전자는 베이지안 회귀 모형의 학습에 사용되며, 후자는 만들어진 모델의 검증에 사용되었다.

## LOAD AWS AND ASOS DATA ##
ASOS_DATA <- read.csv("./INPUT/OBS_ASOS_TIM_20200227022304.csv",stringsAsFactors = F)

colnames(ASOS_DATA) <- c("num","name","date","T","RAIN","WS","WD","RH")

AWS_DATA1 <- read.csv("./INPUT/OBS_AWS_TIM_20200227023747.csv",stringsAsFactors = F)
AWS_DATA2 <- read.csv("./INPUT/OBS_AWS_TIM_20200227024521.csv",stringsAsFactors = F)

AWS_DATA <- rbind(AWS_DATA1,AWS_DATA2)

colnames(AWS_DATA) <- c("num","name","date","T","WD","WS","RAIN","RH")

head(AWS_DATA)
## LOAD AWS AND ASOS DATA ##

관측지점 위치 데이터 불러오기
- 각 지점마다의 위치 자료 (위도/경도/고도) 를 불러온다.

## LOAD LOCATE DATA ##
ASOS_LOCATE <- read.csv("./LOCATE/ASOS/ASOS_LOCATE.csv",stringsAsFactors = F)

colnames(ASOS_LOCATE) <- c("num","sdate","edate","name","address","address2","lat","lon","alt",
                           "dum1","dum2","dum3","dum4")

ASOS_LOCATE_L1 <- ASOS_LOCATE %>%
  dplyr::select(num:name,lat:alt) %>%
  dplyr::mutate(edate = ifelse(edate == "","9999-12-31",edate))


AWS_LOCATE <- read.csv("./LOCATE/AWS/AWS_LOCATE.csv",stringsAsFactors = F)

colnames(AWS_LOCATE) <- c("num","sdate","edate","name","lat","lon","alt",
                          "dum1","dum2","dum3","dum4")

AWS_LOCATE_L1 <- AWS_LOCATE %>%
  dplyr::select(num:name,lat:alt) %>%
  dplyr::mutate(edate = ifelse(edate == "","9999-12-31",edate))

## LOAD LOCATE DATA ##

구하고자 하는 기간에 대한 자료 전처리 작업
- 특정 시점에 대한 베이지안 공간 회귀 모형을 얻고 싶은 경우 기상자료의 날짜를 필터링 해 주어야 한다.
- 예를 들어 2019년 01월 01일 01 시에 대한 베이지안 추정을 하고 싶은 경우 해당 시간에 대한 기상 자료만을 추출한 후 모형에 사용 하여야 한다.
- 따라서 분석을 시작하고자 하는 시간 - 끝 시간을 미리 설정한 후 이를 사전에 전처리 하여 원하는 시간대의 기상 자료만을 추출하고, 동시에 해당 시간에 관측을 수행하고 있던 지역을 찾아 이를 매칭시켜주는 작업이 진행되어야 한다.

## FILTER DATE ##
TS <- as.POSIXct("2019-01-01 01:00")
TE <- as.POSIXct("2019-01-01 01:00")

TARGET_DATE_DATE <- as.Date(substr(TE,1,10))
TARGET_DATE_START <- seq(TS,TE,3600)
TARGET_DATE_END <- seq(TS,TE,3600)
TARGET_DATE_DATE <- as.Date(substr(TARGET_DATE_END,1,10))

ASOS_DATA_TARGET <- ASOS_DATA %>%
  dplyr::filter(date >= TARGET_DATE_START) %>%
  dplyr::filter(date <= TARGET_DATE_END) %>%
  dplyr::filter(!is.na(T)) %>%
  dplyr::group_by(num,name) %>%
  dplyr::summarise(T = mean(T),
                   N = n())

AWS_DATA_TARGET <- AWS_DATA %>%
  dplyr::filter(date >= TARGET_DATE_START) %>%
  dplyr::filter(date <= TARGET_DATE_END) %>%
  dplyr::filter(!is.na(T)) %>%
  dplyr::group_by(num,name) %>%
  dplyr::summarise(T = mean(T),
                   N = n())

ASOS_LOCATE_TARGET <- ASOS_LOCATE_L1 %>%
  dplyr::filter(as.Date(sdate) <= TARGET_DATE_DATE) %>%
  dplyr::filter(as.Date(edate) >= TARGET_DATE_DATE) 

AWS_LOCATE_TARGET <- AWS_LOCATE_L1 %>%
  dplyr::filter(as.Date(sdate) <= TARGET_DATE_DATE) %>%
  dplyr::filter(as.Date(edate) >= TARGET_DATE_DATE)


ASOS_FIN_DATA <- dplyr::inner_join(ASOS_DATA_TARGET,ASOS_LOCATE_TARGET,by=c("num","name"))
AWS_FIN_DATA <- dplyr::inner_join(AWS_DATA_TARGET,AWS_LOCATE_TARGET,by=c("num","name"))
AWS_FIN_DATA <- AWS_FIN_DATA[!duplicated(AWS_FIN_DATA$lat), ]
## FILTER DATE ##

베이지안 모델 작성
- 작성된 자료를 바탕으로 단변량 (종속변수가 하나인 경우) 인 경우에 대한 베이지안 모델을 작성 하였다.
- 모델 작성에 사용된 함수의 하이퍼파라미터 등의 세부사항은 아래의 링크에서 확인이 가능하다.

spLM function | R Documentation

Looks like there are no examples yet.

www.rdocumentation.org

form <- T~lat+lon+alt-1
coords <- as.matrix(cbind(ASOS_FIN_DATA$lat,ASOS_FIN_DATA$lon))
coords2 <- as.matrix(cbind(AWS_FIN_DATA$lat,AWS_FIN_DATA$lon))

X <- as.matrix(cbind(ASOS_FIN_DATA$alt,0))
X2 <- as.matrix(cbind(AWS_FIN_DATA$lat,AWS_FIN_DATA$lon,AWS_FIN_DATA$alt))

starting <- list("phi"=3/0.5, "sigma.sq"=50, "tau.sq"=1)

tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)

priors.2 <- list("beta.Flat", "phi.Unif"=c(3/1, 3/0.1),
                 "sigma.sq.IG"=c(2, 2), "tau.sq.IG"=c(2, 0.1))

cov.model <- "exponential"

n.samples <- 20000

m.1 <- spLM(formula = form, coords=coords, starting=starting, tuning=tuning, data = ASOS_FIN_DATA,
            priors=priors.2, cov.model=cov.model, n.samples=20000,n.report=1000)

베이지안 모델 적용 및 검증
- 작성된 모델을 이용하여, 실제 AWS 관측 지점의 추정 기온을 확인하고 이를 실제 AWS 관측 지점의 기온과 비교 해 본다.

m.1.pred <- spPredict(m.1, pred.covars=as.matrix(X2), pred.coords=coords2,
                      start=0.5*n.samples)

y.hat <- apply(m.1.pred$p.y.predictive.samples, 1, mean)

aaa <- as.data.frame((y.hat))

result <- cbind(AWS_FIN_DATA$T,aaa)
colnames(result) <- c("AWS_T","cal_T")

X1 <- result$AWS_T
Y1 <- result$cal_T

ggplot() +
  coord_fixed(ratio=1) +
  theme_bw() +
  geom_point(aes(X1, Y1)) +
  geom_abline(intercept=0, slope=1, linetype=1, color="black", size=1.0) +
  labs(title = paste0("Bayesian Spatial Regression Models AWS Temperature Result") ) +
  labs(x = expression(paste(bold("AWS Temperature (obs.℃)"))),
       y = expression(paste(bold("AWS Temperature (cal.℃)"))),
       fill = "Count") +
  theme(plot.title=element_text(face="bold", size=20, color="black")) +
  theme(axis.title.x = element_text(face="bold", size=19, colour="black")) +
  theme(axis.title.y = element_text(face="bold", size=19, colour="black", angle=90)) +
  theme(axis.text.x = element_text(face="bold", size=19, colour="black")) +
  theme(axis.text.y = element_text(face="bold", size=19, colour="black")) +
  theme(legend.title=element_text(face="bold", size=14, colour="black")) +
  theme(legend.position=c(0,1), legend.justification=c(0,0.96)) +
  theme(legend.key=element_blank()) +
  theme(legend.text=element_text(size=14, face="bold")) +
  theme(legend.background=element_blank()) +
  #theme(text=element_text(family=font)) +
  theme(plot.margin=unit(c(0,8,0,0),"mm"))

[전체 코드]

library(spatstat)
library(spatial)
library(splancs)
library(gstat)
library(geoR)
library(fields)
library(spBayes)
library(RandomFields)
library(sgeostat)
library(vardiag)
library(spdep)
library(DCluster)
library(spgwr)
library(ade4)
library(mgcv)
library(R2WinBUGS)
library(R2BayesX)
library(dplyr)
library(tidyr)
library(tidyverse)
library(Metrics)
library(ggplot2)
library(Rfast)
library(dplyr)
library(slam)
library(tidyr)
library(readr)
library(data.table)


## LOAD AWS AND ASOS DATA ##
ASOS_DATA <- read.csv("./INPUT/OBS_ASOS_TIM_20200227022304.csv",stringsAsFactors = F)

colnames(ASOS_DATA) <- c("num","name","date","T","RAIN","WS","WD","RH")

AWS_DATA1 <- read.csv("./INPUT/OBS_AWS_TIM_20200227023747.csv",stringsAsFactors = F)
AWS_DATA2 <- read.csv("./INPUT/OBS_AWS_TIM_20200227024521.csv",stringsAsFactors = F)

AWS_DATA <- rbind(AWS_DATA1,AWS_DATA2)

colnames(AWS_DATA) <- c("num","name","date","T","WD","WS","RAIN","RH")

head(AWS_DATA)
## LOAD AWS AND ASOS DATA ##

## LOAD LOCATE DATA ##
ASOS_LOCATE <- read.csv("./LOCATE/ASOS/ASOS_LOCATE.csv",stringsAsFactors = F)

colnames(ASOS_LOCATE) <- c("num","sdate","edate","name","address","address2","lat","lon","alt",
                           "dum1","dum2","dum3","dum4")

ASOS_LOCATE_L1 <- ASOS_LOCATE %>%
  dplyr::select(num:name,lat:alt) %>%
  dplyr::mutate(edate = ifelse(edate == "","9999-12-31",edate))


AWS_LOCATE <- read.csv("./LOCATE/AWS/AWS_LOCATE.csv",stringsAsFactors = F)

colnames(AWS_LOCATE) <- c("num","sdate","edate","name","lat","lon","alt",
                          "dum1","dum2","dum3","dum4")

AWS_LOCATE_L1 <- AWS_LOCATE %>%
  dplyr::select(num:name,lat:alt) %>%
  dplyr::mutate(edate = ifelse(edate == "","9999-12-31",edate))

## LOAD LOCATE DATA ##
##########################################################################################


head(AWS_DATA)

TS <- as.POSIXct("2019-01-01 01:00")
TE <- as.POSIXct("2019-01-01 01:00")

TARGET_DATE_DATE <- as.Date(substr(TE,1,10))
TARGET_DATE_START <- seq(TS,TE,3600)
TARGET_DATE_END <- seq(TS,TE,3600)
TARGET_DATE_DATE <- as.Date(substr(TARGET_DATE_END,1,10))

ASOS_DATA_TARGET <- ASOS_DATA %>%
  dplyr::filter(date >= TARGET_DATE_START) %>%
  dplyr::filter(date <= TARGET_DATE_END) %>%
  dplyr::filter(!is.na(T)) %>%
  dplyr::group_by(num,name) %>%
  dplyr::summarise(T = mean(T),
                   N = n())

AWS_DATA_TARGET <- AWS_DATA %>%
  dplyr::filter(date >= TARGET_DATE_START) %>%
  dplyr::filter(date <= TARGET_DATE_END) %>%
  dplyr::filter(!is.na(T)) %>%
  dplyr::group_by(num,name) %>%
  dplyr::summarise(T = mean(T),
                   N = n())

ASOS_LOCATE_TARGET <- ASOS_LOCATE_L1 %>%
  dplyr::filter(as.Date(sdate) <= TARGET_DATE_DATE) %>%
  dplyr::filter(as.Date(edate) >= TARGET_DATE_DATE) 

AWS_LOCATE_TARGET <- AWS_LOCATE_L1 %>%
  dplyr::filter(as.Date(sdate) <= TARGET_DATE_DATE) %>%
  dplyr::filter(as.Date(edate) >= TARGET_DATE_DATE)


ASOS_FIN_DATA <- dplyr::inner_join(ASOS_DATA_TARGET,ASOS_LOCATE_TARGET,by=c("num","name"))
AWS_FIN_DATA <- dplyr::inner_join(AWS_DATA_TARGET,AWS_LOCATE_TARGET,by=c("num","name"))
AWS_FIN_DATA <- AWS_FIN_DATA[!duplicated(AWS_FIN_DATA$lat), ]
## FILTER DATE ##

form <- T~lat+lon+alt-1
coords <- as.matrix(cbind(ASOS_FIN_DATA$lat,ASOS_FIN_DATA$lon))
coords2 <- as.matrix(cbind(AWS_FIN_DATA$lat,AWS_FIN_DATA$lon))

X <- as.matrix(cbind(ASOS_FIN_DATA$alt,0))
X2 <- as.matrix(cbind(AWS_FIN_DATA$lat,AWS_FIN_DATA$lon,AWS_FIN_DATA$alt))

starting <- list("phi"=3/0.5, "sigma.sq"=50, "tau.sq"=1)

tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1)

priors.2 <- list("beta.Flat", "phi.Unif"=c(3/1, 3/0.1),
                 "sigma.sq.IG"=c(2, 2), "tau.sq.IG"=c(2, 0.1))

cov.model <- "exponential"

n.samples <- 20000

m.1 <- spLM(formula = form, coords=coords, starting=starting, tuning=tuning, data = ASOS_FIN_DATA,
            priors=priors.2, cov.model=cov.model, n.samples=20000,n.report=1000)

m.1.pred <- spPredict(m.1, pred.covars=as.matrix(X2), pred.coords=coords2,
                      start=0.5*n.samples)

y.hat <- apply(m.1.pred$p.y.predictive.samples, 1, mean)

aaa <- as.data.frame((y.hat))

result <- cbind(AWS_FIN_DATA$T,aaa)
colnames(result) <- c("AWS_T","cal_T")

X1 <- result$AWS_T
Y1 <- result$cal_T

ggplot() +
  coord_fixed(ratio=1) +
  theme_bw() +
  geom_point(aes(X1, Y1)) +
  geom_abline(intercept=0, slope=1, linetype=1, color="black", size=1.0) +
  labs(title = paste0("Bayesian Spatial Regression Models AWS Temperature Result") ) +
  labs(x = expression(paste(bold("AWS Temperature (obs.℃)"))),
       y = expression(paste(bold("AWS Temperature (cal.℃)"))),
       fill = "Count") +
  theme(plot.title=element_text(face="bold", size=20, color="black")) +
  theme(axis.title.x = element_text(face="bold", size=19, colour="black")) +
  theme(axis.title.y = element_text(face="bold", size=19, colour="black", angle=90)) +
  theme(axis.text.x = element_text(face="bold", size=19, colour="black")) +
  theme(axis.text.y = element_text(face="bold", size=19, colour="black")) +
  theme(legend.title=element_text(face="bold", size=14, colour="black")) +
  theme(legend.position=c(0,1), legend.justification=c(0,0.96)) +
  theme(legend.key=element_blank()) +
  theme(legend.text=element_text(size=14, face="bold")) +
  theme(legend.background=element_blank()) +
  #theme(text=element_text(family=font)) +
  theme(plot.margin=unit(c(0,8,0,0),"mm"))

참고 문헌

[논문]

없음

[보고서]

없음

[URL]

없음

문의사항

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

sangho.lee.1990@gmail.com

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

saimang0804@gmail.com

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

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

[R] WRF 모델 자료로부터 기상 변수 추출 및 단열선도 시각화 그리고 RadioSonde 패키지 내 단열선도 X축 커스터마이징 (켈빈 to 섭씨) (0)	2021.01.18
[R] 시정거리 자료를 이용한 한반도 지역 등고선 (Contour) 시각화 및 해양 마스킹 (2)	2021.01.17
[R] 구글 번역 API 및 R 프로그램을 이용하여 마인크래프트 1.12.1 모드 파일 자동으로 번역 후 적용하기 (구글 API 신청 방법 포함) (0)	2020.12.28
[R] '끄투코리아' 사이트를 매개체로 한 강화학습을 통한 끝말잇기 봇 프로그램의 학습 알고리즘 구현 및 테스트 (0)	2020.11.28
[R] 특정 날짜의 네이버 뉴스 제목 크롤링 (0)	2020.11.26

[R] 베이지안 선형회귀 모델을 통한 ASOS 종관 기상 관측지점 자료를 이용하여 AWS 관측 지점의 기온 추정

정보

내용

[개요]

[특징]

[기능]

[활용 자료]

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

[사용법]

[사용 OS]

[사용 언어]

소스 코드

[명세]

[전체 코드]

참고 문헌

[논문]

[보고서]

[URL]

문의사항

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

전체 카테고리

블로그 인기글

전체 방문자

티스토리툴바

정보

내용

[개요]

[특징]

[기능]

[활용 자료]

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

[사용법]

[사용 OS]

[사용 언어]

소스 코드

[명세]

[전체 코드]

참고 문헌

[논문]

[보고서]

[URL]

문의사항

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

전체 카테고리

블로그 인기글

최근 글

최근댓글

전체 방문자

티스토리툴바