close
프로필 배경
프로필 로고

해솔 | 해처럼 밝고 소나무처럼 바르게

IDL · 2020. 3. 31. fullscreen 넓게보기

[IDL] 아이디엘 대기 환경 자료 처리 및 가시화 : HDF5 형식인 천리안위성 1A호 (COMS/MI) 기상 위성 자료를 이용한 에어로졸 광학두께 가시화

목차

 정보

  • 업무명    : 대기 환경 자료 처리 및 가시화 : HDF5 형식인 천리안위성 1A호 (COMS/MI) 기상 위성 자료를 이용한 에어로졸 광학두께 가시화

  • 작성자    : 이상호

  • 작성일    : 2020-03-31

  • 설   명    :

  • 수정이력 :

 

 내용

[개요]

  • 안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.

  • 대기 환경 정보는 다양한 차원의 관측 자료 및 분석 자료로 구성되어 있기 때문에 이용자가 적절한 분석 및 가시화 기술이 요구됩니다.

  • 따라서 오늘 포스팅에서는 HDF5 형식인 천리안위성 1A호 (COMS/MI) 기상 위성 자료를 이용한 에어로졸 광학두께 가시화를 소개해 드리고자 합니다.

  • 추가로 대기과학 전공자를 위한 IDL를 소개한 링크를 보내드립니다.

 

[IDL] 아이디엘 대기과학 전공자를 위한 IDL (Interactive Data Language) 소개

정보 업무명 : 대기과학 전공자를 위한 IDL (Interactive Data Language) 소개 작성자 : 이상호 작성일 : 2020-03-23 설 명 : 수정이력 : 내용 [개요] 안녕하세요? 웹 개발 및 연구 개발을 담당하고 있는 해솔입니..

shlee1990.tistory.com

 

[특징]

  • 대기 환경 자료를 이해하기 위해서 자료 처리 및 가시화 기술이 요구되며 이 프로그램은 이러한 목적을 달성하기 위해 고안된 소프트웨어

 

[기능]

  • 위도, 경도, 에어로졸 광학두께 자료 읽기

  • 자료 전처리

  • 가시화

 

[활용 자료]

  • 자료명 : coms_mi_le2_aod_cn_201402010000.h5

  • 자료 종류 : 위/경도, 에어로졸 광학두께

  • 확장자 : h5

  • 영역 : 북반구

  • 기간 : 2014년 02월 01일 00:00 UTC

  • 제공처 : 국가기상위성센터

coms_cn_geos_lonlat.z01
10.00MB
coms_cn_geos_lonlat.zip
6.53MB
coms_mi_le2_aod_cn_201402010000.zip
2.17MB

 

[사용법]

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

  • 소스 코드를 실행 (idl -e Visualization_Using_Aerosol_Data_of_HDF5_Format)

  • 가시화 결과를 확인

 

[사용 OS]

  • Windows 10

 

[사용 언어]

  • IDL v8.5

 

 소스 코드

  • 해당 작업을 하기 위한 컴파일 및 실행 코드는 IDL로 작성되었으며 가시화를 위한 라이브러리는 Coyote's Guide to IDL Programming를 이용하였습니다.

  • 소스 코드는 단계별로 수행하며 상세 설명은 다음과 같습니다.

    • 1 단계는 주 프로그램은 작업 경로 설정, HDF5 파일 읽기, 변수 전처리하여 메모리상에 저장하고 가시화를 위한 초기 설정합니다.

    • 2 단계는 plot 매핑에 따라 영상 장면 표출하여 이미지 형식으로 저장합니다. 이 과정에서 포스트 스크립트 (PS) 형식에서 PNG로 변환합니다.

 

[명세]

  • 전역 변수 설정

    • cd를 통해 작업 디렉터리 설정

    • 배열 정보 (dim)

    • 위/경도 (hdf4_file), 에어로졸 광학두께 (hdf5_file) 파일 설정

cd, 'C:\SYSTEM\PROG\IDL\HDF5'

dim = [1934, 1544]

hdf4_file='INPUT/coms_cn_geos_lonlat.hdf'
hdf5_file = "INPUT/coms_mi_le2_aod_cn_201402010000.h5"

 

  • HDF5 파일 읽기 및 전처리

    • h5d_read를 통해 에어로졸 광학두께 정보 읽기

    • 속성 정보 (valid_min, valid_max, fill, offset, scale)를 통해 에어로졸 광학두께 전처리

