[R] 시정거리 자료를 이용한 한반도 지역 등고선 (Contour) 시각화 및 해양 마스킹

정보

업무명 : 시정거리 자료를 이용한 한반도 지역 등고선 시각화 및 해양 마스킹
작성자 : 박진만
작성일 : 2021-01-17
설 명 :
수정이력 :

내용

[개요]

안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
미세먼지 데이터를 이용하여 한반도 영역의 Contour 그림 시각화 및 해양 마스킹 기능을 소개 합니다

[특징]

등고선 그림 그리기 수행

[기능]

시정 자료의 시각화

[활용 자료]

시정 관측 자료

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

자료처리 : tidyverse
내삽 패키지 : sp / grid
시각화 : ggplot2

[사용법]

소스 코드 참조

[사용 OS]

window 10

[사용 언어]

R 4.0.3

소스 코드

[명세]

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

library(rgdal)
library(ggplot2)
library(rgeos)
library(data.table)
library(dplyr)
library(scales)
library(metR)
library(grid)
#library(sf)
library(stringr)
library(RNetCDF)
library(ggpubr)

shp 파일 로드
- 전 세계 지도 데이터를 불러와 읽는다.
- 지도 형식은 shp 파일이다.

# shp 파일 로드 #
shape <- readOGR(dsn = path.expand("./input/geo/WorldCountries.shp"),
                 layer = "WorldCountries")

world = fortify(shape)

기상관측 데이터 및 한반도 영역 고도 데이터 읽어오기
- 기상관측 데이터의 경우 기상자료개방포털에서 구할 수 있다.
- 고도 데이터의 경우 미리 전처리 된 데이터를 읽는다. 이에 대해 자료 및 처리기법을 확인하고 싶은 경우 아래의 링크에서 확인 가능하다.

[재능상품] 특정 지점의 위치 좌표가 주어졌을 때 해당 좌표의 고도 구하기

정보 업무명 : 특정 지저의 위치 좌표가 주어졌을 때 해당 좌표의 고도 구하기 작성자 : 박진만 작성일 : 2020-12-05 설 명 : 수정이력 : 내용 [개요] 안녕하세요? 웹 개발 및 연구 개발을 담

shlee1990.tistory.com

#####
data = read.csv("./input/df_L2.csv")

station_data = fread("./input/station.txt",sep=" ")

hg_data = fread("./input/고도.csv")

unique(hg_data$lat)

hg_data = hg_data %>%
  dplyr::mutate(lat = round(lat,2),
                lon = round(lon,2)) %>%
  dplyr::distinct() %>%
  dplyr::arrange(lat,lon)

#####
data_L1 = inner_join(data,station_data,by="station")

head(data_L1)

#####

선택된 날짜에 대한 내삽/외삽 수행
- 기상 관측 지역은 한정되어 있으므로, 이를 격자화 시켜주기 위함.
- 결과적으로 내삽 및 외삽을 수행하여 그리드화 시켜준다.

dfData_res = data.frame()
count = 0

for (i in c("evis","vis","rh","pm25")) {
  
  count = count + 1
  select_col = c("lon","lat",i)
  
  dfData <- data_L1 %>%
    dplyr::filter(date == "201901130600") %>%
    dplyr::select(select_col)
  
  colnames(dfData)[1:3] = c("lon","lat","val")
  
  nNewLon = seq(from = xRange[1], to = xRange[2], by = 0.01)
  nNewLat = seq(from = yRange[1], to = yRange[2], by = 0.01)
  
  dfTmp = akima::interp(
    dfData$lon
    , dfData$lat
    , dfData$val
    , xo = nNewLon
    , yo = nNewLat
    , linear = FALSE
    , extrap = TRUE
  )
  
  
  # Expand Points to Grid
  spGeoData = expand.grid(nNewLon, nNewLat)
  
  as.data.frame(expand.grid(nNewLon, nNewLat))
  
  coordinates(spGeoData) = ~ Var1 + Var2
  gridded(spGeoData) = TRUE
  
  # L2 Processing Using L1 Data Frame
  dfDataL2 = data.frame(
    nLon = spGeoData$Var1
    , nLat = spGeoData$Var2
    , nVal = c(dfTmp$z)
  )
  
  # dplyr::glimpse(dfDataL2)
  colnames(dfDataL2)<- c("lon","lat","val")

  if(count == 1){
    dfData_res = dfDataL2
    colnames(dfData_res)[dim(dfData_res)[2]] = i
  } else {
    dfData_res = cbind(dfData_res,dfDataL2$val)
    colnames(dfData_res)[dim(dfData_res)[2]] = i
  }
  
}

