[R] "ncdf4" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (2)

정보
내용
소스 코드
참고 문헌
문의사항

정보

업무명 : "ncdf4" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (2)
작성자 : 이상호
작성일 : 2020-02-13
설 명 :
수정이력 :

내용

[개요]

안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
이전 포스팅에서 NetCDF 파일을 R로 직접 변환하고 데이터 프레임으로 변환하는 방법을 설명했습니다. 이 방법은 간단하나 단점이 있습니다. 즉 Netcdf 파일에서 배열의 잘못된 행렬을 읽을 수 없다는 것입니다. 이는 다른 행렬을 제거하여 배열의 1번째 행렬로부터 데이터 프레임을 변환합니다.

[R] NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (1)

정보 업무명 : NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (1) 작성자 : 이상호 작성일 : 2020-02-13 설 명 : 수정이력 : 내용 [개요] 안녕하세요? 기상..

shlee1990.tistory.com

이 포스팅는 접근을 확장하고 이전에서의 문제점을 해결하기 위한 것입니다. 즉 최근 자주 사용하는 패키지 (ncdf4)를 통해 NetCDF 파일을 데이터 프레임 형식으로 변환하는 과정을 소개해 드리고자 합니다.

[R] "RNetCDF" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (3)

정보 업무명 : "RNetCDF" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (3) 작성자 : 이상호 작성일 : 2020-02-13 설 명 : 수정이력 : 내용 [개..

shlee1990.tistory.com

크게 주요 3 단계로 구성하였습니다. 우선 NetCDF 파일에 포함 된 메타 데이터를 읽고 파일에 저장된 데이터구조를 확인하는 방법을 보여 드리겠습니다. 2 단계의 경우 데이터를 추출하고 마지막 (3 단계)에서는 데이터 프레임으로 변환 및 가시화하겠습니다.

[특징]

R에서 NetCDF 파일을 처리하기 위해서 데이터 프레임으로 변환 처리가 요구되며 이 프로그램은 이러한 목적을 달성하기 위한 소프트웨어

[기능]

ncdf4 패키지를 통해 파일 읽기
배열을 데이터 프레임으로 변환
해수면 온도 가시화

[활용 자료]

자료명 : sst.wkmean.1990-present.nc
자료 종류 : 최적 내삽 해수면 온도
확장자 : NetCDF
영역 : 전지구
기간 : 1989년 12월 31일 - 2020년 02월 20일
시간 해상도 : 일 1개 (21시)
제공처 : ESRL | Physical sciences Division

ESRL : PSD : NOAA Optimum Interpolation (OI) Sea Surface Temperature (SST) V2

www.esrl.noaa.gov

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

없음

[사용법]

소스 코드 참조

[사용 OS]

Windows 10

[사용 언어]

R v3.6.2
R Studio v1.2.5033

소스 코드

[명세]

전역 설정
- 최대 10 자리 설정
- 메모리 해제
- 영어 인코딩 설정
- 폰트 설정

# Set Option
memory.limit(size = 9999999999999)
options(digits = 10)
Sys.setlocale("LC_TIME", "english")
font = "Palatino Linotype"

라이브러리 읽기

# Library Load
library(extrafont)
library(ncdf4)
library(tidyverse)
library(ncdump)
library(RNetCDF)
library(tidyverse)
library(lubridate)
library(gganimate)
library(insol)
library(spData)
library(raster)

메타 데이터 이해
- ncdump 패키지에서 NetCDF 함수를 사용하여 파일의 메타 데이터를 탐색
- 이러한 메타 데이터의 데이터 그룹 중 dimension에는 위도, 경도 및 시간과 같은 벡터 형식 변수가 포함
- 추가로 메타 데이터에서 variable 배열입니다. 배열은 해수면 온도 (sst)를 행렬로서 포함 하였다. 즉 매일마다 데이터이므로 각 행렬은 의미합니다.

sFileDirName = Sys.glob("INPUT/sst.wkmean.1990-present.nc")

liNcInfo = ncdump::NetCDF(sFileDirName)

liNcInfo

> liNcInfo

$dimension
# A tibble: 4 x 7
  name    len unlim group_index group_id    id create_dimvar
  <chr> <int> <lgl>       <int>    <int> <int> <lgl>        
1 lat     180 FALSE           1   196608     0 TRUE         
2 lon     360 FALSE           1   196608     1 TRUE         
3 time   1574 TRUE            1   196608     2 TRUE         
4 nbnds     2 FALSE           1   196608     3 FALSE        

$unlimdims
# A tibble: 1 x 2
     id units                       
  <int> <chr>                       
1     2 days since 1800-1-1 00:00:00

$dimvals
# A tibble: 2,116 x 2
      id  vals
   <int> <dbl>
 1     0  89.5
 2     0  88.5
 3     0  87.5
 4     0  86.5
 5     0  85.5
 6     0  84.5
 7     0  83.5
 8     0  82.5
 9     0  81.5
10     0  80.5
# ... with 2,106 more rows

$groups
# A tibble: 1 x 6
      id name  ndims nvars natts fqgn 
   <int> <chr> <int> <int> <int> <chr>
1 196608 ""        4     5    11 ""   

$file
# A tibble: 1 x 10
  filename          writable     id safemode format    is_GMT ndims natts unlimdimid nvars
  <chr>             <lgl>     <int> <lgl>    <chr>     <lgl>  <dbl> <dbl>      <dbl> <dbl>
1 INPUT/sst.wkmean~ FALSE    196608 FALSE    NC_FORMA~ FALSE      4    11          3     2

$variable
# A tibble: 2 x 18
  name  ndims natts prec  units longname group_index storage shuffle compression unlim
  <chr> <int> <int> <chr> <chr> <chr>          <int>   <dbl> <lgl>   <lgl>       <lgl>
1 sst       3    16 short "deg~ Weekly ~           1       2 FALSE   NA          TRUE 
2 time~     2     1 doub~ ""    Time Bo~           1       2 FALSE   NA          TRUE 
# ... with 7 more variables: make_missing_value <lgl>, missval <dbl>, hasAddOffset <lgl>,
#   addOffset <dbl>, hasScaleFact <lgl>, scaleFact <dbl>, id <dbl>

$vardim
# A tibble: 5 x 2
     id dimids
  <dbl>  <int>
1     2      1
2     2      0
3     2      2
4     4      3
5     4      2

$attribute
[1] "NetCDF attributes:"
[1] "Global"
[1] "\n"
# A tibble: 1 x 11
  title Conventions history comments platform source institution References NCO  
  <chr> <chr>       <chr>   <chr>    <chr>    <chr>  <chr>       <chr>      <chr>
1 NOAA~ CF-1.0      Create~ "Data d~ Model    NCEP ~ National C~ https://w~ 4.0.0
# ... with 2 more variables: dataset_title <chr>, source_url <chr>
[1] "\n"
[1] "Variable attributes:"
[1] "variable attributes: sst"

attr(,"class")
[1] "NetCDF" "list"

NetCDF 파일 설정 및 읽기
- ncdf4::nc_open를 통해 파일 읽기 및 각종 변수에 대한 세부 정보를 확인

liNcOpen = ncdf4::nc_open(sFileDirName)

liNcOpen

> liNcOpen

File INPUT/sst.wkmean.1990-present.nc (NC_FORMAT_CLASSIC):

     2 variables (excluding dimension variables):
        short sst[lon,lat,time]   
            long_name: Weekly Mean of Sea Surface Temperature
            unpacked_valid_range: -5
             unpacked_valid_range: 40
            actual_range: -1.79999995231628
             actual_range: 36.1599998474121
            units: degC
            add_offset: 0
            scale_factor: 0.00999999977648258
            missing_value: 32767
            precision: 2
            least_significant_digit: 2
            var_desc: Sea Surface Temperature
            dataset: NOAA Optimum Interpolation (OI) SST V2
            level_desc: Surface
            statistic: Weekly Mean
            parent_stat: Individual obs
            standard_name: sea_surface_temperature
            valid_range: -500
             valid_range: 4000
        double time_bnds[nbnds,time]   
            long_name: Time Boundaries

     4 dimensions:
        lat  Size:180
            units: degrees_north
            long_name: Latitude
            actual_range: 89.5
             actual_range: -89.5
            standard_name: latitude
            axis: Y
            coordinate_defines: center
        lon  Size:360
            units: degrees_east
            long_name: Longitude
            actual_range: 0.5
             actual_range: 359.5
            standard_name: longitude
            axis: X
            coordinate_defines: center
        time  Size:1574   *** is unlimited ***
            units: days since 1800-1-1 00:00:00
            long_name: Time
            actual_range: 69395
             actual_range: 80406
            delta_t: 0000-00-07 00:00:00
            avg_period: 0000-00-07 00:00:00
            standard_name: time
            axis: T
            bounds: time_bnds
        nbnds  Size:2
[1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named nbnds BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"

    11 global attributes:
        title: NOAA Optimum Interpolation (OI) SST V2
        Conventions: CF-1.0
        history: Created 10/2002 by RHS
        comments: Data described in  Reynolds, R.W., N.A. Rayner, T.M.
Smith, D.C. Stokes, and W. Wang, 2002: An Improved In Situ and Satellite
SST Analysis for Climate, J. Climate
        platform: Model
        source: NCEP Climate Modeling Branch
        institution: National Centers for Environmental Prediction
        References: https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html
        NCO: 4.0.0
        dataset_title: NOAA Optimum Interpolation (OI) SST V2
        source_url: http://www.emc.ncep.noaa.gov/research/cmb/sst_analysis/

변수 추출
- ncdf4::ncvar_get를 통해 벡터 및 배열 형식으로 변수 (위도, 경도, 시간, 해수면 온도) 추출

# spatial components
nLon = ncdf4::ncvar_get(liNcOpen, "lon")
nLat = ncdf4::ncvar_get(liNcOpen, "lat")

# temporal component
nTime = ncdf4::ncvar_get(liNcOpen, "time")

# sea surface temperature
nVal = ncdf4::ncvar_get(liNcOpen, "sst")

시간 변환
- 시간 변수 (nTime)는 줄리안 데이 (Julian Day)로서 날짜의 시작 나타냄
- 그러나 메타 데이터에는 줄리안 데이를 그레고리안 날짜 (Gregorian Calendar)로 변환하는 데 사용할 수있는 시작 데이터 정보가 제공 (days since 1800-1-1 00:00:00).
- 원래 날짜를 알면이 요일을 캘린더로 변환 할 수 있습니다. 메타 데이터를 살펴보면 캘린더의 원래 날짜가로 지정된 것으로 나타났습니다 julian_day_unit: days since 1950-01-01 00:00:00. 여기 또 다른 도전이 있습니다.
- 즉 시간 변수는 줄리안 데이이나 원래 날짜는 그레고안 날짜형식임
- 따라서 시간을 공통 형식으로 표준화 필요
  - insol::JDymd를 통해 시작 시간 (1800-1-1 00:00:00)을 줄리안 데이터로 변환
  - 시간 변수에 시작 시간 추가하고 insol::JD를 통해 줄리안 데이를 그레고리안 날짜로 변환

# "days since 1800-1-1 00:00:00"
liNcInfo$unlimdims$units

# convert time original (to) to julian 
nTo = insol::JDymd(year = 1800, month = 1, day = 1)

# add the original time to the extracted time
nJd = nTo + nTime

#convert the julian day to gregorian calender
dtDate = insol::JD(nJd, inverse = TRUE)

dplyr::tbl_df(
   data.frame(nTime, nJd, dtDate)
)

데이터 프레임 (Data Frame) 변환
- 배열의 행렬을 데이터 프레임으로 변환할 때 반복문 및 noncompliance::expand.grid.DT 함수를 이용할 수 있음

반복문
- 시간 변수에 대한 반복문을 통해 위/경도 및 해수면 온도를 추가
- 이 과정에서 시간, 위도, 경도, 해수면 온도 순으로 변환되나 다소 오랜 시간 소요

dfData = NULL

for (iCount in 1:length(dtDate)) {
# for (iCount in 1:1) {
   
   dfTmpData = data.frame(
      nLon = c(nLon)
      , nLat = rep(nLat, each = length(nLon))
      , nVal = c(nVal[ , ,iCount])
      , dtDate = dtDate[iCount]
   )
   
      dfData = dplyr::bind_rows(dfData, dfTmpData)
}


dplyr::tbl_df(dfData)

noncompliance::expand.grid.DT 함수
- 앞서 반복문과 동일한 결과를 반환될 뿐만 아니라 벡터화를 사용하기 때문에 속도 개선

iTime = length(dtDate)
# iTime = 1

dfData = data.frame(
   noncompliance::expand.grid.DT(dtDate[1:iTime], nLat, nLon)
   , c(nVal[ , , 1:iTime])
)

colnames(dfData) = c("dtDate", "nLat", "nLon", "nVal")

dplyr::tbl_df(dfData)

Data Frame을 이용한 L1 전처리
- 전체 기간 (1989년 12월 31일 - 2020년 02월 20일), 위도, 경도에 따른 평균 수행

# L1 Processing Using Data Frame
dfDataL1 = dfData %>%
   dplyr::group_by(nLon, nLat) %>%
   dplyr::summarise(nMeanVal = mean(nVal, na.rm = TRUE))

dplyr::tbl_df(dfDataL1)

Data Frame L1을 이용한 L2 전처리
- 위도 변수 (nLon)가 0~360을 -180~180으로 변환

# L2 Processing Using Data Frame L1
dfDataL2 = dfDataL1 %>%
   dplyr::mutate(nLon180 = metR::ConvertLongitude(nLon, from = 360))

summary(dfDataL2)

가시화
- 데이터 세트에서 정보는 여전히 숨겨져 있기 때문에 가시화 필요
- 가시화는 일반적인 정적 탐색 데이터 분석에서 웹 브라우저의 동적 대화식 데이터 시각화에 이르기까지 다양함
- 특히 R의 기본 plot으로 여러 미학적 측면을 제어 할 수 있으나 Hadley Wickham (2016)이 개발한 ggplot2는 새로운 방법으로 시각화하기 때문에 이를 사용

가시화를 위한 초기 설정
ggplot2를 이용한 가시화

# Visualization Using ggplot2
ggplot() +
   theme_bw() +
   geom_tile(data = dfDataL2, aes(x = nLon, y = nLat, fill = nMeanVal)) +
   metR::geom_text_contour(data = dfDataL2, aes(x = nLon, y = nLat, z = nMeanVal), stroke = 0.2, check_overlap = TRUE, rotate = TRUE, na.rm = TRUE) +
   metR::geom_contour2(data = dfDataL2, aes(x = nLon, y = nLat, z = nMeanVal), color = "black", alpha = 0.3) +
   scale_fill_gradientn(colours = cbOcean, limits=c(-5, 35), breaks = seq(-5, 35, 10), na.value = cbOcean[length(cbOcean)]) +
   geom_sf(data = mapData, fill = "grey100", col = "black") +
   metR::scale_x_longitude(expand = c(0, 0), breaks = seq(-180, 180, 60), limits = c(-180, 180)) +
   metR::scale_y_latitude(expand = c(0, 0), breaks = seq(-90, 90, 30), limits = c(-90, 90)) +
   labs(
      x = ""
      , y = ""
      , fill = "Mean Sea Surface Temperature [℃]"
      , colour = ""
      , title  = "NOAA Optimum Interpolation (OI) SST V2"
      , subtitle = "Period : December 31, 1989 - February 20, 2020"
      , caption = "Source : NCEP Climate Modeling Branch"
   ) +
   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.title = element_text(face = "bold", size = 14, colour = "white")
      , legend.position = c(0, 1)
      , legend.justification = c(0, 0.96)
      , legend.key = element_blank()
      , legend.text = element_text(size = 14, face = "bold", colour = "white")
      , legend.background = element_blank()
      , text = element_text(family = font)
      , plot.margin = unit(c(0, 8, 0, 0), "mm")
   ) +
   ggsave(filename = paste0("FIG2/Mean_Sea_Surface_Temperature.png"), width = 12, height = 8, dpi = 600)

[전체]

	#======================================================================================================================================
	# Routine : Main R program
	#
	# Purpose : Visualization and Conversion Processing into Data Frame Using NOAA Optimum Interpolation Data in NetCDF Format (2)
	#
	# Author : MS. Sang-Ho Lee
	#
	# Revisions: V1.0 March 07, 2019 First release (MS. Sang-Ho Lee)
	#======================================================================================================================================

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

	# Library Load
	library(extrafont)
	library(ncdf4)
	library(tidyverse)
	library(ncdump)
	library(RNetCDF)
	library(tidyverse)
	library(lubridate)
	library(gganimate)
	library(insol)
	library(spData)
	library(raster)

	sFileDirName = Sys.glob("INPUT/sst.wkmean.1990-present.nc")

	liNcInfo = ncdump::NetCDF(sFileDirName)

	liNcInfo

	liNcOpen = ncdf4::nc_open(sFileDirName)

	liNcOpen

	# spatial components
	nLon = ncdf4::ncvar_get(liNcOpen, "lon")
	nLat = ncdf4::ncvar_get(liNcOpen, "lat")

	# temporal component
	nTime = ncdf4::ncvar_get(liNcOpen, "time")

	# sea surface temperature
	nVal = ncdf4::ncvar_get(liNcOpen, "sst")


	# "days since 1800-1-1 00:00:00"
	liNcInfo$unlimdims$units

	# convert time original (to) to julian
	nTo = insol::JDymd(year = 1800, month = 1, day = 1)

	# add the original time to the extracted time
	nJd = nTo + nTime

	#convert the julian day to gregorian calender
	dtDate = insol::JD(nJd, inverse = TRUE)

	dplyr::tbl_df(
	data.frame(nTime, nJd, dtDate)
	)


	dfData = NULL

	# for (iCount in 1:length(dtDate)) {
	for (iCount in 1:1) {

	dfTmpData = data.frame(
	dtDate = dtDate[iCount]
	, nLat = rep(nLat, each = length(nLon))
	, nLon = c(nLon)
	, nVal = c(nVal[ , ,iCount])
	)

	dfData = dplyr::bind_rows(dfData, dfTmpData)
	}


	dplyr::tbl_df(dfData)


	iTime = length(dtDate)
	# iTime = 1

	dfData = data.frame(
	noncompliance::expand.grid.DT(dtDate[1:iTime], nLat, nLon)
	, c(nVal[ , , 1:iTime])
	)

	colnames(dfData) = c("dtDate", "nLat", "nLon", "nVal")

	dplyr::tbl_df(dfData)

	# L1 Processing Using Data Frame
	dfDataL1 = dfData %>%
	dplyr::group_by(nLon, nLat) %>%
	dplyr::summarise(nMeanVal = mean(nVal, na.rm = TRUE))

	dplyr::tbl_df(dfDataL1)

	# L2 Processing Using Data Frame L1
	dfDataL2 = dfDataL1 %>%
	dplyr::mutate(nLon180 = metR::ConvertLongitude(nLon, from = 360))

	summary(dfDataL2)


	# Set Value for Visualization
	font = "Palatino Linotype"
	mapData = spData::world
	cbOcean = oce::oceColors9A(120)

	# Visualization Using ggplot2
	ggplot() +
	theme_bw() +
	geom_tile(data = dfDataL2, aes(x = nLon180, y = nLat, fill = nMeanVal)) +
	metR::geom_text_contour(data = dfDataL2, aes(x = nLon180, y = nLat, z = nMeanVal), stroke = 0.2, check_overlap = TRUE, rotate = TRUE, na.rm = TRUE) +
	metR::geom_contour2(data = dfDataL2, aes(x = nLon180, y = nLat, z = nMeanVal), color = "black", alpha = 0.3) +
	scale_fill_gradientn(colours = cbOcean, limits=c(-5, 35), breaks = seq(-5, 35, 10), na.value = cbOcean[length(cbOcean)]) +
	geom_sf(data = mapData, fill = "grey100", col = "black") +
	metR::scale_x_longitude(expand = c(0, 0), breaks = seq(-180, 180, 60), limits = c(-180, 180)) +
	metR::scale_y_latitude(expand = c(0, 0), breaks = seq(-90, 90, 30), limits = c(-90, 90)) +
	labs(
	x = ""
	, y = ""
	, fill = "Mean Sea Surface Temperature [℃]"
	, colour = ""
	, title = "NOAA Optimum Interpolation (OI) SST V2"
	, subtitle = "Period : December 31, 1989 - February 20, 2020"
	, caption = "Source : NCEP Climate Modeling Branch"
	) +
	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.title = element_text(face = "bold", size = 14, colour = "white")
	, legend.position = c(0, 1)
	, legend.justification = c(0, 0.96)
	, legend.key = element_blank()
	, legend.text = element_text(size = 14, face = "bold", colour = "white")
	, legend.background = element_blank()
	, text = element_text(family = font)
	, plot.margin = unit(c(0, 8, 0, 0), "mm")
	) +
	ggsave(filename = paste0("FIG2/Mean_Sea_Surface_Temperature.png"), width = 12, height = 8, dpi = 600)