file_id = h5f_open(hdf5_file)
data_id = h5d_open(file_id, '/Product/Aerosol_Optical_Depth')
data = h5d_read(data_id)
h5f_close, file_id
data = float(data)

valid_min = 0.000000
valid_max = 5.00000
fill = 32767
offset = 0.000000
scale = 0.00488759

loc1 = where( ((data lt valid_min[0]) or (data gt valid_max[0]) or (data eq fill[0])), nloc1, complement=loc2, ncomplement=nloc2)
if (nloc1 gt 0) then data[loc1] = !values.f_nan
if (nloc2 gt 0) then data[loc2] = scale[0] * (data[loc2] - offset[0])
val2D = data

 

  • HDF4 파일 읽기 및 전처리

    • hdf_sd_getdata를 통해 위/경도 정보 읽기
    • 비정상적인 값 (-999.0)에 대한 위/경도 전처리
sd_id=hdf_sd_start(hdf4_file, /read)
sds_index = hdf_sd_nametoindex(sd_id, 'Lat')
sds_id = hdf_sd_select(sd_id, sds_index)
hdf_sd_getdata, sds_id, data
data = float(data)

loc1 = where( (data eq -999.0), nloc1, complement=loc2, ncomplement=nloc2)
if (nloc1 gt 0) then data[loc1] = !values.f_nan
if (nloc2 gt 0) then data[loc2] = data[loc2]
lat2D = data

sd_id=hdf_sd_start(hdf4_file, /read)
sds_index = hdf_sd_nametoindex(sd_id, 'Lon')
sds_id = hdf_sd_select(sd_id, sds_index)
hdf_sd_getdata, sds_id, data
data = float(data)

loc1 = where( (data eq -999.0), nloc1, complement=loc2, ncomplement=nloc2)
if (nloc1 gt 0) then data[loc1] = !values.f_nan
if (nloc2 gt 0) then data[loc2] = data[loc2]
lon2D = data

hdf_sd_endaccess, sds_id
hdf_sd_end, sd_id

 

  • 요약 통계량 확인

    • 최소값 (min), 평균값 (mean), 최대값 (max) 계산

    • help를 통해 자료형 및 배열 정보 확인

print, min(lon2D, /nan), mean(lon2D, /nan), max(lon2D, /nan)
print, min(lat2D, /nan), mean(lat2D, /nan), max(lat2D, /nan)
print, min(val2D, /nan), mean(val2D, /nan), max(val2D, /nan)

help, lon2D, lat2D, val2D

 

etc-image-0

 

  • 가시화

    • 데이터 세트에서 정보는 여전히 숨겨져 있기 때문에 가시화 필요

    • 가시화는 일반적인 정적 탐색 데이터 분석에서 웹 브라우저의 동적 대화식 데이터 시각화에 이르기까지 다양함

 

  • 에어로졸 광학두께를 위한 설정 및 가시화

    • 관심 영역 (ROI)와 중심 위/경도 (cnLon, cnLat) 설정

    • 이미지 저장 파일명 (psName) 및 제목명 (mainName)

    • fnMakePsPlot를 통해 이미지 저장

ROI = [-20, 60, 79.9, 200] 
cnLat = 0
cnLon = 128.2

dtDateYmd = STRMID(hdf5_file, 25, 12)
psName = "Img01_" + dtDateYmd
mainName = "HDF5 Aerosol Optical Depth : " + dtDateYmd

fnMakePsPlot, lon2D, lat2D, val2D, dim, cnLat, cnLon, ROI, psName, mainName, 0, 2, 0.1, 33

 

Img01_201402010000.png

 

  • 사용자 편의성 함수 정의

    • 전달 인자

      • 2차원 위도 (lon2D), 2차원 경도 (lat2D), 2차원 값 (val2D), 배열 정보 (dim), 중심 위도 (cnLat), 중심 경도 (cnLon), 관심 영역 (ROI), 이미지 저장 파일명 (psName), 제목명 (mainName), 최소값 (zmin), 최대값 (zmax), 등고선 간격 (interval), 컬러 테이블 (color_table)

    • 기능

      • plot를 통해 위/경도에 따른 값 매핑

      • 포스트 스크립트 (PS)를 PNG로 변환

