정보
-
업무명 : KLAPS 수치예측 모델 자료를 이용하여 내삽 방법 (Multi level B-Spline Approximation, Kriging)에 따른 전처리 및 가시화
-
작성자 : 이상호
-
작성일 : 2020-04-13
-
설 명 :
-
수정이력 :
내용
[개요]
-
안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
-
이전 포스팅에서 역거리 가중치 (Inverse Distance Weighting) 및 선형 내삽 (Linear Interpolation) 방법을 설명했습니다. 이 방법은 간단하나 단점이 있습니다. 즉 내삽의 정확성이 다소 부정확하는 것입니다.
[R] KLAPS 수치예측 모델 자료를 이용하여 내삽 방법 (Inverse Distance Weighting, Linear Interpolation)에 따른 전처리 및 가시화
정보 업무명 : KLAPS 수치예측 모델 자료를 이용하여 내삽 방법 (Inverse Distance Weighting, Linear Interpolation)에 따른 전처리 및 가시화 작성자 : 이상호 작성일 : 2020-01-03 설 명 : 수정이력 : 내용 [특..
shlee1990.tistory.com
-
오늘 포스팅은 최근에 사용되는 내삽 방법 (Multi level B-Spline Approximation, Kriging)을 통해 전처리 및 가시화를 소개해 드리고자 합니다.
[특징]
-
KLAPS 수치예측 모델 자료 및 내삽 방법을 이해하기 위해서 가시화가 요구되며 이 프로그램은 이러한 목적을 달성하기 위한 소프트웨어
[기능]
-
KLAPS 수치예측 모델 자료를 이용한 전처리 및 가시화
-
다단계 B-Spline 근사화 (Multi level B-Spline Approximation, MBA)를 통해 전처리 및 가시화
-
정규 크리킹 (Kriging)을 통해 전처리 및 가시화
[활용 자료]
-
자료명 : 초단기예보모델 (KLAPS)
-
자료 종류 : 1100 hPa에서 지표면 온도
-
영역 : 한반도
-
확장자 : NetCDF
-
기간 : 2019년 07월 03일 09시 00분 (KST)
-
자료 설명 : 한반도 영역을 대상으로 3차원 예측자료를 생산하며 주로 한반도 영역의 3차원 분석/예측에 이용
[자료 처리 방안 및 활용 분석 기법]
-
다단계 B-Spline 근사화 (Multi level B-Spline Approximation)를 통해 3D 내삽
-
R에서 규칙 및 불규칙 격자에 대해 데이터 내삽 방법로서 MBA 패키지 사용
-
최근 많은 이미지 영상과 관련하여 다양한 연구 분야에서 데이터 보간에 대한 요구는 증가 추세에 있으며 이에 대해 많은 방법이 제안되어있다.
-
이 중 B-Spline 근사화는 약간의 오차를 함유하지만 자연스럽고 매끄러운 형상 변화가 가능한 방법이다.
-
또한 이러한 오차 감소 기법으로서 다단계 B-Spline 근사화 방법이 제안되었다. 이 방법은 B-spline 근사화의 장점 (자연그러움, 매끄러움)을 유지하면서 데이터 보간을 해준다.
-
그러나 정밀도 향상에 따라 필요한 계산량이 기존의 B-Spline 근사화에 비해 크게 증가한다. 이는 사용상의 가장 큰 병목이 되나 최근에 많은 계산 및 시간 절감이 되었다.
-
[사용법]
-
입력 자료를 동일 디렉터리 위치
-
소스 코드를 실행
(Rscript Visualization_and_Preprocessing_According_to_Kriging_and_MBA_Interpolation_
Method_using_KLAPS_numerical_Prediction_Model_Data.R) -
가시화 결과를 확인
[사용 OS]
-
Windows 10
[사용 언어]
-
R v3.6.3
-
R Studio v1.2.5033
소스 코드
[명세]
-
전역 설정
-
최대 10 자리 설정
-
메모리 해제
-
영어 인코딩 설정
-
폰트 설정
-
# Set Option
options(digits = 10)
memory.limit(size = 9999999999999)
Sys.setlocale("LC_TIME", "english")
font = "Palatino Linotype"
-
라이브러리 읽기
# Library Load
library(RNetCDF)
library(tidyverse)
library(lubridate)
library(data.table)
library(colorRamps)
library(lubridate)
library(extrafont)
library(ggrepel)
library(scales)
library(sf)
library(gstat)
library(sp)
library(metR)
library(akima)
library(stringr)
library(noncompliance)
library(gstat)
library(automap)
-
위/경도 파일 읽기
# File Geo Read
dfGeo = data.table::fread("INPUT/klaps_reanal_lation.txt", sep = "\t", header = FALSE)
colnames(dfGeo) = c("lon", "lat")
dplyr::tbl_df(dfGeo)
-
NetCDF 파일 검색
# File NetCDF Read
sInputFileDirName = "INPUT/*.nc"
sFileDirName = Sys.glob(sInputFileDirName)
length(sFileDirName)
sFileDirName
-
NetCDF 파일 열기
# NetCDF File Open
ncFile = open.nc(sFileDirName[iCount])
-
출력/그림 파일명 설정
# Split File Name
sFileDirNameSplit = unlist(str_split(string = sFileDirName[iCount], pattern = "_|/|\\."))
# Get Date Time
dtDateTimeUtc = lubridate::ymd_hm(sFileDirNameSplit[5])
# Convert UTC to KST
dtDateTimeKst = with_tz(dtDateTimeUtc, tzone = "Asia/Seoul")
sDateTimeKst = format(dtDateTimeKst, "%Y%m%d%H%M")
# Set Output File Name
sOutputFileDirName = paste0("OUTPUT/", sFileDirNameSplit[2], "_", sFileDirNameSplit[4], "_", sDateTimeKst, ".csv")
# Set Figure File Name
sFigureType = "Ori"
sFigureFileDirName = paste0("FIG/", sFileDirNameSplit[2], "_", sFileDirNameSplit[4], "_", sFigureType, "_", sDateTimeKst, ".png")
-
NetCDF 파일 읽기
# NetCDF File Read
ncData = read.nc(ncFile)
# print(ncData)
netcdf klps_lc05_anal_201907030000 {
dimensions:
record = UNLIMITED ; // (1 currently)
n_valtimes = 1 ;
data_variables = 46 ;
namelen = 132 ;
charsPerLevel = 10 ;
x = 235 ;
y = 283 ;
levels_1 = 1 ;
levels_2 = 2 ;
levels_3 = 3 ;
levels_22 = 22 ;
levels_23 = 23 ;
variables:
float gh(record, levels_22, y, x) ;
gh:long_name = "Geopotential height" ;
gh:units = "m" ;
gh:udunits = "meters" ;
gh:uiname = "geoPotHt" ;
gh:valid_range = 0.f, 20000.f ;
gh:_FillValue = -99999.f ;
gh:_n3D = 22 ;
gh:levels = "MB 1100-50 by 50" ;
char ghLevels(levels_22, charsPerLevel) ;
char ghInventory(n_valtimes, levels_22) ;
float rh(record, levels_23, y, x) ;
rh:long_name = "Relative Humidity" ;
rh:units = "%" ;
rh:udunits = "percent" ;
rh:uiname = "rh" ;
rh:valid_range = 0.f, 100.f ;
rh:_FillValue = -99999.f ;
rh:_n3D = 23 ;
rh:levels = "SFC MB 1100-50 by 50" ;
char rhLevels(levels_23, charsPerLevel) ;
char rhInventory(n_valtimes, levels_23) ;
float t(record, levels_23, y, x) ;
t:long_name = "Temperature" ;
t:units = "K" ;
t:udunits = "degree_Kelvin" ;
t:uiname = "T" ;
t:valid_range = 180.f, 330.f ;
t:_FillValue = -99999.f ;
t:_n3D = 23 ;
t:levels = "SFC MB 1100-50 by 50" ;
char tLevels(levels_23, charsPerLevel) ;
char tInventory(n_valtimes, levels_23) ;
float uw(record, levels_23, y, x) ;
uw:long_name = "u wind component" ;
uw:units = "m/s" ;
uw:udunits = "meter/sec" ;
uw:uiname = "uWind" ;
uw:valid_range = -150.f, 150.f ;
uw:_FillValue = -99999.f ;
uw:_n3D = 23 ;
uw:levels = "SFC MB 1100-50 by 50" ;
char uwLevels(levels_23, charsPerLevel) ;
char uwInventory(n_valtimes, levels_23) ;
float vw(record, levels_23, y, x) ;
vw:long_name = "v wind component" ;
vw:units = "m/s" ;
vw:udunits = "meter/sec" ;
vw:uiname = "vWind" ;
vw:valid_range = -150.f, 150.f ;
vw:_FillValue = -99999.f ;
vw:_n3D = 23 ;
vw:levels = "SFC MB 1100-50 by 50" ;
char vwLevels(levels_23, charsPerLevel) ;
char vwInventory(n_valtimes, levels_23) ;
float pvv(record, levels_23, y, x) ;
pvv:long_name = "Pressure vertical velocity" ;
pvv:units = "Pa/s" ;
pvv:udunits = "pascal/second" ;
pvv:uiname = "Pvv" ;
pvv:valid_range = -2.5f, 2.5f ;
pvv:_FillValue = -99999.f ;
pvv:_n3D = 23 ;
pvv:levels = "SFC MB 1100-50 by 50" ;
char pvvLevels(levels_23, charsPerLevel) ;
char pvvInventory(n_valtimes, levels_23) ;
float sh(record, levels_22, y, x) ;
sh:long_name = "specific humidity" ;
sh:units = "kg/kg" ;
sh:udunits = "" ;
sh:uiname = "sh" ;
sh:valid_range = 0.f, 0.1f ;
sh:_FillValue = -99999.f ;
sh:_n3D = 22 ;
sh:levels = "MB 1100-50 by 50" ;
char shLevels(levels_22, charsPerLevel) ;
char shInventory(n_valtimes, levels_22) ;
float rr(record, levels_22, y, x) ;
rr:long_name = "LAPS radar reflectivity" ;
rr:units = "dBZ" ;
rr:udunits = "" ;
rr:uiname = "LapsRRef" ;
rr:valid_range = -20.f, 80.f ;
rr:_FillValue = -99999.f ;
rr:_n3D = 22 ;
rr:levels = "MB 1100-50 by 50" ;
char rrLevels(levels_22, charsPerLevel) ;
char rrInventory(n_valtimes, levels_22) ;
float ptyp(record, levels_22, y, x) ;
ptyp:long_name = "precipitation type" ;
ptyp:units = "none" ;
ptyp:udunits = "" ;
ptyp:uiname = "PrecipType" ;
ptyp:valid_range = 0.f, 20.f ;
ptyp:_FillValue = -99999.f ;
ptyp:_n3D = 22 ;
ptyp:levels = "MB 1100-50 by 50" ;
char ptypLevels(levels_22, charsPerLevel) ;
char ptypInventory(n_valtimes, levels_22) ;
float ctyp(record, levels_22, y, x) ;
ctyp:long_name = "cloud type" ;
ctyp:units = "none" ;
ctyp:udunits = "" ;
ctyp:uiname = "CldType" ;
ctyp:valid_range = 0.f, 20.f ;
ctyp:_FillValue = -99999.f ;
ctyp:_n3D = 22 ;
ctyp:levels = "MB 1100-50 by 50" ;
char ctypLevels(levels_22, charsPerLevel) ;
char ctypInventory(n_valtimes, levels_22) ;
float ccpc(record, levels_22, y, x) ;
ccpc:long_name = "fractional cloud cover pressure coord" ;
ccpc:units = "fractional" ;
ccpc:udunits = "" ;
ccpc:uiname = "FracCldCvr" ;
ccpc:valid_range = 0.f, 100.f ;
ccpc:_FillValue = -99999.f ;
ccpc:_n3D = 22 ;
ccpc:levels = "MB 1100-50 by 50" ;
char ccpcLevels(levels_22, charsPerLevel) ;
char ccpcInventory(n_valtimes, levels_22) ;
float cw(record, levels_22, y, x) ;
cw:long_name = "cloud liquid water" ;
cw:units = "grams/meter**3" ;
cw:udunits = "kilogram/meters3" ;
cw:uiname = "CldH20" ;
cw:valid_range = 0.f, 100.f ;
cw:_FillValue = -99999.f ;
cw:_n3D = 22 ;
cw:levels = "MB 1100-50 by 50" ;
char cwLevels(levels_22, charsPerLevel) ;
char cwInventory(n_valtimes, levels_22) ;
float cice(record, levels_22, y, x) ;
cice:long_name = "cloud ice" ;
cice:units = "grams/meter**3" ;
cice:udunits = "kilogram/meters3" ;
cice:uiname = "CldIce" ;
cice:valid_range = 0.f, 100.f ;
cice:_FillValue = -99999.f ;
cice:_n3D = 22 ;
cice:levels = "MB 1100-50 by 50" ;
char ciceLevels(levels_22, charsPerLevel) ;
char ciceInventory(n_valtimes, levels_22) ;
float hyc(record, levels_22, y, x) ;
hyc:long_name = "hydrometeor concentration" ;
hyc:units = "grams/meter**3" ;
hyc:udunits = "kilogram/meters3" ;
hyc:uiname = "HmConc" ;
hyc:valid_range = 0.f, 100.f ;
hyc:_FillValue = -99999.f ;
hyc:_n3D = 22 ;
hyc:levels = "MB 1100-50 by 50" ;
char hycLevels(levels_22, charsPerLevel) ;
char hycInventory(n_valtimes, levels_22) ;
float p(record, levels_2, y, x) ;
p:long_name = "pressure" ;
p:units = "Pa" ;
p:udunits = "pascal" ;
p:uiname = "atmP" ;
p:valid_range = 0.f, 110000.f ;
p:_FillValue = -99999.f ;
p:_n3D = 1 ;
p:levels = "SFC FH 1500" ;
char pLevels(levels_2, charsPerLevel) ;
char pInventory(n_valtimes, levels_2) ;
float mslp(record, levels_1, y, x) ;
mslp:long_name = "Mean Sea Level Pressure" ;
mslp:units = "Pa" ;
mslp:udunits = "pascal" ;
mslp:uiname = "MSL" ;
mslp:valid_range = 80000.f, 110000.f ;
mslp:_FillValue = -99999.f ;
mslp:_n3D = 0 ;
mslp:levels = "SFC" ;
char mslpLevels(levels_1, charsPerLevel) ;
char mslpInventory(n_valtimes, levels_1) ;
float dpt(record, levels_1, y, x) ;
dpt:long_name = "surface dewpoint temperature" ;
dpt:units = "K" ;
dpt:udunits = "degree_Kelvin" ;
dpt:uiname = "Td" ;
dpt:valid_range = 180.f, 330.f ;
dpt:_FillValue = -99999.f ;
dpt:_n3D = 0 ;
dpt:levels = "SFC" ;
char dptLevels(levels_1, charsPerLevel) ;
char dptInventory(n_valtimes, levels_1) ;
float pot(record, levels_1, y, x) ;
pot:long_name = "potential temperature" ;
pot:units = "K" ;
pot:udunits = "degree_Kelvin" ;
pot:uiname = "Tp" ;
pot:valid_range = 180.f, 330.f ;
pot:_FillValue = -99999.f ;
pot:_n3D = 0 ;
pot:levels = "SFC" ;
char potLevels(levels_1, charsPerLevel) ;
char potInventory(n_valtimes, levels_1) ;
float cc(record, levels_2, y, x) ;
cc:long_name = "cloud ceiling" ;
cc:units = "K" ;
cc:udunits = "degree_Kelvin" ;
cc:uiname = "CldCeil" ;
cc:valid_range = 0.f, 20000.f ;
cc:_FillValue = -99999.f ;
cc:_n3D = 0 ;
cc:levels = "SFC MSL" ;
char ccLevels(levels_2, charsPerLevel) ;
char ccInventory(n_valtimes, levels_2) ;
float tadv(record, levels_1, y, x) ;
tadv:long_name = "temperature advection" ;
tadv:units = "K/s" ;
tadv:udunits = "Kelvin/second" ;
tadv:uiname = "Tadv" ;
tadv:valid_range = -0.02f, 0.02f ;
tadv:_FillValue = -99999.f ;
tadv:_n3D = 0 ;
tadv:levels = "SFC" ;
char tadvLevels(levels_1, charsPerLevel) ;
char tadvInventory(n_valtimes, levels_1) ;
float ept(record, levels_1, y, x) ;
ept:long_name = "equivalent potential temperature" ;
ept:units = "K" ;
ept:udunits = "degree_Kelvin" ;
ept:uiname = "EqPot" ;
ept:valid_range = 180.f, 330.f ;
ept:_FillValue = -99999.f ;
ept:_n3D = 0 ;
ept:levels = "SFC" ;
char eptLevels(levels_1, charsPerLevel) ;
char eptInventory(n_valtimes, levels_1) ;
float ilw(record, levels_1, y, x) ;
ilw:long_name = "integrated liquid water" ;
ilw:units = "grams/meter**3" ;
ilw:udunits = "kilogram/meters3" ;
ilw:uiname = "IntLiqH20" ;
ilw:valid_range = 0.f, 2000.f ;
ilw:_FillValue = -99999.f ;
ilw:_n3D = 0 ;
ilw:levels = "SFC" ;
char ilwLevels(levels_1, charsPerLevel) ;
char ilwInventory(n_valtimes, levels_1) ;
float mcon(record, levels_1, y, x) ;
mcon:long_name = "moisture convergence" ;
mcon:units = "grams/meters**3" ;
mcon:udunits = "kilogram/meters3" ;
mcon:uiname = "moistConv" ;
mcon:valid_range = -0.01f, 0.01f ;
mcon:_FillValue = -99999.f ;
mcon:_n3D = 0 ;
mcon:levels = "SFC" ;
char mconLevels(levels_1, charsPerLevel) ;
char mconInventory(n_valtimes, levels_1) ;
float pta(record, levels_1, y, x) ;
pta:long_name = "potential temperature advection" ;
pta:units = "K/s" ;
pta:udunits = "Kelvin/second" ;
pta:uiname = "PotTadv" ;
pta:valid_range = -0.02f, 0.02f ;
pta:_FillValue = -99999.f ;
pta:_n3D = 0 ;
pta:levels = "SFC" ;
char ptaLevels(levels_1, charsPerLevel) ;
char ptaInventory(n_valtimes, levels_1) ;
float sli(record, levels_1, y, x) ;
sli:long_name = "lifted index" ;
sli:units = "K" ;
sli:udunits = "degree_Kelvin" ;
sli:uiname = "LftInd" ;
sli:valid_range = -20.f, 20.f ;
sli:_FillValue = -99999.f ;
sli:_n3D = 0 ;
sli:levels = "SFC" ;
char sliLevels(levels_1, charsPerLevel) ;
char sliInventory(n_valtimes, levels_1) ;
float ki(record, levels_1, y, x) ;
ki:long_name = "K index" ;
ki:units = "K" ;
ki:udunits = "degree_Kelvin" ;
ki:uiname = "KInd" ;
ki:valid_range = -20.f, 20.f ;
ki:_FillValue = -99999.f ;
ki:_n3D = 1 ;
ki:levels = "EA" ;
char kiLevels(levels_1, charsPerLevel) ;
char kiInventory(n_valtimes, levels_1) ;
float ttot(record, levels_1, y, x) ;
ttot:long_name = "Total totals index" ;
ttot:units = "K" ;
ttot:udunits = "degree_Kelvin" ;
ttot:uiname = "TotInd" ;
ttot:valid_range = -20.f, 20.f ;
ttot:_FillValue = -99999.f ;
ttot:_n3D = 1 ;
ttot:levels = "EA" ;
char ttotLevels(levels_1, charsPerLevel) ;
char ttotInventory(n_valtimes, levels_1) ;
float shwlt(record, levels_1, y, x) ;
shwlt:long_name = "Showalter index" ;
shwlt:units = "K" ;
shwlt:udunits = "degree_Kelvin" ;
shwlt:uiname = "TotInd" ;
shwlt:valid_range = -20.f, 20.f ;
shwlt:_FillValue = -99999.f ;
shwlt:_n3D = 1 ;
shwlt:levels = "850MB" ;
char shwltLevels(levels_1, charsPerLevel) ;
char shwltInventory(n_valtimes, levels_1) ;
float hidx(record, levels_1, y, x) ;
hidx:long_name = "Heat index" ;
hidx:units = "K" ;
hidx:udunits = "degree_Kelvin" ;
hidx:uiname = "HeatInd" ;
hidx:valid_range = 270.f, 330.f ;
hidx:_FillValue = -99999.f ;
hidx:_n3D = 0 ;
hidx:levels = "SFC" ;
char hidxLevels(levels_1, charsPerLevel) ;
char hidxInventory(n_valtimes, levels_1) ;
float cssi(record, levels_1, y, x) ;
cssi:long_name = "colorado severe storm index" ;
cssi:units = "none" ;
cssi:udunits = "" ;
cssi:uiname = "SevStInd" ;
cssi:valid_range = -20000.f, 20000.f ;
cssi:_FillValue = -99999.f ;
cssi:_n3D = 0 ;
cssi:levels = "SFC" ;
char cssiLevels(levels_1, charsPerLevel) ;
char cssiInventory(n_valtimes, levels_1) ;
float pbe(record, levels_1, y, x) ;
pbe:long_name = "positive buoyant energy" ;
pbe:units = "J/kg" ;
pbe:udunits = "joule/Kilogram" ;
pbe:uiname = "PosBuoyE" ;
pbe:valid_range = 0.f, 6000.f ;
pbe:_FillValue = -99999.f ;
pbe:_n3D = 0 ;
pbe:levels = "SFC" ;
char pbeLevels(levels_1, charsPerLevel) ;
char pbeInventory(n_valtimes, levels_1) ;
float nbe(record, levels_1, y, x) ;
nbe:long_name = "negative buoyant energy" ;
nbe:units = "J/kg" ;
nbe:udunits = "joule/Kilogram" ;
nbe:uiname = "NegBuoyE" ;
nbe:valid_range = 0.f, 6000.f ;
nbe:_FillValue = -99999.f ;
nbe:_n3D = 0 ;
nbe:levels = "SFC" ;
char nbeLevels(levels_1, charsPerLevel) ;
char nbeInventory(n_valtimes, levels_1) ;
float vis(record, levels_1, y, x) ;
vis:long_name = "visibility" ;
vis:units = "m" ;
vis:udunits = "meters" ;
vis:uiname = "Vis" ;
vis:_FillValue = -99999.f ;
vis:valid_range = 0.f, 100000.f ;
vis:_n3D = 0 ;
vis:levels = "SFC" ;
char visLevels(levels_1, charsPerLevel) ;
char visInventory(n_valtimes, levels_1) ;
float ccov(record, levels_1, y, x) ;
ccov:long_name = "LAPS cloud cover" ;
ccov:units = "%" ;
ccov:udunits = "percent" ;
ccov:uiname = "LapsCldCvr" ;
ccov:valid_range = 0.f, 100.f ;
ccov:_FillValue = -99999.f ;
ccov:_n3D = 0 ;
ccov:levels = "SFC" ;
char ccovLevels(levels_1, charsPerLevel) ;
char ccovInventory(n_valtimes, levels_1) ;
float cb(record, levels_1, y, x) ;
cb:long_name = "LAPS cloud base" ;
cb:units = "m" ;
cb:udunits = "meters" ;
cb:uiname = "LapsCldBase" ;
cb:valid_range = 0.f, 10000.f ;
cb:_FillValue = -99999.f ;
cb:_n3D = 0 ;
cb:levels = "SFC" ;
char cbLevels(levels_1, charsPerLevel) ;
char cbInventory(n_valtimes, levels_1) ;
float ctop(record, levels_1, y, x) ;
ctop:long_name = "LAPS cloud top" ;
ctop:units = "m" ;
ctop:udunits = "meters" ;
ctop:uiname = "LapsCldTop" ;
ctop:valid_range = 0.f, 50000.f ;
ctop:_FillValue = -99999.f ;
ctop:_n3D = 0 ;
ctop:levels = "SFC" ;
char ctopLevels(levels_1, charsPerLevel) ;
char ctopInventory(n_valtimes, levels_1) ;
float mret(record, levels_1, y, x) ;
mret:long_name = "maximum radar echo tops" ;
mret:units = "m" ;
mret:udunits = "meters" ;
mret:uiname = "MaxRdrEcho" ;
mret:valid_range = 0.f, 50000.f ;
mret:_FillValue = -99999.f ;
mret:_n3D = 0 ;
mret:levels = "SFC" ;
char mretLevels(levels_1, charsPerLevel) ;
char mretInventory(n_valtimes, levels_1) ;
float llr(record, levels_1, y, x) ;
llr:long_name = "low level reflectivity" ;
llr:units = "dBZ" ;
llr:udunits = "" ;
llr:uiname = "ReflectLL" ;
llr:valid_range = -20.f, 80.f ;
llr:_FillValue = -99999.f ;
llr:_n3D = 0 ;
llr:levels = "SFC" ;
char llrLevels(levels_1, charsPerLevel) ;
char llrInventory(n_valtimes, levels_1) ;
float heli(record, levels_1, y, x) ;
heli:long_name = "helicity" ;
heli:units = "m2/s2" ;
heli:udunits = "meter2/second2" ;
heli:uiname = "hel" ;
heli:valid_range = 0.f, 1000.f ;
heli:_FillValue = -99999.f ;
heli:_n3D = 0 ;
heli:levels = "SFC" ;
char heliLevels(levels_1, charsPerLevel) ;
char heliInventory(n_valtimes, levels_1) ;
float tpw(record, levels_1, y, x) ;
tpw:long_name = "integrated total precipitable water" ;
tpw:units = "m" ;
tpw:udunits = "meters" ;
tpw:uiname = "TotPrecipH20" ;
tpw:valid_range = 0.f, 0.01f ;
tpw:_FillValue = -99999.f ;
tpw:_n3D = 0 ;
tpw:levels = "SFC" ;
char tpwLevels(levels_1, charsPerLevel) ;
char tpwInventory(n_valtimes, levels_1) ;
float s1hr(record, levels_1, y, x) ;
s1hr:long_name = "LAPS 60 minute snow accum." ;
s1hr:units = "m" ;
s1hr:udunits = "meters" ;
s1hr:uiname = "LapsSnow60" ;
s1hr:valid_range = 0.f, 1.f ;
s1hr:_FillValue = -99999.f ;
s1hr:_n3D = 0 ;
s1hr:levels = "SFC" ;
char s1hrLevels(levels_1, charsPerLevel) ;
char s1hrInventory(n_valtimes, levels_1) ;
float stot(record, levels_1, y, x) ;
stot:long_name = "storm total snow accumulation" ;
stot:units = "m" ;
stot:udunits = "meters" ;
stot:uiname = "StmTotSnow" ;
stot:valid_range = 0.f, 10.f ;
stot:_FillValue = -99999.f ;
stot:_n3D = 0 ;
stot:levels = "SFC" ;
char stotLevels(levels_1, charsPerLevel) ;
char stotInventory(n_valtimes, levels_1) ;
float pc(record, levels_1, y, x) ;
pc:long_name = "LAPS 60 minute precip. accum." ;
pc:units = "m" ;
pc:udunits = "meters" ;
pc:uiname = "LapsPrecip60" ;
pc:valid_range = 0.f, 0.1f ;
pc:_FillValue = -99999.f ;
pc:_n3D = 0 ;
pc:levels = "SFC" ;
char pcLevels(levels_1, charsPerLevel) ;
char pcInventory(n_valtimes, levels_1) ;
float stpa(record, levels_1, y, x) ;
stpa:long_name = "storm total precip. accum." ;
stpa:units = "m" ;
stpa:udunits = "meters" ;
stpa:uiname = "StmTotPrecip" ;
stpa:valid_range = 0.f, 2.f ;
stpa:_FillValue = -99999.f ;
stpa:_n3D = 0 ;
stpa:levels = "SFC" ;
char stpaLevels(levels_1, charsPerLevel) ;
char stpaInventory(n_valtimes, levels_1) ;
float spt(record, levels_1, y, x) ;
spt:long_name = "LAPS precip type" ;
spt:units = "none" ;
spt:udunits = "" ;
spt:uiname = "LapsPrecipTyp" ;
spt:valid_range = 0.f, 20.f ;
spt:_FillValue = -99999.f ;
spt:_n3D = 0 ;
spt:levels = "SFC" ;
char sptLevels(levels_1, charsPerLevel) ;
char sptInventory(n_valtimes, levels_1) ;
float ptt(record, levels_1, y, x) ;
ptt:long_name = "LAPS precip type - LL radar ref modified" ;
ptt:units = "none" ;
ptt:udunits = "" ;
ptt:uiname = "LapsPrecipTypLL" ;
ptt:valid_range = 0.f, 7.f ;
ptt:_FillValue = -99999.f ;
ptt:_n3D = 0 ;
ptt:levels = "SFC" ;
char pttLevels(levels_1, charsPerLevel) ;
char pttInventory(n_valtimes, levels_1) ;
float fd(record, levels_1, y, x) ;
fd:long_name = "Fire Index" ;
fd:units = "none" ;
fd:udunits = "" ;
fd:uiname = "FInd" ;
fd:valid_range = 0.f, 20.f ;
fd:_FillValue = -99999.f ;
fd:_n3D = 0 ;
fd:levels = "SFC" ;
char fdLevels(levels_1, charsPerLevel) ;
char fdInventory(n_valtimes, levels_1) ;
float scp(record, levels_1, y, x) ;
scp:long_name = "snow cover percentage" ;
scp:units = "%" ;
scp:udunits = "percent" ;
scp:uiname = "SnCvrPer" ;
scp:valid_range = 0.f, 100.f ;
scp:_FillValue = -99999.f ;
scp:_n3D = 0 ;
scp:levels = "SFC" ;
char scpLevels(levels_1, charsPerLevel) ;
char scpInventory(n_valtimes, levels_1) ;
float smc(record, levels_3, y, x) ;
smc:long_name = "soil moisture" ;
smc:units = "m/m" ;
smc:udunits = "" ;
smc:uiname = "soilM" ;
smc:valid_range = 0.f, 1.f ;
smc:_FillValue = -99999.f ;
smc:_n3D = 0 ;
smc:levels = "BLS 1 2 3" ;
char smcLevels(levels_3, charsPerLevel) ;
char smcInventory(n_valtimes, levels_3) ;
int valtimeMINUSreftime(n_valtimes) ;
valtimeMINUSreftime:units = "seconds" ;
double valtime(record) ;
valtime:long_name = "valid time" ;
valtime:units = "seconds since (1970-1-1 00:00:00.0)" ;
double reftime ;
reftime:long_name = "reference time" ;
reftime:units = "seconds since (1970-1-1 00:00:00.0)" ;
char origin(namelen) ;
char model(namelen) ;
float staticTopo(y, x) ;
staticTopo:units = "meters" ;
staticTopo:long_name = "Topography" ;
staticTopo:_FillValue = -99999.f ;
float staticCoriolis(y, x) ;
staticCoriolis:units = "/second" ;
staticCoriolis:long_name = "Coriolis parameter" ;
staticCoriolis:_FillValue = -99999.f ;
float staticSpacing(y, x) ;
staticSpacing:units = "meters" ;
staticSpacing:long_name = "Grid spacing" ;
staticSpacing:_FillValue = -99999.f ;
// global attributes:
:cdlDate = "20020221" ;
:depictorName = "klps@19022779" ;
:projIndex = 3 ;
:projName = "LAMBERT_CONFORMAL" ;
:centralLat = 38.f ;
:centralLon = 126.f ;
:rotation = 30.f ;
:xMin = -0.04383551f ;
:xMax = 0.04585243f ;
:yMax = -0.614959f ;
:yMin = -0.7228721f ;
:lat00 = 31.34451f ;
:lon00 = 119.7993f ;
:latNxNy = 44.28384f ;
:lonNxNy = 133.6194f ;
:dxKm = 5.13378f ;
:dyKm = 5.125592f ;
:latDxDy = 38.02223f ;
:lonDxDy = 126.1543f ;
}
-
NetCDF 파일에서 1000 hPa의 지표면 온도 가져오기
# Get Variable
nTemp = ncData[["t"]][ , , 1]
-
Data Frame 형태로 초기화
# Set Data Frame
dfData = data.frame(c(nTemp)) %>%
dplyr::bind_cols(dfGeo) %>%
dplyr::rename(
nTemp = c.nTemp.
, nLon = lon
, nLat = lat
)
# Print Data Frame
dplyr::tbl_df(dfData)
-
Data Frame를 통해 L1 전처리
-
위도, 경도, 온도에 대한 최대값/최소값 설정
-
# L1 Processing Using Data Frame
dfDataL1 = dfData %>%
dplyr::filter(
dplyr::between(nLat, -90.0, 90.0)
, dplyr::between(nLon, -180.0, 180.0)
, dplyr::between(nTemp, 180.0, 330.0)
) %>%
dplyr::select(nLon, nLat, nTemp) %>%
dplyr::rename(nVal = nTemp)
# Print Data Frame
dplyr::tbl_df(dfDataL1)
-
다단계 B-Spline 근사화를 통해 L2 전처리
-
내삽을 위한 위/경도 및 해상도 설정 (dfNewData)
-
MBA::mba.points를 통해 1차원 변수 (위도, 경도) 내삽
-
자세한 내용은 상기 자료 처리 방안 및 활용 분석 기법에서 참조
-
#========================================================================
# 다단계 B-Spline 근사화 (Multi level B-Spline Approximation, MBA)
#========================================================================
sFigureType = "MBA"
# 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))
dfNewData = noncompliance::expand.grid.DT(
seq(from = xRange[1], to = xRange[2], by = 0.01)
, seq(from = yRange[1], to = yRange[2], by = 0.01)
, col.names = c("lon", "lat")
)
print(dfNewData)
dfInterData = MBA::mba.points(dfDataL1, dfNewData)
dfDataL2 = dfInterData %>%
as.data.frame() %>%
dplyr::rename(
xAxis = xyz.est.x
, yAxis = xyz.est.y
, zAxis = xyz.est.z
)
dplyr::tbl_df(dfDataL2)
-
가시화
-
데이터 세트에서 정보는 여전히 숨겨져 있기 때문에 가시화 필요
-
가시화는 일반적인 정적 탐색 데이터 분석에서 웹 브라우저의 동적 대화식 데이터 시각화에 이르기까지 다양함
-
특히 R의 기본 plot으로 여러 미학적 측면을 제어 할 수 있으나 Hadley Wickham (2016)이 개발한 ggplot2는 새로운 방법으로 시각화하기 때문에 이를 사용
-
-
가시화를 위한 초기 설정
-
컬러바 설정 (cbMatlab)
-
대한민국, 북한, 일본 Shp 파일 설정 (mapKor, mapPrk, mapJpn)
-
그림 파일명 설정 (sFigureFileDirName)
-
-
ggplot2를 이용한 지표면 온도 맵핑
# Set Value for Visualization
cbMatlab = colorRamps::matlab.like(11)
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")
sFigureFileDirName = paste0("FIG/", sFileDirNameSplit[2], "_", sFileDirNameSplit[4], "_", sFigureType, "_", sDateTimeKst, ".png")
# Map Visualization Using ggplot2
ggplot() +
coord_fixed(ratio = 1.1) +
theme_bw() +
geom_tile(data = dfDataL2, aes(x = xAxis, y = yAxis, fill = zAxis)) +
scale_fill_gradientn(colours = cbMatlab, limits=c(285, 305), breaks = seq(285, 305, 5), na.value = cbMatlab[length(cbMatlab)]) +
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 = sFigureFileDirName, width = 8, height = 10, dpi = 600)
-
정규 크리킹을 통해 L2 전처리
-
내삽을 위한 위/경도 및 해상도 설정 (spNewData)
-
gstat::krige를 통해 1차원 변수 (위도, 경도) 크리킹 내삽
-
이 과정에서 신규 격자를 기준으로 4개 자료 (nmax = 4)에 대해서만 수행
-
-
#=================================================
# Basic Operations-Krige Interpolation (Krige)
#=================================================
# set Spatial Points Data Frame
spData = dfDataL1
coordinates(spData) = ~ nLon + nLat
sFigureType = "Krige"
# 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 = noncompliance::expand.grid.DT(
seq(from = xRange[1], to = xRange[2], by = 0.01)
, seq(from = yRange[1], to = yRange[2], by = 0.01)
, col.names = c("lon", "lat")
)
coordinates(spNewData) = ~ lon + lat
gridded(spNewData) = TRUE
print(spNewData)
# fitting variogram model
variogram = automap::autofitVariogram(nVal ~ 1, spData)
plot(variogram)
# Apply Krige Model for the Data
dfInterData = gstat::krige(
formula = nVal ~ 1
, locations = spData
, newdata = spNewData
, model = variogram$var_model
, nmax = 4
)
dfDataL2 = dfInterData %>%
as.data.frame() %>%
dplyr::rename(
xAxis = lon
, yAxis = lat
, zAxis = var1.pred
)
dplyr::tbl_df(dfDataL2)
-
automap::autofitVariogram를 통해 최적의 베리오그램 모형 도출
-
ggplot2를 이용한 지표면 온도 맵핑
# Set Value for Visualization
cbMatlab = colorRamps::matlab.like2(11)
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")
sFigureFileDirName = paste0("FIG/", sFileDirNameSplit[2], "_", sFileDirNameSplit[4], "_", sFigureType, "_", sDateTimeKst, ".png")
# Map Visualization Using ggplot2
ggplot() +
coord_fixed(ratio = 1.1) +
theme_bw() +
geom_tile(data = dfDataL2, aes(x = xAxis, y = yAxis, fill = zAxis)) +
scale_fill_gradientn(colours = cbMatlab, limits=c(288, 300), breaks = seq(288, 300, 4), na.value = cbMatlab[length(cbMatlab)]) +
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 = sFigureFileDirName, width = 8, height = 10, dpi = 600)
[전체]
참고 문헌
[논문]
- 없음
[보고서]
- 없음
[URL]
- 없음
문의사항
[기상학/프로그래밍 언어]
- sangho.lee.1990@gmail.com
[해양학/천문학/빅데이터]
- saimang0804@gmail.com
본 블로그는 파트너스 활동을 통해 일정액의 수수료를 제공받을 수 있음
'프로그래밍 언어 > R' 카테고리의 다른 글
| [R] R을 이용한 수치해석 : 2020년 대학수학능력시험 (수능) 가형 기출문제 (0) | 2020.05.18 |
|---|---|
| [R] 서울대 통계 연구소 R을 이용한 빅데이터 분석 교육 연수 (0) | 2020.05.12 |
| [R] R을 이용한 통계 분석 및 데이터 시각화 : dplyr (2) | 2020.04.09 |
| [R] R을 이용한 통계 분석 및 데이터 시각화 : readr (0) | 2020.04.09 |
| [R] R을 이용한 통계 분석 및 데이터 시각화 : tidyr (0) | 2020.04.08 |
최근댓글