view raw Visualization_and_Conversion_Processing_into_Data_Frame_Using_NOAA_Optimum_Interpolation_Data_in_NetCDF_Format_(2).R hosted with ❤ by GitHub

참고 문헌

[논문]

없음

[보고서]

없음

[URL]

없음

문의사항

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

sangho.lee.1990@gmail.com

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

saimang0804@gmail.com

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

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

[R] 육지 및 해양에 대한 마스킹 (Masking) 처리 및 가시화 (0)	2020.03.08
[R] "RNetCDF" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (3) (0)	2020.03.08
[R] "raster" 패키지 및 NetCDF 형식인 NOAA 최적 내삽 자료를 이용하여 데이터 프레임 (Data Frame)으로 변환 처리 및 가시화 (1) (0)	2020.03.07
[R] sf 패키지를 이용한 지도 투영법 소개 및 가시화 (0)	2020.03.06
[R] 시/공간 자료 처리 및 분석 (Handling and Analyzing Spatio-Temporal Data) 패키지 소개 (0)	2020.03.05

내 블로그 - 관리자 홈 전환	`Q` `Q`
새 글 쓰기	`W` `W`

글 수정 (권한 있는 경우)	`E` `E`
댓글 영역으로 이동	`C` `C`

이 페이지의 URL 복사	`S` `S`
맨 위로 이동	`T` `T`
티스토리 홈 이동	`H` `H`
단축키 안내	`Shift` + `/` `⇧` + `/`