;====================================================================
; Subroutine : Function Make Postscript to Png Graph
;====================================================================
pro fnMakePsPlot, lon2D, lat2D, val2D, dim, cnLat, cnLon, ROI, psName, mainName, zmin, zmax, interval, color_table

    set_plot, "ps"

    device, filename = psName + ".ps", decomposed = 0, bits = 8, /color, xsize = 16, ysize = 11, /inches, font_size = 13, /Helvetica

    !p.font = 0 & !p.charsize = 2.0 & !p.charthick = 1.6  & !p.multi = [0,1,1] & !p.background = 255

    xvert = [0.1, 0.1, -0.1, -0.1 ,0.1]
    yvert = [-0.1, 0.1, 0.1, -0.1, -0.1]

    usersym, xvert, yvert, /fill

    start_color = 0 & end_color = 255 & colorn = end_color - start_color + 1

    cgloadct, color_table

    latmin = ROI(0) & latmax= ROI(1) & lonmin = ROI(2) & lonmax= ROI(3)

    cgMAP_SET, cnLat, cnLon, /CYLINDRICAL $
      , limit = [latmin, lonmin, latmax, lonmax] $
      , position = [0.05, 0.0, 0.85, 0.90], /noborder, /ISOTROPIC

    for i = 0L, dim[0] - 1 do begin
      for j = 0L, dim[1] - 1 do begin

        nLon = lon2D[i, j]
        nLat = lat2D[i, j]
        nVal = val2D[i, j]

        if (lonmin gt nLon or nLon gt lonmax or finite(nLon, /nan) eq 1) then continue
        if (latmin gt nLat or nLat gt latmax or finite(nLat, /nan) eq 1) then continue
        if (finite(nVal, /nan) eq 1) then continue

        plots, nLon, nLat, psym = 8, symsize = 2, color = BYTSCL(nVal, zmin, zmax)

      endfor
    endfor

    cgColorbar, NColors = 255, BOTTOM = 1 $
      , Divisions = 5, COLOR = "black" $
      , Position = [0.1, 0.90, 0.8, 0.94] $
      , Range = [zmin, zmax] $
      , CHARSIZE = 1.5,CHARthick = 2 $
      , VERTICAL = 1, RIGHT = 1

  ;  cgLoadct, 2

  ;  contour, val2D, lon, lat, /overplot, C_THICK = 2 $
  ;    , LEVELS = zmin + findgen(10) * interval $
  ;    , C_LABELS = zmin + findgen(10) * interval $
  ;    , C_CHARSIZE = 1.5, C_CHARTHICK = 2 $
  ;    , C_COLORS = BYTSCL(zmin + findgen(10) * interval, zmin, zmax)

    cgMAP_CONTINENTS, /coast, /countries, COLOR = "black"

    lats = [-90:90:15]
    lons = [-180:360:20]

    lats_names = strarr(n_elements(lats))
    lons_names = strarr(n_elements(lons))

    for i = 0L, n_elements(lats) - 1 do begin
      if (lats[i] gt 0) then begin
        lats_names[i] = textoidl(string(lats[i]) + "\circN")
      endif else if (lats[i] eq 0) then begin
        lats_names[i] = textoidl(string(lats[i]) + "\circ")
      endif else begin
        lats_names[i] = textoidl(string(-lats[i]) + "\circS")
      endelse
    endfor

    for i = 0L,n_elements(lons) - 1  do begin
      if (lons[i] gt 0) then begin
        lons_names[i] = textoidl(string(lons[i]) + "\circE")
      endif else if (lons[i] eq 0) then begin
        lons_names[i] = textoidl(string(lons[i]) + "\circ")
      endif else begin
        lons_names[i] = textoidl(string(-lons[i]) + "\circW")
      endelse
    endfor

    cgMap_grid, color = "black", charsize = 1.8, thick = 2.5, lats = lats,  latnames = lats_names, lons = lons, lonnames = lons_names $
               , linestyle = 1, bthick = 300, /box_axes; , /no_grid

    cgText, 0.45, 0.95, mainName, /normal, charsize = 2, CHARTHICK = 2, color = "black", alignment = 0.5

    device, /close_file

    com = "convert -flatten -background white " + psName + ".ps" + " " + "./FIG/" + file_basename(psName, ".ps") + ".png"
    spawn, /hide, com

    file_delete, ps_name + ".ps"

    return