내삽/외삽 결과 (외삽을 같이 수행하는 경우 값이 정확하지 않을 수 있다.)

그리드화 된 데이터의 후처리
- 시각화를 위해 해양 지역 (고도가 0인 지역) 의 경우 마스킹을 수행한다.
- 자료 표출 범위를 나타내기 위해 물리적으로 틀린 값이 결과로 나오는 경우 이를 제거한다.

dfDataL3 = dfData_res %>%
  dplyr::mutate(visdiff = evis - vis) %>%
  dplyr::mutate(evis = ifelse(evis > 20,20.0,evis)) %>%
  dplyr::mutate(evis = ifelse(evis < 1, 1.0, evis)) %>%
  dplyr::mutate(vis = ifelse(vis > 20,20.0,evis)) %>%
  dplyr::mutate(vis = ifelse(vis < 1, 1.0, evis)) %>%
  dplyr::mutate(visdiff = ifelse(visdiff > 5,5,visdiff)) %>%
  dplyr::mutate(visdiff = ifelse(visdiff < -5, -5, visdiff)) %>%
  dplyr::mutate(rh = ifelse(rh > 80,80,rh)) %>%
  dplyr::mutate(rh = ifelse(rh < 20,20, rh)) %>%
  dplyr::mutate(pm25 = ifelse(pm25 > 100,100,pm25)) %>%
  dplyr::mutate(pm25 = ifelse(pm25 < 0, 0, pm25)) %>%
  dplyr::arrange(lat,lon) %>%
  dplyr::select(-lat,-lon)

summary(dfDataL3)

dfDataL4 = cbind(dfDataL3,hg_data)

dfDataL5 = dfDataL4 %>%
  dplyr::mutate(evis = ifelse(high == 0, NA, evis)) %>%
  dplyr::mutate(vis = ifelse(high == 0, NA, vis)) %>%
  dplyr::mutate(rh = ifelse(high == 0, NA, rh)) %>%
  dplyr::mutate(pm25 = ifelse(high == 0, NA, pm25)) %>%
  dplyr::mutate(visdiff = ifelse(high == 0, NA, visdiff))

시정 시각화
- 처리된 결과를 이용하여 매핑을 수행한다.
- ggplot2 패키지를 이용함.

ggplot() +
  coord_fixed() +
  theme_bw() +
  geom_polygon(data=world, aes(x=long,y=lat,group=group), fill="white") + # 지도 외곽선 내부색깔
  geom_raster(data=dfDataL5,aes(lon, lat, fill = vis),interpolate = TRUE) +
  scale_fill_gradientn(colours = rev(c("blueviolet","blue","green","yellow","orange","red","darkred")), # 컬러바 옵션 조정
                       breaks = seq(0,20,2), limits = c(0,20),na.value = "white") +
  geom_contour(data = dfDataL5, aes(lon, lat, z = vis),colour='black',breaks = seq(0,20,2)) +
  geom_text_contour(breaks = seq(0,20,2),data = dfDataL5, aes(lon, lat, z = vis),colour='black') +
  geom_path(data=world, aes(x=long,y=lat,group=group), colour='black') + # 지도 외곽선 색깔
  geom_point(data = dfStation,aes(x=lon,y=lat),color = "black", size = 2.0) + 
  theme(plot.title=element_text(face="bold", size=12, 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=20, colour="black")) +
  theme(legend.key.heigh= unit(4,'cm')) +
  theme(legend.key.width= unit(2,'cm')) +
  theme(legend.key=element_blank()) +
  theme(legend.text=element_text(size=15, face="bold")) +
  theme(legend.background=element_blank()) +
  scale_x_continuous(limits = c(minlon,maxlon),expand = c(0,0)) +
  scale_y_continuous(limits = c(minlat,40),expand = c(0,0)) +
  ggsave("./vis.png",dpi = 150,height = 14,width = 14)

