반응형
정보
-
업무명 : 아이디엘 히마와리 8호 (Himawari-8/AHI) 기상 위성 자료를 이용한 RGB 합성영상 가시화
-
작성자 : 이상호
-
작성일 : 2019-08-30
-
설 명 :
-
수정이력 :
-
2020-04-09 : 소스 코드 명세 추가
-
내용
[개요]
-
안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
-
RGB 합성영상은 빛의 삼원색인 빨강, 초록, 파랑의 색을 조합하여 색을 재현하는 기술입니다.
-
주로 위성분야에서는 위성에서 관측한 채널의 반사율 또는 밝기온도의 차이에 색을 조합하여 영상을 합성합니다.
-
이는 기존의 흑백영상에 비해 다양한 색깔로 대기의 상태를 표현하기 때문에 효과적인 위성 분석뿐만 아니라 구름의 변화 및 물리적인 특성 변화를 쉽게 파악할 수 있습니다.
-
오늘 포스팅에서는 히마와리 8호 (Himawari-8/AHI) 기상 위성 자료를 이용한 RGB 합성영상 가시화를 소개해 드리고자 합니다.
[특징]
-
시시각각 변화하는 기상 변화를 시각화하기 위해서 RGB 합성 기술이 요구되며 이 프로그램은 이러한 목적을 달성하기 위해 고안된 소프트웨어
[기능]
-
히마와리 8호 AHI 센서 (Himawari-8/AHI)의 RGB (Red: CH03, Green: CH02 ,Blue: CH01) 채널을 입력자료를 이용하여 단계별로 RGB 합성 및 RGB 천연색 가시화
-
채널별 DN (Digital Number)값은 선형변환 계수를 통해 물리적인 복사량으로 환산
-
복사량은 다음과 같은 식1에 의거하여 대기상단에서의 반사율로 변환하여 RGB 합성
-
RGB 천연색 표출하기 위해 RGB 합성에서 레일라이 산란 보정
[사용법]
-
입력 자료를 해당 디렉터리에 저장
-
소스 코드를 실행 (idl -e RGB_Composite_Image_Using_Himawari8_AHI_Data.pro)
-
단계별로 가시화 결과 (RGB 합성영상, 레일라이 산란 보정, RGB 천연색)확인
[사용 OS]
-
Linux
-
Windows 10
[사용 언어]
-
IDL v8.5
소스 코드
-
해당 작업을 하기 위한 컴파일 및 실행 코드는 IDL로 작성되었으며 가시화를 위한 라이브러리는 Coyote's Guide to IDL Programming를 이용하였습니다.
-
소스 코드는 단계별로 수행하며 상세 설명은 다음과 같습니다.
-
1 단계는 주 프로그램은 작업 경로 설정 및 전역 변수를 설정하여 메모리상에 저장하고 가시화를 위한 초기 설정합니다.
-
2 단계는 RGB 합성 영상 장면 표출하여 이미지 형식으로 저장합니다. 이 과정에서 포스트 스크립트 (PS) 형식에서 PNG로 변환합니다.
-
[명세]
-
작업 환경 설정
cd, 'D:\할 일\Program\각종 자료\Aerosol mapping+RGB\modis hdf_ps_RGB x\Ozone'
-
전역 변수 설정
-
배열 정보 정의 (dims)
-
입력 자료 설정 (datpath, cofpath, calc, llpath)
-
날짜 및 시간 (time)
-
입력 자료 관련 정보로서 서두 (head), 해상도 (res), 영역 (area)
-
dims = [1375L, 1375L] ; Full disk_8km
datpath = 'INPUT'
cofpath = 'INPUT/'
calc = cofpath+path_sep()+'calibration_coeff.txt' ; Coeff
llpath = '.'
time = '201508010000' ; case day
head = 'himawari8_ahi_le1b_'
res = '8km'
area = '_f_'
-
히마와리 8호 이진 자료에 대한 입력 자료 설정
ch01_file = datpath+path_sep()+head+'b01'+'_'+res+area+time+'.bin'
ch02_file = datpath+path_sep()+head+'b02'+'_'+res+area+time+'.bin'
ch03_file = datpath+path_sep()+head+'b03'+'_'+res+area+time+'.bin'
ch04_file = datpath+path_sep()+head+'b04'+'_'+res+area+time+'.bin'
-
채널별 디지털값 (count01, count02, count03, count04) , 위도 (lat), 경도 (lon), 태양 천정각 (SZA)에 대한 배열 설정
count01 = intarr(dims[0], dims[1]) & count02 = intarr(dims[0], dims[1]) ; Count
count03 = intarr(dims[0], dims[1]) & count04 = intarr(dims[0], dims[1])
lat = fltarr(dims[0], dims[1]) ; latitude
lon = fltarr(dims[0], dims[1]) ; longitude
SZA = fltarr(dims[0], dims[1]) ; solor zenith angle
-
히마와리 8호 이진 자료 (count01, count02, count03, count04) 및 위/경도 (lat, lon) 읽기
;******************************************************************
; Read Himawari data
;******************************************************************
openr, lun, ch01_file, /get_lun & readu, lun, count01 & free_lun, lun
openr, lun, ch02_file, /get_lun & readu, lun, count02 & free_lun, lun
openr, lun, ch03_file, /get_lun & readu, lun, count03 & free_lun, lun
openr, lun, ch04_file, /get_lun & readu, lun, count04 & free_lun, lun
openr, lun, 'INPUT/lat.dat', /get_lun & readf, lun, lat & free_lun, lun
openr, lun, 'INPUT/lon.dat', /get_lun & readf, lun, lon & free_lun, lun
-
채널별 디지털값 (Digital Number) 및 선형변환 계수 (gain_vis, off_vis)를 통해 물리적인 복사량 (복사휘도)으로 환산
-
이상치 및 최대/최소값에 대한 NaN 처리
-
;******************************************************************
; Himawari-8 AHI DN ==> Narrowband Radiance, IR
;******************************************************************
READCOL, calc, format='f,f,f,f', wave_vis, gain_vis, off_vis, rad_to_vis
rad01 = (gain_vis[0]*(count01)) + off_vis[0]
rad02 = (gain_vis[1]*(count02)) + off_vis[1]
rad03 = (gain_vis[2]*(count03)) + off_vis[2]
rad04 = (gain_vis[3]*(count04)) + off_vis[3]
bad01 = where(count01 LT 0 OR count01 GE 2048) & rad01[bad01] = !values.f_nan
bad02 = where(count02 LT 0 OR count02 GE 2048) & rad02[bad02] = !values.f_nan
bad03 = where(count03 LT 0 OR count03 GE 2048) & rad03[bad03] = !values.f_nan
bad04 = where(count04 LT 0 OR count04 GE 2048) & rad04[bad04] = !values.f_nan
loc1 = where(lat eq -999.0, nloc1) & if (nloc1 gt 0) then lat[loc1] = !values.f_nan
loc2 = where(lon eq -999.0, nloc2) & if (nloc2 gt 0) then lon[loc2] = !values.f_nan
-
태양 천정각 (Solar Zenith Angle, SZA), 태양 방위각 (Solar Azimuth Angle, SAA), 위성 천정각 (Viewing Zenith Angle, VZA), 위성 방위각 (Viewing Azimuth Angle, VAA), 산란각 (Scattering Angle, SA), 태양 반짝임 (Glint Angle, GA) 계산
;******************************************************************
; SZA, SAA, VZA, VAA, SA, RAA
;******************************************************************
year = long(strmid(time, 0, 4))
month = long(strmid(time, 4, 2))
day = long(strmid(time, 6, 2))
mon = [31,28,31,30,31,30,31,31,30,31,30,31]
if ((((year MOD 4) eq 0) and ((year mod 100) ne 0)) or (year mod 400) eq 0) then mon[1]=29
jd = 0
for i = 0, 11 do begin
for j = 0, mon[i]-1 do begin
jd = jd+1
if ((i+1 eq month) and (j+1 eq day)) then jday=jd
endfor
endfor
; input parameter
time2 = long(strmid(time, 8, 6)+'00') ; HH MM SS
clat = -0.031184d & clon = 140.204602d
; SZA = solar zenith angle
HH = fix( FLOAT(TIME2) / 10000.0)
MM = fix((FLOAT(TIME2)/100) - FLOAT(HH) * 100)
SS = TIME2 - HH * 10000L - MM * 100L
tu = FLOAT(HH) + FLOAT(MM)/60.0 + FLOAT(SS)/3600.0
pi = !pi
rad = !pi / 180.
glat = clat
glon = clon
HH = fix( FLOAT(TIME2) / 10000.0)
MM = fix((FLOAT(TIME2)/100) - FLOAT(HH) * 100)
SS = TIME2 - HH * 10000l - MM * 100l
DTIME = FLOAT(HH) + FLOAT(MM)/60.0 + FLOAT(0)/3600.0
H_ANG = -1.0 * !PI * (DTIME-12.0) / 12.0
nc = 0
nd = 0
DECLIN = -23.45 * COS(FLOAT(jday)*2.*!PI/365. + 10.*2.*!PI/365.)
na = where(lat ge -90. and lat le 90. and lat gt declin,na,$
complement=b,ncomplement=nb)
gd = where(lat ge -90. and lat le 90.) ;,complement=ng)
g=gd
betaa=lat
azimuth=lat
betaa(*)=0.
mu_not = lat
mu_not(*,*)=0.
gd = where(lat ge -90. and lat le 90.) ;,complement=ng)
ng = where(lat lt -90. or lat gt 90.) ;,complement=ng)
MU_NOT(gd) = (SIN(LAT(gd)*RAD)*SIN(DECLIN*RAD) + $
COS(LAT(gd)*RAD)*COS(DECLIN*RAD)*COS(H_ANG-LON(gd)*RAD))
solzen = acos(mu_not < 1.) / rad
mu_not = 0
nc = 0
nd = 0
a = 0
; VZA (viewing zenith angle)
betaa=lat
theta=lat
betaa(*)=0.
betaa(g) = ACOS(COS(LAT(g)*RAD) * COS((GLON - LON(g))*RAD))
theta(g) = asin((42164.*sin(betaa(g)) / $
sqrt(1.8084e9 - 5.3725e8*cos(betaa(g)))) < 1) / rad
; SAA (solar azimuth angle)
a = where(lat ge -90. and lat le 90. and lat gt declin,na,$
complement=b,ncomplement=nb)
g=gd
gd = 0
az_not = lat
az_not(*,*) = 0.
if (na gt 0) then az_not(a) = ACOS( COS((LAT(a)-DECLIN)*RAD) * COS(H_ANG-LON(a)*RAD))
if (nb gt 0) then az_not(b) = ACOS( COS((LAT(b)-DECLIN)*RAD) * $
COS(LON(b)*RAD-H_ANG))
if (na gt 0) then AZ_NOT(a) = SIN(H_ANG-LON(a)*RAD) / SIN(az_not(a))
if (nb gt 0) then AZ_NOT(b) = SIN(LON(b)*RAD-H_ANG) / SIN(az_not(b))
az_not(g) = (az_not(g) > (-1.)) < 1.
AZ_NOT(g) = ASIN( AZ_NOT(g) )
if (na gt 0) then AZ_NOT(a) = !PI - AZ_NOT(a)
if (nb gt 0) then AZ_NOT(b) = !PI + AZ_NOT(b)
AZ_NOT(g) = AZ_NOT(g) /RAD
nc = 0
nd = 0
a = 0
if (nb gt 0) then c = where(az_not(b) lt 180.,nc,complement=d,ncomplement=nd)
if (nc gt 0) then az_not(b(c)) = 180. - az_not(b(c))
if (nd gt 0) then az_not(b(d)) = 540. - az_not(b(d))
; RAA (relative azimuth angle)
; VAZ (viewing azimuth angle)
; SA (scattering angle)
betaa=lat
betaa(*)=0.
betaa(g) = ACOS(COS(LAT(g)*RAD) * COS((GLON - LON(g))*RAD))
a = where(lat gt (-90) and lat le 90. and lat gt glat,na,$
complement=b,ncomplement=nb)
azimuth=lat
azimuth(*,*)=0.
if (na gt 0) then AZIMUTH(a) = SIN((GLON-LON(a))*RAD) / SIN(BETAa(a))
if (nb gt 0) then AZIMUTH(b) = SIN((LON(b)-GLON)*RAD) / SIN(BETAa(b))
azimuth(g) = (azimuth(g) > (-1.)) < 1.
azimuth(g) = asin(azimuth(g))
if (na gt 0) then azimuth(a) = !pi - azimuth(a)
a = 0
if (b(0) ge 0) then azimuth(b) = !pi + azimuth(b)
azimuth(g) = azimuth(g) / rad
if (nb gt 0) then c = where(azimuth(b) lt 180.,nc)
if (nb gt 0) then d = where(azimuth(b) ge 180.,nd)
if (nb gt 0 and nc gt 0) then azimuth(b(c)) = 180. - azimuth(b(c))
if (nb gt 0 and nd gt 0) then azimuth(b(d)) = 540. - azimuth(b(d))
imgsza=solzen
imgvza=theta
imgazi=(az_not-azimuth)
idx1=where (imgazi lt 0.,n1)
if (n1 gt 0 ) then imgazi(idx1)=imgazi(idx1)+2.*!pi
idx2=where (imgazi gt (2.*!pi),n2)
if (n2 gt 0 ) then imgazi(idx2)=imgazi(idx2)-2.*!pi
rad = !dtor
cosga =-cos(imgsza*rad)*cos(imgvza*rad) -sin(imgsza*rad)*sin(imgvza*rad)*cos(imgazi*rad)
cosga = (cosga > (-1.)) < 1.
sca_ang=acos(cosga) / rad
;====================================
; LSH
;=====================================
SZA = solzen
SAA = az_not
VZA = theta
VAA = azimuth
RAA = fltarr(dims[0], dims[1])
SA = fltarr(dims[0], dims[1])
GA = fltarr(dims[0], dims[1])
RAA = abs(SAA-VAA-180.0)
loc1 = where(RAA gt 360.0, nloc1) & if(nloc1 gt 0) then RAA[loc1]=RAA[loc1] mod 360.0
loc2 = where(RAA gt 180.0, nloc2) & if(nloc2 gt 0) then RAA[loc2]=360.0 - RAA[loc2]
; SA (Scatteriing Angle)
SA = -cos(SZA*!dtor)*cos(VZA*!dtor) $
+sin(SZA*!dtor)*sin(VZA*!dtor) $
*cos(RAA*!dtor)
SA = acos(SA)*!radeg
; GA (Glint Angle)
GA = cos(SZA*!dtor)*cos(VZA*!dtor) $
+sin(SZA*!dtor)*sin(VZA*!dtor) $
*cos(RAA*!dtor)
GA = acos(GA)*!radeg
-
이상치 및 최대/최소값에 대한 NaN 처리
;******************************************************************
; NaN
;******************************************************************
loc1 = where(lat eq -999.0, nloc1) & if (nloc1 gt 0) then lat[loc1] = !values.f_nan
loc1 = where(lon eq -999.0, nloc1) & if (nloc1 gt 0) then lon[loc1] = !values.f_nan
loc1 = where(sza lt 0 or sza gt 80, nloc1) & if (nloc1 gt 0) then sza[loc1] = !values.f_nan
loc1 = where(vza lt 0 or vza gt 80, nloc1) & if (nloc1 gt 0) then vza[loc1] = !values.f_nan
-
채널별 협대역 복사휘도 (rad01-rad02)를 이용하여 협대역 반사율 (ref01-ref04)로 변환
-
이상치 및 최대/최소값에 대한 NaN 처리
-
;******************************************************************
; Narrowband Radiance ==> Narrowband Reflectance
;******************************************************************
var = !pi/(cos(SZA*!dtor))
ref01 = ((var*rad01)/2015.3606)
ref02 = ((var*rad02)/1891.1653)
ref03 = ((var*rad03)/1631.5726)
ref04 = ((var*rad04)/971.8778)
loc1 = where(ref01 lt 0 or ref01 gt 1, nloc1) & if (nloc1 gt 0) then ref01[loc1] = !values.f_nan
loc1 = where(ref02 lt 0 or ref02 gt 1, nloc1) & if (nloc1 gt 0) then ref02[loc1] = !values.f_nan
loc1 = where(ref03 lt 0 or ref03 gt 1, nloc1) & if (nloc1 gt 0) then ref03[loc1] = !values.f_nan
-
RGB 합성 영상을 위한 설정 및 가시화
-
24 비트 (bits=24) 설정
-
plots를 통해 위/경도 (lon, lat)에 따른 값 (val) 매핑
-
포스트 스크립트 (PS)를 PNG로 변환
-
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> True RGB
;******************************************************************
ps_name = 'Himawari_True_RGB_'+time
set_plot, 'ps'
device, filename=ps_name+'.ps', decomposed=1, bits=24, /color, xsize=12, ysize=13, /inches, $
/schoolbook, /bold
!p.font=0 & !p.charsize=2.0 & !p.charthick=1.6 & !p.multi=[0,1,1]
scaledBand = scalemodis(ref03, ref02, ref01)
red = fltarr(dims[0], dims[1])
green = fltarr(dims[0], dims[1])
blue = fltarr(dims[0], dims[1])
red = reform(scaledBand[*,*,0])
green = reform(scaledBand[*,*,1])
blue = reform(scaledBand[*,*,2])
loc1 = where( (red eq 0) and (green eq 0) and (blue eq 0), nloc1 )
if (nloc1 gt 0) then begin
red[loc1] = 255
green[loc1] = 255
blue[loc1] = 255
endif
region = [-90, 90, -180, 180] ; All earth
; region = [ 32, 44, 122, 134] ; Hanbando
; region = [ 10, 100, 50, 132] ; East asia
; region = [ 32, 120, 42, 132] ; GOCI (1)
; region = [ 25, 110, 50, 150] ; GOCI (2)
latmin=region[0] & latmax=region[1]
lonmin=region[2] & lonmax=region[3]
clat = -0.007606d & clon = 140.714793d
cgmap_set, /satellite, sat_p=[42159.934825,0,0], clat, clon, charsize=1.0, limit=[-90,-180,90,180], $
position=[0.10,0.15,0.9,0.95], /isotropic, color='white', xmargin=0, ymargin=0, /noerase
usersym, [0.1,0.1,-0.1,-0.1,0.1], [-0.1,0.1,0.1,-0.1,-0.1], /fill
for j=0, dims[1]-1 do begin
for i=0, dims[0]-1 do begin
val = red[i,j] + (green[i,j]*(2L^8)) + (blue[i,j]*(2L^16))
plots, lon[i,j], lat[i,j], psym=8, symsize=1.0, color=val
endfor
endfor
cgmap_gshhs, 'gshhg'+path_sep()+'gshhs_i.b', color='black', /outline, level=4, thick=1.0
cgmap_continents, /continents, /hires, color='black', mlinethick=0.001, thick=1.0, limit=[-90,-180,90,180]
lats = [-90, -60, -30, 0, 30, 60, 90]
; lats_names = ['',textoidl('60\circS', font=0), textoidl('30\circS', font=0), ' ', textoidl('30\circN', font=0),textoidl('60\circN', font=0), '']
lats_names = ['', '', '', '', '', '', '']
lons = [10, 80, 110, 140, 170, -160, -10]
; lons_names = ['',textoidl('80\circE'), textoidl('110\circE'), textoidl('140\circE'), textoidl('170\circE'), textoidl('160\circW'),'']
lons_names = ['', '', '', '', '', '', '']
cgmap_grid, color='black', /label, lats=lats, latlabel=170, lonlabel=-5, latnames=lats_names, lons=lons, lonnames=lons_names, clip_text=0, linestyle=1, noclip=0, /horizon, charsize=1.7
cgtext, 0.5, 0.965, ps_name, charsize=2, charthick=1.6, alignment=0.5, color='black', font=0, /normal
device, /close_file
com = 'convert -flatten -background white '+ps_name+'.ps'+' OUTPUT/'+file_basename(ps_name, '.ps')+'.png'
spawn, com, /hide
file_delete, ps_name+'.ps'
-
히마와리 8호 (Himawari-8/AHI) 자료를 이용한 RGB 합성영상
-
레일라이 산란 보정을 위한 설정 및 가시화
-
quick_Ray_crr_single를 통해 채널별 반사율을 보정
-
plots를 통해 위/경도에 따른 값 매핑
-
포스트 스크립트 (PS)를 PNG로 변환
-
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> RGB Rayleigh Correction
;******************************************************************
ps_name = 'Himawari_RGB_Rayleigh_Correction_'+time
wvldesnm=470.63 ; input wvl in nm
quick_Ray_crr_single, VZA, SZA, RAA, SA, wvldesnm, ref01, ref01_ray
wvldesnm=510.00 ; input wvl in nm
quick_Ray_crr_single, VZA, SZA, RAA, SA, wvldesnm, ref02, ref02_ray
wvldesnm=649.14 ; input wvl in nm
quick_Ray_crr_single, VZA, SZA, RAA, SA, wvldesnm, ref03, ref03_ray
wvldesnm=856.70 ; input wvl in nm
quick_Ray_crr_single, VZA, SZA, RAA, SA, wvldesnm, ref04, ref04_ray
loc1 = where(ref01_ray lt 0 or ref01_ray gt 1, nloc1) & if (nloc1 gt 0) then ref01_ray[loc1] = !values.f_nan
loc1 = where(ref02_ray lt 0 or ref02_ray gt 1, nloc1) & if (nloc1 gt 0) then ref02_ray[loc1] = !values.f_nan
loc1 = where(ref03_ray lt 0 or ref03_ray gt 1, nloc1) & if (nloc1 gt 0) then ref03_ray[loc1] = !values.f_nan
loc1 = where(ref04_ray lt 0 or ref04_ray gt 1, nloc1) & if (nloc1 gt 0) then ref04_ray[loc1] = !values.f_nan
;
set_plot, 'ps'
device, filename=ps_name+'.ps', decomposed=1, bits=24, /color, xsize=12, ysize=13, /inches, $
/schoolbook, /bold
!p.font=0 & !p.charsize=2.0 & !p.charthick=1.6 & !p.multi=[0,1,1]
scaledBand = scalemodis(ref03_ray, ref02_ray, ref01_ray)
red = fltarr(dims[0], dims[1])
green = fltarr(dims[0], dims[1])
blue = fltarr(dims[0], dims[1])
red = reform(scaledBand[*,*,0])
green = reform(scaledBand[*,*,1])
blue = reform(scaledBand[*,*,2])
loc1 = where( (red eq 0) and (green eq 0) and (blue eq 0), nloc1 )
if (nloc1 gt 0) then begin
red[loc1] = 255
green[loc1] = 255
blue[loc1] = 255
endif
region = [-90, 90, -180, 180] ; All earth
; region = [ 32, 44, 122, 134] ; Hanbando
; region = [ 10, 100, 50, 132] ; East asia
; region = [ 32, 120, 42, 132] ; GOCI (1)
; region = [ 25, 110, 50, 150] ; GOCI (2)
latmin=region[0] & latmax=region[1]
lonmin=region[2] & lonmax=region[3]
clat = -0.007606d & clon = 140.714793d
cgmap_set, /satellite, sat_p=[42159.934825,0,0], clat, clon, charsize=1.0, limit=[-90,-180,90,180], $
position=[0.10,0.15,0.9,0.95], /isotropic, color='white'
usersym, [0.1,0.1,-0.1,-0.1,0.1], [-0.1,0.1,0.1,-0.1,-0.1], /fill
for j=0, dims[1]-1 do begin
for i=0, dims[0]-1 do begin
val = red[i,j] + (green[i,j]*(2L^8)) + (blue[i,j]*(2L^16))
plots, lon[i,j], lat[i,j], psym=8, symsize=1.0, color=val
endfor
endfor
cgmap_gshhs, 'gshhg'+path_sep()+'gshhs_i.b', color='black', /outline, level=4, thick=1.0
cgmap_continents, /continents, /hires, color='black', mlinethick=0.001, thick=1.0, limit=[-90,-180,90,180]
lats = [-90, -60, -30, 0, 30, 60, 90]
; lats_names = ['',textoidl('60\circS', font=0), textoidl('30\circS', font=0), ' ', textoidl('30\circN', font=0),textoidl('60\circN', font=0), '']
lats_names = ['', '', '', '', '', '', '']
lons = [10, 80, 110, 140, 170, -160, -10]
; lons_names = ['',textoidl('80\circE'), textoidl('110\circE'), textoidl('140\circE'), textoidl('170\circE'), textoidl('160\circW'),'']
lons_names = ['', '', '', '', '', '', '']
cgmap_grid, color='black', /label, lats=lats, latlabel=170, lonlabel=-5, latnames=lats_names, lons=lons, lonnames=lons_names, clip_text=0, linestyle=1, noclip=0, /horizon, charsize=1.7
cgtext, 0.5, 0.965, ps_name, charsize=2, charthick=1.6, alignment=0.5, color='black', font=0, /normal
device, /close_file
com = 'convert -flatten -background white '+ps_name+'.ps'+' OUTPUT/'+file_basename(ps_name, '.ps')+'.png'
spawn, com, /hide
file_delete, ps_name+'.ps'
-
히마와리 8호 (Himawari-8/AHI) 자료를 이용한 레일라이 산란 보정
-
RGB 천연색을 위한 설정 및 가시화
-
레일라이 산란 보정된 채널 02 및 채널 04를 통해 하이브리드 채널 02로 변환
-
plots를 통해 위/경도에 따른 값 매핑
-
포스트 스크립트 (PS)를 PNG로 변환
-
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> Hybrid RGB
;******************************************************************
ps_name = 'Himawari_Hybrid_RGB_'+time
set_plot, 'ps'
device, filename=ps_name+'.ps', decomposed=1, bits=24, /color, xsize=12, ysize=13, /inches, $
/schoolbook, /bold
!p.font=0 & !p.charsize=2.0 & !p.charthick=1.6 & !p.multi=[0,1,1]
factor = 0.16
hydro_ref02_ray = ((1-factor)*ref02_ray) + (factor*ref04_ray)
scaledBand = scalemodis(ref03_ray, hydro_ref02_ray, ref01_ray)
red = fltarr(dims[0], dims[1])
green = fltarr(dims[0], dims[1])
blue = fltarr(dims[0], dims[1])
red = reform(scaledBand[*,*,0])
green = reform(scaledBand[*,*,1])
blue = reform(scaledBand[*,*,2])
loc1 = where( (red eq 0) and (green le 64) and (blue eq 0), nloc1 )
if (nloc1 gt 0) then begin
red[loc1] = 255
green[loc1] = 255
blue[loc1] = 255
endif
region = [-90, 90, -180, 180] ; All earth
; region = [ 32, 44, 122, 134] ; Hanbando
;region = [ 10, 100, 50, 132] ; East asia
; region = [ 32, 120, 42, 132] ; GOCI (1)
; region = [ 25, 110, 50, 150] ; GOCI (2)
latmin=region[0] & latmax=region[1]
lonmin=region[2] & lonmax=region[3]
clat = -0.007606d & clon = 140.714793d
cgmap_set, /satellite, sat_p=[42159.934825,0,0], clat, clon, charsize=1.0, limit=[latmin,lonmin,latmax,lonmax], $
position=[0.10,0.15,0.9,0.95], /isotropic, color='white'
usersym, [0.1,0.1,-0.1,-0.1,0.1], [-0.1,0.1,0.1,-0.1,-0.1], /fill
for j=0, dims[1]-1 do begin
for i=0, dims[0]-1 do begin
val = red[i,j] + (green[i,j]*(2L^8)) + (blue[i,j]*(2L^16))
plots, lon[i,j], lat[i,j], psym=8, symsize=1.0, color=val
endfor
endfor
; stop
; tvimage, rgb, /true , /overplot, /order
cgmap_gshhs, 'gshhg'+path_sep()+'gshhs_i.b', color='black', /outline, level=4, thick=1.0
cgmap_continents, /continents, /hires, color='black', mlinethick=0.001, thick=1.0, limit=[-90,-180,90,180]
lats = [-90, -60, -30, 0, 30, 60, 90]
; lats_names = ['',textoidl('60\circS', font=0), textoidl('30\circS', font=0), ' ', textoidl('30\circN', font=0),textoidl('60\circN', font=0), '']
lats_names = ['', '', '', '', '', '', '']
lons = [10, 80, 110, 140, 170, -160, -10]
; lons_names = ['',textoidl('80\circE'), textoidl('110\circE'), textoidl('140\circE'), textoidl('170\circE'), textoidl('160\circW'),'']
lons_names = ['', '', '', '', '', '', '']
cgmap_grid, color='black', /label, lats=lats, latlabel=170, lonlabel=-5, latnames=lats_names, lons=lons, lonnames=lons_names, clip_text=0, linestyle=1, noclip=0, /horizon, charsize=1.7
cgtext, 0.5, 0.965, ps_name, charsize=2, charthick=1.6, alignment=0.5, color='black', font=0, /normal
device, /close_file
com = 'convert -flatten -background white '+ps_name+'.ps'+' OUTPUT/'+file_basename(ps_name, '.ps')+'.png'
spawn, com, /hide
file_delete, ps_name+'.ps'
-
히마와리 8호 (Himawari-8/AHI) 자료를 이용한 RGB 천연색
-
가시화에 필요한 라이브러리
pro quick_Ray_crr_single, xvza_mod03, xsza_mod03, xraa_mod03, mscatt, wvl_desnm, refl_inp, refl_out
; This code do a quick Rayleigh correction for MODIS TOA reflectance at red, green, and blue bands.
; Program written by Myeong-Jae Jeong (MJ) at GWNU, Korea
; Version 1.2 May 3, 2013
; All Copyright (C) is reserved.
; Inputs
; xvza_mod03: satellite viewing zenith angle
; xsza_mod03: solar zenith angle
; xraa_mod03: relative azimuth angle
; mscatt : scattering angle
; wvldesnm : center wavelength of input band
; refl_inp : input TOA reflectance
; Output
; refl_out : Rayleigh scattering effect removed reflectance
fillvnan=!values.f_nan
; quick and dirty Rayleigh Correction...
;compute scattering angle
mthet=xvza_mod03 & mthet0=xsza_mod03 & mphi=xraa_mod03 ; MODIS-obs geometry (2-D)
;get_scatt_angle_only, mthet0, mthet, mphi, mscatt
scat_ang=mscatt ; scattering angle
wvl_desum=wvl_desnm/1000.0 ; wvl unit conversion (nm --> um)
;Rayleigh Phase function & Rayleigh ;reflectance
;get_rayleigh_od, wvl_des, ray_od_des, depol_des
rexp=-4.08 & ray_od_test=0.008735*wvl_desum^rexp ; based on Iqbal(1983)
ray_od_des=ray_od_test
pha_ray=0.75*(1.0+cos(!dtor*scat_ang)*cos(!dtor*scat_ang))
ssa_ray=1.0 ; SSA for Air moleculs
;mu_0r=(cos(!dtor*mthet0)+0.15*(93.885-mthet0)^(-1.253))
;mu_vr=(cos(!dtor*mthet)+0.15*(93.885-mthet)^(-1.253))
mu_0=cos(!dtor*mthet0) & mu_v=cos(!dtor*mthet)
mu_0r=mu_0*1.0D+0 & mu_vr=mu_v*1.0D+0
; Calculates Rayleigh reflectance based on single-scattering approximation.
;refl_ray_sscat=ray_od_des*ssa_ray*pha_ray/(4.0*mu_0*mu_v) ; approx
cfim=1.0D+0
refl_ray=ssa_ray*pha_ray/(4.0D+0*(mu_0r+mu_vr))*(1.0D+0 - exp(-1.0D+0*ray_od_des*(1.0D+0/mu_vr + 1.0D+0/mu_0r)))/cfim
refl_out = refl_inp - refl_ray
;save memory
refl_ray=fillvnan
return
end
;+
; NAME:
; SCALEMODIS
;
; PURPOSE:
;
; MODIS corrected reflectance images often appear drab when initially processed
; and displayed on a computer using BYTSCL. In fact, the resulting true-color
; images look nothing like the images you can find on the MODIS Rapid Response
; web page (http://rapidfire.sci.gsfc.nasa.gov/gallery/). After poking around on
; the Internet for awhile, I discovered that the Rapid Response Team doesn't use
; BYTSCL to prepare the images. Rather, they selectively scale portions of the
; reflectance image, using a piecewise scaling with different slopes. This program
; implements this Rapid Response Team algorithm.
;
; AUTHOR:
;
; FANNING SOFTWARE CONSULTING
; David Fanning, Ph.D.
; 1645 Sheely Drive
; Fort Collins, CO 80526 USA
; Phone: 970-221-0438
; E-mail: davidf@dfanning.com
; Coyote's Guide to IDL Programming: http://www.dfanning.com/
;
; CATEGORY:
;
; Graphics
;
; CALLING SEQUENCE:
;
; scaledBand = ScaleModis(red, green, blue)
;
; ARGUMENTS:
;
; red: A two-dimensional array representing the corrected reflectance
; values of a MODIS image. This is a required parameter. If the
; green and blue parameters are also used, this parameter will
; represent the red band of a RGB 24-bit scaled image that is returned.
;
; green: If the three parameters--red, green, and blue--are present, the returned
; array is a 24-bit true-color image, scaled appropriately. This parameter
; is used as the green band in such an image. The parameter is a two-dimensional
; array of corrected reflectance values.
;
; blue: If the three parameters--red, green, and blue--are present, the returned
; array is a 24-bit true-color image, scaled appropriately. This parameter
; is used as the blue band in such an image. The parameter is a two-dimensional
; array of corrected reflectance values.
;
; KEYWORD PARAMETERS:
;
; RANGE: A two-dimensional array that the input bands are first scaled into, prior to
; the differential scaling using the MODIS Rapid Response algorithm. The default
; input range is [-0.01, 1.10]. These values will be used to set the MIN and MAX
; keywords for the BYTSCL command in the initial scaling of the input bands.
;
; CLOUD: The MODIS Rapid Response team uses a slightly different scaling algorithm when
; the idea is to emphasize clouds in a MODIS scene. Set this keyword to use the
; alternate cloud scaling algorithm.
;
; OUTPUTS:
;
; scaledBand: If a single 2D array is passed as the argument, then scaledBand is the scaled
; 2D output array. If all three arguments are passed to the program, then scaledBand
; is a scaled 24-bit image that represents a true-color or false color representation
; of the three input bands.
;
; MODIFICATION HISTORY:
;
; Written by: David W. Fanning, July 2009, using the IDL programs MODIS_FALSE_COLOR and
; and SCALE_IMAGE for inspiration. I found these programs on the Internet when poking
; around MODIS web pages. I suspect, but I am not sure, these programs were originally
; written by Liam Gumley.
; Minor changes to the ScaleIt function to be sure partitioning is done correctly. 5 Aug 2009. DWF.
;-
;
;******************************************************************************************;
; Copyright (c) 2009, by Fanning Software Consulting, Inc. ;
; All rights reserved. ;
; ;
; Redistribution and use in source and binary forms, with or without ;
; modification, are permitted provided that the following conditions are met: ;
; ;
; * Redistributions of source code must retain the above copyright ;
; notice, this list of conditions and the following disclaimer. ;
; * Redistributions in binary form must reproduce the above copyright ;
; notice, this list of conditions and the following disclaimer in the ;
; documentation and/or other materials provided with the distribution. ;
; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
; contributors may be used to endorse or promote products derived from this ;
; software without specific prior written permission. ;
; ;
; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
;******************************************************************************************;
FUNCTION ScaleModis_ScaleIt, image, input, output
; This function performs the actual piecewise scaling and returns
; the differentially scaled image.
ON_Error, 2 ; Return to caller.
; Create the output array.
scaled = image * 0B
; Check the input vector lengths to be sure they are the same.
inum = N_Elements(input)
onum = N_Elements(output)
IF inum NE onum THEN Message, 'Scaling vectors must be the same length.'
; Partition the image into input values.
h = Histogram(Value_Locate(input[1:inum-2], image) + 2, REVERSE_INDICES=ri)
; Differentially scale the image.
FOR index=0,N_Elements(input)-2 DO BEGIN
x1 = input[index]
x2 = input[index+1]
y1 = output[index]
y2 = output[index+1]
; Find slope and intercept for scaling.
m = (y2 - y1) / Float((x2 - x1))
b = y2 - (m * x2)
; Get the indices you need to scale.
IF ri[index] NE ri[index+1] THEN BEGIN
indices = ri[ri[index]:ri[index+1]-1]
; Scale the image.
scaled[indices] = Byte(m * image[indices] + b)
ENDIF
ENDFOR
RETURN, scaled
END ; ---------------------------------------------------------------------------------------
FUNCTION ScaleModis, red, grn, blu, RANGE=range, CLOUD=cloud
; Error handling. Return to caller with error condition.
On_Error, 2
; Check keywords.
IF N_Elements(range) EQ 0 THEN range = [-0.01, 1.10]
IF N_Elements(range) EQ 1 THEN Message, 'RANGE must be a two-element array.'
cloud = Keyword_Set(cloud)
; Set up input and output vectors for Rapid Response scaling. Image pixels
; between adjacent values in the input vector are linearly scaled to the corresponding
; values in the output vector. In other words, pixels with values between 0 and 30 in
; the byte scaled input image are scaled into the range 0 to 110 in the output image, and
; so on. Each portion of the input image is scaled differently.
IF cloud THEN BEGIN
input = [0, 25, 55, 100, 255]
output = [0, 90, 140, 175, 255]
ENDIF ELSE BEGIN
input = [0, 30, 60, 120, 190, 255]
output = [0, 110, 160, 210, 240, 255]
ENDELSE
; Proceed differently, depending upon the the number of positional parameters.
CASE N_Params() OF
0: Message, 'Must pass a MODIS corrected reflectance 2D image as an argument.'
1: BEGIN
ndims = Size(red, /N_DIMENSIONS)
IF ndims NE 2 THEN Message, 'Input image must be a 2D array.'
scaled = ScaleModis_ScaleIt(BytScl(red, MIN=range[0], MAX=range[1]), input, output)
END
2: Message, 'Either 1 or 3 positional arguments are required in this function.'
3: BEGIN
ndims = Size(red, /N_DIMENSIONS)
IF ndims NE 2 THEN Message, 'Input image must be a 2D array.'
dims = Size(red, /DIMENSIONS)
IF Total(dims EQ Size(grn, /DIMENSIONS)) NE 2 THEN Message, 'Input images must have the same dimensions'
IF Total(dims EQ Size(blu, /DIMENSIONS)) NE 2 THEN Message, 'Input images must have the same dimensions'
scaled = BytArr(dims[0], dims[1], 3)
scaled[0,0,0] = ScaleModis_ScaleIt(BytScl(red, MIN=range[0], MAX=range[1]), input, output)
scaled[0,0,1] = ScaleModis_ScaleIt(BytScl(grn, MIN=range[0], MAX=range[1]), input, output)
scaled[0,0,2] = ScaleModis_ScaleIt(BytScl(blu, MIN=range[0], MAX=range[1]), input, output)
END
ENDCASE
RETURN, scaled
END ; ---------------------------------------------------------------------------------------
end
[전체]
참고 문헌
[논문]
- 없음
[보고서]
- 없음
[URL]
- 없음
문의사항
[기상학/프로그래밍 언어]
- sangho.lee.1990@gmail.com
[해양학/천문학/빅데이터]
- saimang0804@gmail.com
반응형
최근댓글