end

 

[전체] 

pro Visualization_Using_Aerosol_Data_of_HDF5_Format
cd, 'C:\SYSTEM\PROG\IDL\HDF5'
dim = [1934, 1544]
hdf4_file='INPUT/coms_cn_geos_lonlat.hdf'
hdf5_file = "INPUT/coms_mi_le2_aod_cn_201402010000.h5"
file_id = h5f_open(hdf5_file)
data_id = h5d_open(file_id, '/Product/Aerosol_Optical_Depth')
data = h5d_read(data_id)
h5f_close, file_id
data = float(data)
valid_min = 0.000000
valid_max = 5.00000
fill = 32767
offset = 0.000000
scale = 0.00488759
loc1 = where( ((data lt valid_min[0]) or (data gt valid_max[0]) or (data eq fill[0])), nloc1, complement=loc2, ncomplement=nloc2)
if (nloc1 gt 0) then data[loc1] = !values.f_nan
if (nloc2 gt 0) then data[loc2] = scale[0] * (data[loc2] - offset[0])
val2D = data
sd_id=hdf_sd_start(hdf4_file, /read)
sds_index = hdf_sd_nametoindex(sd_id, 'Lat')
sds_id = hdf_sd_select(sd_id, sds_index)
hdf_sd_getdata, sds_id, data
data = float(data)
loc1 = where( (data eq -999.0), nloc1, complement=loc2, ncomplement=nloc2)
if (nloc1 gt 0) then data[loc1] = !values.f_nan
if (nloc2 gt 0) then data[loc2] = data[loc2]
lat2D = data
sd_id=hdf_sd_start(hdf4_file, /read)
sds_index = hdf_sd_nametoindex(sd_id, 'Lon')
sds_id = hdf_sd_select(sd_id, sds_index)
hdf_sd_getdata, sds_id, data
data = float(data)
loc1 = where( (data eq -999.0), nloc1, complement=loc2, ncomplement=nloc2)
if (nloc1 gt 0) then data[loc1] = !values.f_nan
if (nloc2 gt 0) then data[loc2] = data[loc2]
lon2D = data
hdf_sd_endaccess, sds_id
hdf_sd_end, sd_id
dims = size(var, /dimensions)
print, min(lon2D, /nan), mean(lon2D, /nan), max(lon2D, /nan)
print, min(lat2D, /nan), mean(lat2D, /nan), max(lat2D, /nan)
print, min(val2D, /nan), mean(val2D, /nan), max(val2D, /nan)
help, lon2D, lat2D, val2D
ROI = [-20, 60, 79.9, 200]
cnLat = 0
cnLon = 128.2
dtDateYmd = STRMID(hdf5_file, 25, 12)
psName = "Img01_" + dtDateYmd
mainName = "HDF5 Aerosol Optical Depth : " + dtDateYmd
fnMakePsPlot, lon2D, lat2D, val2D, dim, cnLat, cnLon, ROI, psName, mainName, 0, 2, 0.1, 33
end
;====================================================================
; Subroutine : Function Make Postscript to Png Graph
;====================================================================
pro fnMakePsPlot, lon2D, lat2D, val2D, dim, cnLat, cnLon, ROI, psName, mainName, zmin, zmax, interval, color_table
set_plot, "ps"
device, filename = psName + ".ps", decomposed = 0, bits = 8, /color, xsize = 16, ysize = 11, /inches, font_size = 13, /Helvetica
!p.font = 0 & !p.charsize = 2.0 & !p.charthick = 1.6 & !p.multi = [0,1,1] & !p.background = 255
xvert = [0.1, 0.1, -0.1, -0.1 ,0.1]
yvert = [-0.1, 0.1, 0.1, -0.1, -0.1]
usersym, xvert, yvert, /fill
start_color = 0 & end_color = 255 & colorn = end_color - start_color + 1
cgloadct, color_table
latmin = ROI(0) & latmax= ROI(1) & lonmin = ROI(2) & lonmax= ROI(3)
cgMAP_SET, cnLat, cnLon, /CYLINDRICAL $
, limit = [latmin, lonmin, latmax, lonmax] $
, position = [0.05, 0.0, 0.85, 0.90], /noborder, /ISOTROPIC
for i = 0L, dim[0] - 1 do begin
for j = 0L, dim[1] - 1 do begin
nLon = lon2D[i, j]
nLat = lat2D[i, j]
nVal = val2D[i, j]
if (lonmin gt nLon or nLon gt lonmax or finite(nLon, /nan) eq 1) then continue
if (latmin gt nLat or nLat gt latmax or finite(nLat, /nan) eq 1) then continue
if (finite(nVal, /nan) eq 1) then continue
plots, nLon, nLat, psym = 8, symsize = 2, color = BYTSCL(nVal, zmin, zmax)
endfor
endfor
cgColorbar, NColors = 255, BOTTOM = 1 $
, Divisions = 5, COLOR = "black" $
, Position = [0.1, 0.90, 0.8, 0.94] $
, Range = [zmin, zmax] $
, CHARSIZE = 1.5,CHARthick = 2 $
, VERTICAL = 1, RIGHT = 1
; cgLoadct, 2
; contour, val2D, lon, lat, /overplot, C_THICK = 2 $
; , LEVELS = zmin + findgen(10) * interval $
; , C_LABELS = zmin + findgen(10) * interval $
; , C_CHARSIZE = 1.5, C_CHARTHICK = 2 $
; , C_COLORS = BYTSCL(zmin + findgen(10) * interval, zmin, zmax)
cgMAP_CONTINENTS, /coast, /countries, COLOR = "black"
lats = [-90:90:15]
lons = [-180:360:20]
lats_names = strarr(n_elements(lats))
lons_names = strarr(n_elements(lons))
for i = 0L, n_elements(lats) - 1 do begin
if (lats[i] gt 0) then begin
lats_names[i] = textoidl(string(lats[i]) + "\circN")
endif else if (lats[i] eq 0) then begin
lats_names[i] = textoidl(string(lats[i]) + "\circ")
endif else begin
lats_names[i] = textoidl(string(-lats[i]) + "\circS")
endelse
endfor
for i = 0L,n_elements(lons) - 1 do begin
if (lons[i] gt 0) then begin
lons_names[i] = textoidl(string(lons[i]) + "\circE")
endif else if (lons[i] eq 0) then begin
lons_names[i] = textoidl(string(lons[i]) + "\circ")
endif else begin
lons_names[i] = textoidl(string(-lons[i]) + "\circW")
endelse
endfor
cgMap_grid, color = "black", charsize = 1.8, thick = 2.5, lats = lats, latnames = lats_names, lons = lons, lonnames = lons_names $
, linestyle = 1, bthick = 300, /box_axes; , /no_grid
cgText, 0.45, 0.95, mainName, /normal, charsize = 2, CHARTHICK = 2, color = "black", alignment = 0.5
device, /close_file
com = "convert -flatten -background white " + psName + ".ps" + " " + "./FIG/" + file_basename(psName, ".ps") + ".png"
spawn, /hide, com
file_delete, ps_name + ".ps"
return
end
view raw Visualization_Using_Ozone_Data_of_HDF5_Format.pro hosted with ❤ by GitHub

 

 참고 문헌

