반응형

     정보

    • 업무명     : NetCDF 형식인 NPP/CERES SSF 1deg 기상위성 자료를 이용한 가시화

    • 작성자     : 이상호

    • 작성일     : 2019-12-30

    • 설   명      :

    • 수정이력 :

     

     내용

    [특징]

    • NetCDF 형태인 기상위성 자료를 이해하기 위해서 가시화 도구가 필요하며 이 프로그램은 이러한 목적을 달성하기 위한 소프트웨어

     

    [기능]

    • NPP/CERES SSF 1deg 기상위성 자료를 이용하여 전처리

    • 2015년 07월 - 2016년 12월에 대한 요약 통계량 (평균, 최대값, 최소값, 자료 개수) 및 평균 강수량 가시화

     

    [활용 자료]

    • 위성명 : NPP 기상위성

    • 센서명 : CERES

    • 자료 레벨 : SSF 1deg

    • 자료 종류 : 날짜, 시간, 위도, 경도, 강수량, 지표면 특성, 운량

    • 영역 : 호주

    • 해상도 : 100 km

    • 확장자 : NetCDF

    • 기간 : 2015년 07월 15일 - 2016년 12월 21일

     

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

    • 없음

     

    [사용법]

    • 입력 자료를 동일 디렉터리 위치

    • 소스 코드를 실행 (Rscript Processing_Using_NetCDF_Format_of_NPP_CERES_SSF_1deg_Data.R)

    • 가시화 결과 확인

     

    [사용 OS]

    • Windows10

     

    [사용 언어]

    • R v3.6.2

    • R Studio v1.2.5033

     

     소스 코드

    [명세]

    • 메모리 해제

    # Set Option
    memory.limit(size = 9999999999999)

     

    • 라이브러리 읽기

    # Library Load
    library(RNetCDF)
    library(tidyverse)
    library(metR)
    library(colorRamps)
    library(ggrepel)
    library(extrafont)
    library(sf)

     

    • 파일 읽기

    # File Read
    sInputFileDirName = "INPUT_SSF1deg/*.nc"
    sFileDirName = Sys.glob(sInputFileDirName)
    
    length(sFileDirName)
    sFileDirName

     

     

    • NetCDF 파일 읽기

    # NetCDF File Open
    ncFile = open.nc(sFileDirName[iCount])

     

    • 출력을 위한 파일명 초기화

    # Set Output File Name
    sFileDirNameSplit = unlist(str_split(string = sFileDirName[iCount], pattern = "_|/|\\."))
    sOutputFileDirName = paste0('OUTPUT/', sFileDirNameSplit[5], "_", sFileDirNameSplit[7], ".OUT")

     

     

    • NetCDF 파일 읽기

    # NetCDF File Read
    ncData = read.nc(ncFile)
    
    # print.nc(ncFile)
    netcdf classic {
    dimensions:
    	time = UNLIMITED ; // (29240974 currently)
    	The_8_most_prevalent_surface_types = 8 ;
    	Conditions_clear__lower__upper__upper_over_lower = 4 ;
    variables:
    	NC_DOUBLE time(time) ;
    		NC_CHAR time:long_name = "time" ;
    		NC_CHAR time:units = "days since 1970-01-01 00:00:00" ;
    		NC_DOUBLE time:_FillValue = 1.79769313486232e+308 ;
    		NC_DOUBLE time:valid_range = 0, 39412.5 ;
    	NC_FLOAT lon(time) ;
    		NC_CHAR lon:long_name = "longitude" ;
    		NC_CHAR lon:units = "degrees_east" ;
    		NC_FLOAT lon:_FillValue = 3.40282346638529e+38 ;
    		NC_FLOAT lon:valid_range = -180, 180 ;
    	NC_FLOAT lat(time) ;
    		NC_CHAR lat:long_name = "latitude" ;
    		NC_CHAR lat:units = "degrees_north" ;
    		NC_FLOAT lat:_FillValue = 3.40282346638529e+38 ;
    		NC_FLOAT lat:valid_range = -90, 90 ;
    	NC_DOUBLE Time_of_observation(time) ;
    		NC_CHAR Time_of_observation:orig_name = "Time of observation" ;
    		NC_CHAR Time_of_observation:units = "day" ;
    		NC_CHAR Time_of_observation:format = "F18.9" ;
    		NC_DOUBLE Time_of_observation:_FillValue = 1.79769313486232e+308 ;
    		NC_DOUBLE Time_of_observation:valid_range = 2440000, 2480000 ;
    	NC_FLOAT Longitude_of_CERES_FOV_at_surface(time) ;
    		NC_CHAR Longitude_of_CERES_FOV_at_surface:orig_name = "Longitude of CERES FOV at surface" ;
    		NC_CHAR Longitude_of_CERES_FOV_at_surface:units = "degrees" ;
    		NC_CHAR Longitude_of_CERES_FOV_at_surface:format = "F18.9" ;
    		NC_FLOAT Longitude_of_CERES_FOV_at_surface:_FillValue = 3.40282346638529e+38 ;
    		NC_FLOAT Longitude_of_CERES_FOV_at_surface:valid_range = 0, 360 ;
    	NC_FLOAT Colatitude_of_CERES_FOV_at_surface(time) ;
    		NC_CHAR Colatitude_of_CERES_FOV_at_surface:orig_name = "Colatitude of CERES FOV at surface" ;
    		NC_CHAR Colatitude_of_CERES_FOV_at_surface:units = "degrees" ;
    		NC_CHAR Colatitude_of_CERES_FOV_at_surface:format = "F18.9" ;
    		NC_FLOAT Colatitude_of_CERES_FOV_at_surface:_FillValue = 3.40282346638529e+38 ;
    		NC_FLOAT Colatitude_of_CERES_FOV_at_surface:valid_range = 0, 180 ;
    	NC_SHORT Surface_type_index(The_8_most_prevalent_surface_types, time) ;
    		NC_CHAR Surface_type_index:orig_name = "Surface type index" ;
    		NC_CHAR Surface_type_index:units = "N/A" ;
    		NC_CHAR Surface_type_index:format = "I10" ;
    		NC_SHORT Surface_type_index:_FillValue = 32767 ;
    		NC_SHORT Surface_type_index:valid_range = 1, 20 ;
    	NC_FLOAT Precipitable_water(time) ;
    		NC_CHAR Precipitable_water:orig_name = "Precipitable water" ;
    		NC_CHAR Precipitable_water:units = "centimeters" ;
    		NC_CHAR Precipitable_water:format = "F18.9" ;
    		NC_FLOAT Precipitable_water:_FillValue = 3.40282346638529e+38 ;
    		NC_FLOAT Precipitable_water:valid_range = 0.00100000004749745, 10 ;
    	NC_FLOAT Clear_layer_overlap_percent_coverages(Conditions_clear__lower__upper__upper_over_lower, time) ;
    		NC_CHAR Clear_layer_overlap_percent_coverages:orig_name = "Clear/layer/overlap percent coverages" ;
    		NC_CHAR Clear_layer_overlap_percent_coverages:units = "N/A" ;
    		NC_CHAR Clear_layer_overlap_percent_coverages:format = "F18.9" ;
    		NC_FLOAT Clear_layer_overlap_percent_coverages:_FillValue = 3.40282346638529e+38 ;
    		NC_FLOAT Clear_layer_overlap_percent_coverages:valid_range = 0, 100 ;
    
    // global attributes:
    		NC_CHAR :Conventions = "CF-1.0" ;
    		NC_CHAR :Subsetter_title = "ASDC CERES Subset" ;
    		NC_CHAR :Subsetter_version = "2.9.b1" ;
    		NC_CHAR :Subsetter_institution = "Atmospheric Science Data Center (ASDC) http://eosweb.larc.nasa.gov" ;
    		NC_CHAR :Subsetter_history = "2017-12-29T02:00:57 -0500 SubsetCeresSsf" ;
    		NC_CHAR :Subsetter_temporalFilter = "2015-07-15T00:00:00.000000Z to 2017-02-15T23:59:59.999999Z" ;
    		NC_CHAR :Subsetter_spatialFilter = "POLYGON ((112.6593017578125 -44.285888671875, 112.6593017578125 -10.535888671875, 153.9678955078125 -10.535888671875, 153.9678955078125 -44.285888671875, 112.6593017578125 -44.285888671875))" ;
    		NC_CHAR :Subsetter_parameterFilter = "none" ;
    		NC_CHAR :platform = "NPP" ;
    		NC_CHAR :history = "Fri Dec 29 02:17:45 2017: ncrcat -o CER_SSF_NPP-FM5-VIIRS_Edition1A_2015071503-2016122118_Subset.nc" ;
    		NC_INT :nco_input_file_number = 4151 ;
    		NC_CHAR :nco_input_file_list = "CER_SSF_NPP-FM5-VIIRS_Edition1A_100101.2015071503_Subset.nc
    

     

    • NetCDF 파일 갯수 및 Time 단위 가져오기

    # Get Time Unit
    nFileNumber = dim.inq.nc(ncFile, "time")$length
    sTimeUnit = att.get.nc(ncFile, "time", "units")

     

     

    • NetCDF 변수 가져오기

    # Get Variable
    nDateTime = utcal.nc(sTimeUnit, ncData$time)
    nLon  = ncData[["lon"]]
    nLat  = ncData[["lat"]]
    nTpw  = ncData[["Precipitable_water"]]
    nSurfaceType = ncData[["Surface_type_index"]][1, ]
    nClearFraction = ncData[["Clear_layer_overlap_percent_coverages"]][1, ]

     

    • NetCDF 변수들을 Data Frame 형태로 초기화

    # Set Data Frame
    dfData = data.frame(nDateTime, nLat, nLon, nTpw, nSurfaceType, nClearFraction)
       
    # Print Data Frame
    dplyr::tbl_df(dfData)

     

     

    • Data Frame를 통해 L1 전처리

      • 각 변수에 따라 최소값/최대값 설정

    # L1 Processing Using Data Frame
    dfDataL1 = dfData %>%
       dplyr::filter(
          between(nLat, -90.0, 90.0)
          , between(nLon, -180.0, 180.0)
          , between(nTpw, 0, 10)
          , between(nSurfaceType, 1, 20)
          , between(nClearFraction, 0, 100)
       )
       
    # Print Data Frame
    dplyr::tbl_df(dfDataL1)

     

     

    • L1 자료를 이용해서 L2 전처리

      • NA 제거

    # Delete NA Using L1 Data Frame
    dfDataL2 = na.omit(dfData)
    
    # Print Data Frame
    dplyr::tbl_df(dfDataL2)

     

     

    • L2 자료를 출력

      • 옵션 설명

        • sep : 구분자

        • file : 출력 파일명

        • append : 이어쓰기 여부

        • row.names : 행 이름 포함 여부 

        • col.names : 열 이름 포함 여부

    # Write Using L2 Data Frame
    write.table(
       dfDataL2
       , sep = " "
       , file = sOutputFileDirName
       , append = FALSE
       , row.names = FALSE
       , col.names = FALSE
    )

     

     

    • L2 자료를 이용해서 L3 전처리

      • 공간 해상도 일치를 위해서 위/경도에 대한 정수 처리

      • 해양 영역에 대한 마스킹 처리

    # L3 Processing Using L2 Data Frame
    dfDataL3 = dfDataL2 %>%
       dplyr::mutate(
          nRefLat = round(nLat, 0)
          , nRefLon = round(nLon, 0)
          , isLand = metR::MaskLand(dfDataL3$nRefLon, dfDataL3$nRefLat, mask = "world")
       ) %>%
       dplyr::filter(isLand == TRUE)
    
    # Write Using L3 Data Frame
    dplyr::glimpse(dfDataL3)

     

     

    • L3 자료를 이용해서 L4 전처리

      • 2015년 07월 - 2016년 12월에 대한 가강수량의 요약 통계량 (평균, 최대값, 최소값, 개수)

    # Mean For the Period 201507-201612 Using L3 Data Frame
    dfDataL4 = dfDataL3 %>%
       dplyr::group_by(nRefLat, nRefLon) %>%
       dplyr::summarise(
          nMeanTpw = mean(nTpw, na.rm = TRUE)
          , nMaxTpw = max(nTpw, na.rm = TRUE)
          , nMinTpw = min(nTpw, na.rm = TRUE)
          , iNumber = n()
       )
    
    # Write Using L4 Data Frame
    dplyr::tbl_df(dfDataL4)

     

     

    • 가시화를 위한 설정

      • 폰트

      • 컬러바

      • 호주에 대한 Shape 지도

    # Set Value for Visualization
    font = "Palatino Linotype"
    cbMatlab =  colorRamps::matlab.like(11)
    mapAus = read_sf("INPUT/gadm36_AUS_shp/gadm36_AUS_1.shp")

     

    • L4 자료를 이용한 가시화

      • 등고선

      • 호주에 대한 Shape

    # Visualization Using ggplot2
    ggplot() +
       coord_fixed(ratio = 1.1) +
       theme_bw() +
       geom_tile(data = dfDataL4, aes(x = nRefLon, y = nRefLat, fill = nMeanTpw)) +
       geom_text_contour(data = dfDataL4, aes(x = nRefLon, y = nRefLat, z = nMeanTpw), stroke = 0.2) +
       geom_contour(data = dfDataL4, aes(x = nRefLon, y = nRefLat, z = nMeanTpw), color = "black", alpha = 0.3) +
       scale_fill_gradientn(colours = cbMatlab, limits=c(0, 5), breaks = seq(0, 5, 1), na.value = cbMatlab[length(cbMatlab)]) +
       theme(plot.title = element_text(face = "bold", size = 18, color = "black")) +
       theme(axis.title.x = element_text(face = "bold", size = 18, colour = "black")) +
       theme(axis.title.y = element_text(face = "bold", size=18, colour = "black", angle=90)) +
       theme(axis.text.x  = element_text(face = "bold", size=18, colour = "black")) +
       theme(axis.text.y  = element_text(face = "bold", size=18, colour = "black")) +
       metR::scale_x_longitude(expand = c(0, 0), breaks = seq(110, 160, 10), limits = c(110, 160)) +
       metR::scale_y_latitude(expand = c(0, 0), breaks = seq(-40, -10, 5), limits = c(-40, -10)) +
       geom_sf(data = mapAus, color = "black", fill = NA) +
       theme(legend.position = c(1, 1), legend.justification = c(1, 1)) +
       theme(legend.key=element_blank()) +
       theme(legend.text=element_text(size = 14, face = "bold")) +
       theme(legend.title=element_text(face = "bold", size = 14, colour = "black")) +
       labs(x = ""
            , y = ""
            , fill = ""
            , colour = ""
            , title  = ""
       ) +
       theme(legend.background=element_blank()) +
       theme(text=element_text(family = font)) +
       theme(plot.margin = unit(c(0, 8, 0, 0), "mm")) +
       ggsave(filename = paste0('FIG/Tpw.png'), width = 8, height = 10, dpi = 600)

     

    그림. NPP/CERES SSF 1deg 자료를 이용한 가시화

     

    [전체]

     

     참고 문헌

    [논문]

    • 없음

    [보고서]

    • 없음

    [URL]

    • 없음

     

     문의사항

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

    • sangho.lee.1990@gmail.com

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

    • saimang0804@gmail.com
    반응형
    • 네이버 블러그 공유하기
    • 네이버 밴드에 공유하기
    • 페이스북 공유하기
    • 카카오스토리 공유하기