[전체 코드]

library(rgdal)
library(ggplot2)
library(rgeos)
library(data.table)
library(dplyr)
library(scales)
library(metR)
library(grid)
#library(sf)
library(stringr)
library(RNetCDF)
library(ggpubr)


###########################################################

# shp 파일 로드 #
shape <- readOGR(dsn = path.expand("./input/geo/WorldCountries.shp"),
                 layer = "WorldCountries")

world = fortify(shape)

###########################################################

#####
data = read.csv("./input/df_L2.csv")

station_data = fread("./input/station.txt",sep=" ")

hg_data = fread("./input/고도.csv")

unique(hg_data$lat)

hg_data = hg_data %>%
  dplyr::mutate(lat = round(lat,2),
                lon = round(lon,2)) %>%
  dplyr::distinct() %>%
  dplyr::arrange(lat,lon)

#####
data_L1 = inner_join(data,station_data,by="station")

head(data_L1)

#####

dfStation = data_L1 %>%
  dplyr::filter(date == "201901011300") 

################################################################################
colnames(data_L1)

dfData_res = data.frame()
count = 0

for (i in c("evis","vis","rh","pm25")) {
  
  count = count + 1
  select_col = c("lon","lat",i)
  
  dfData <- data_L1 %>%
    dplyr::filter(date == "201901130600") %>%
    dplyr::select(select_col)
  
  colnames(dfData)[1:3] = c("lon","lat","val")
  
  nNewLon = seq(from = xRange[1], to = xRange[2], by = 0.01)
  nNewLat = seq(from = yRange[1], to = yRange[2], by = 0.01)
  
  dfTmp = akima::interp(
    dfData$lon
    , dfData$lat
    , dfData$val
    , xo = nNewLon
    , yo = nNewLat
    , linear = FALSE
    , extrap = TRUE
  )
  
  
  # Expand Points to Grid
  spGeoData = expand.grid(nNewLon, nNewLat)
  
  as.data.frame(expand.grid(nNewLon, nNewLat))
  
  coordinates(spGeoData) = ~ Var1 + Var2
  gridded(spGeoData) = TRUE
  
  # L2 Processing Using L1 Data Frame
  dfDataL2 = data.frame(
    nLon = spGeoData$Var1
    , nLat = spGeoData$Var2
    , nVal = c(dfTmp$z)
  )
  
  # dplyr::glimpse(dfDataL2)
  colnames(dfDataL2)<- c("lon","lat","val")

  if(count == 1){
    dfData_res = dfDataL2
    colnames(dfData_res)[dim(dfData_res)[2]] = i
  } else {
    dfData_res = cbind(dfData_res,dfDataL2$val)
    colnames(dfData_res)[dim(dfData_res)[2]] = i
  }
  
}