[논문]

  • 없음

[보고서]

  • 없음

[URL]

  • 없음

 

 문의사항

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

  • sangho.lee.1990@gmail.com

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

  • saimang0804@gmail.com
저작자표시 비영리 변경금지

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

[IDL] 아이디엘 대기 환경 자료 처리 및 가시화 : 태양 위치 정보를 이용하여 태양 천정각 (Solar Zenith Angle) 및 태양 방위각 (Solar Azimuth Angle) 가시화  (0) 2020.04.01
[IDL] 아이디엘 대기 환경 자료 처리 및 가시화 : Grib 형식인 ECMWF Interim 기상 자료를 이용한 가시화  (0) 2020.03.31
[IDL] 아이디엘 대기 환경 자료 처리 및 가시화 : HDF 형식인 전지구 오존량 자료를 이용한 가시화  (0) 2020.03.31
[IDL] 아이디엘 대기 환경 자료 처리 및 가시화 : NetCDF 형식인 대기 온도 자료를 이용한 가시화  (0) 2020.03.31
[IDL] 아이디엘 대기 환경 자료 처리 및 가시화 : 아스키 (ASCII) 형식인 전지구 오존량 자료를 이용한 가시화  (0) 2020.03.31
프로그래밍 언어/IDL 관련 글
더 보기

티스토리툴바