정보

    • 업무명    : 베이지안 선형회귀 모델을 통한 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 ##

    종관 기상 관측지점의 자료 예시
    AWS 기상 관측 지점의 자료 예시

     

    • 베이지안 모델 작성

      • 작성된 자료를 바탕으로 단변량 (종속변수가 하나인 경우) 인 경우에 대한 베이지안 모델을 작성 하였다.

      • 모델 작성에 사용된 함수의 하이퍼파라미터 등의 세부사항은 아래의 링크에서 확인이 가능하다.

     

    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
    • 네이버 블러그 공유하기
    • 네이버 밴드에 공유하기
    • 페이스북 공유하기
    • 카카오스토리 공유하기