dfDataL3 = dfData_res %>%
  dplyr::mutate(visdiff = evis - vis) %>%
  dplyr::mutate(evis = ifelse(evis > 20,20.0,evis)) %>%
  dplyr::mutate(evis = ifelse(evis < 1, 1.0, evis)) %>%
  dplyr::mutate(vis = ifelse(vis > 20,20.0,evis)) %>%
  dplyr::mutate(vis = ifelse(vis < 1, 1.0, evis)) %>%
  dplyr::mutate(visdiff = ifelse(visdiff > 5,5,visdiff)) %>%
  dplyr::mutate(visdiff = ifelse(visdiff < -5, -5, visdiff)) %>%
  dplyr::mutate(rh = ifelse(rh > 80,80,rh)) %>%
  dplyr::mutate(rh = ifelse(rh < 20,20, rh)) %>%
  dplyr::mutate(pm25 = ifelse(pm25 > 100,100,pm25)) %>%
  dplyr::mutate(pm25 = ifelse(pm25 < 0, 0, pm25)) %>%
  dplyr::arrange(lat,lon) %>%
  dplyr::select(-lat,-lon)

summary(dfDataL3)

dfDataL4 = cbind(dfDataL3,hg_data)

dfDataL5 = dfDataL4 %>%
  dplyr::mutate(evis = ifelse(high == 0, NA, evis)) %>%
  dplyr::mutate(vis = ifelse(high == 0, NA, vis)) %>%
  dplyr::mutate(rh = ifelse(high == 0, NA, rh)) %>%
  dplyr::mutate(pm25 = ifelse(high == 0, NA, pm25)) %>%
  dplyr::mutate(visdiff = ifelse(high == 0, NA, visdiff)) 
  
  
  
  ggplot() +
  coord_fixed() +
  theme_bw() +
  geom_polygon(data=world, aes(x=long,y=lat,group=group), fill="white") + # 지도 외곽선 내부색깔
  geom_raster(data=dfDataL5,aes(lon, lat, fill = vis),interpolate = TRUE) +
  scale_fill_gradientn(colours = rev(c("blueviolet","blue","green","yellow","orange","red","darkred")), # 컬러바 옵션 조정
                       breaks = seq(0,20,2), limits = c(0,20),na.value = "white") +
  geom_contour(data = dfDataL5, aes(lon, lat, z = vis),colour='black',breaks = seq(0,20,2)) +
  geom_text_contour(breaks = seq(0,20,2),data = dfDataL5, aes(lon, lat, z = vis),colour='black') +
  geom_path(data=world, aes(x=long,y=lat,group=group), colour='black') + # 지도 외곽선 색깔
  geom_point(data = dfStation,aes(x=lon,y=lat),color = "black", size = 2.0) + 
  theme(plot.title=element_text(face="bold", size=12, 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=20, colour="black")) +
  theme(legend.key.heigh= unit(4,'cm')) +
  theme(legend.key.width= unit(2,'cm')) +
  theme(legend.key=element_blank()) +
  theme(legend.text=element_text(size=15, face="bold")) +
  theme(legend.background=element_blank()) +
  scale_x_continuous(limits = c(minlon,maxlon),expand = c(0,0)) +
  scale_y_continuous(limits = c(minlat,40),expand = c(0,0)) +
  ggsave("./vis.png",dpi = 150,height = 14,width = 14)

참고 문헌

[논문]

없음

[보고서]

없음

[URL]

없음

문의사항

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

sangho.lee.1990@gmail.com

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

saimang0804@gmail.com

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

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

[R] 끄투 자동 게임 수행 프로그램 (test) (11)	2021.01.22
[R] WRF 모델 자료로부터 기상 변수 추출 및 단열선도 시각화 그리고 RadioSonde 패키지 내 단열선도 X축 커스터마이징 (켈빈 to 섭씨) (0)	2021.01.18
[R] 베이지안 선형회귀 모델을 통한 ASOS 종관 기상 관측지점 자료를 이용하여 AWS 관측 지점의 기온 추정 (0)	2021.01.17
[R] 구글 번역 API 및 R 프로그램을 이용하여 마인크래프트 1.12.1 모드 파일 자동으로 번역 후 적용하기 (구글 API 신청 방법 포함) (0)	2020.12.28
[R] '끄투코리아' 사이트를 매개체로 한 강화학습을 통한 끝말잇기 봇 프로그램의 학습 알고리즘 구현 및 테스트 (0)	2020.11.28

정보

내용