정보
-
업무명 : 아이디엘 히마와리 8호 (Himawari-8/AHI) 기상 위성 자료를 이용한 RGB 합성 및 천리안 1호 (COMS/MI) 기상 위성 자료로부터 에어로졸 광학두께 (AOD) 중첩
-
작성자 : 이상호
-
작성일 : 2019-08-29
-
설 명 :
-
수정이력 :
-
2020-07-09 : 내용 보완
-
내용
[개요]
-
안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
-
오늘 포스팅은 히마와리 8호 (Himawari-8/AHI) 기상 위성 자료를 이용한 RGB 합성 및 천리안 1호 (COMS/MI) 기상 위성 자료로부터 에어로졸 광학두께 (AOD) 중첩를 소개해 드리고자 합니다.
[특징]
-
천리안 위성 1호 (COMS/MI) 및 히마와리 위성 8호 (Himawari-8/AHI) 자료에 대한 합성 및 중첩 기술이 요구되며 이 프로그램은 이러한 목적을 달성하기 위해 고안된 소프트웨어
[기능]
-
Himawari-8/AHI 자료에 대한 원시 자료 읽기 및 반사율로 변환
-
RGB 합성 영상을 위한 태양 천정각 (SZA), 태양 방위각 (SAA), 위성 천정각 (VZA), 위성 방위각 (VAA), 상대방위각 (RAA) 계산
-
각 채널별 반사율을 이용하여 적색 (Red), 녹색 (Green), 청색 (Blue)으로 변환
-
COMS/MI 자료에 대한 원시자료 읽기 및 에어로졸 광학두께로 변환
-
RGB 합성 후에 에어로졸 광학두께 매핑
[활용 자료]
-
자료명 : himawari8_ahi_le1b_[b01, b02, b03, b04]_8km_f_201508010000.bin
-
자료 종류 : 기본영상 가시채널 4종 (CH01-04)
-
확장자 : bin
-
영역 : 전구
-
기간 : 2015년 01월 01일 00:00 UTC
-
공간 해상도 : 8 km
himawari8_ahi_le1b_b01_8km_f_201508010000.bin
3.61MB
himawari8_ahi_le1b_b02_8km_f_201508010000.bin
3.61MB
himawari8_ahi_le1b_b03_8km_f_201508010000.bin
3.61MB
himawari8_ahi_le1b_b04_8km_f_201508010000.bin
3.61MB
-
물리적 값으로 산출하기 위한 디지털 값에 대한 변환 계수
-
히마와리위성 8호에 대한 위/경도 자료
-
천리안위성 1A에 대한 위/경도 및 에어로졸 광학두께 자료
ComsData.zip
8.69MB
ComsData.z01
10.00MB
[자료 처리 방안 및 활용 분석 기법]
-
없음
[사용법]
-
천리안 위성 1호 및 히마와리 위성 8호 자료를 동일 디렉터리에 저장
-
소스 코드를 실행 (idl -e RGB_Composite_Image_and_AOD_Using_Himawari8_Data.pro)
-
가시화 결과를 확인
[사용 OS]
-
Linux (CentOS v7.3)
-
Windows 10
[사용 언어]
-
IDL v8.5
소스 코드
-
해당 작업을 하기 위한 컴파일 및 실행 코드는 IDL로 작성되었으며 시각화를 위한 라이브러리는 Coyote Library를 이용하였습니다.
-
작업 환경의 경우 3개 디렉터리로 구성되며 주 프로그램은 RGB_Composite_Image_and_AOD_ Using_Himawari8_Data.pro입니다.
-
소스 코드는 단계별로 수행하며 상세 설명은 다음과 같습니다.
-
1 단계는 기본 정보 설정
-
[2 단계] 입력 자료 읽기 및 물리적 변수로 환산
-
[3 단계] 각도 정보 계산
-
[4 단계] 채널별 반사율로 환산
-
[5 단계] RGB 합성 영상
-
[6 단계] RGB 합성 영상 (레일라이 산란 보정)
-
[7 단계] RGB 합성 영상 (하이브리드)
-
[8 단계] RGB 합성 영상 (하이브리드) + AOD 중첩
-
[명세]
-
[1 단계] 기본 정보 설정
-
dim : 배열
-
datpath, cofpath : 입력자료 경로
-
count01-04, lat, lon, SZA : 변수 선언
-
세부 내용은 변수명을 참조하시기 바랍니다.
-
dim = [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_'
lon_file = datpath + path_sep() + 'lon' + '.dat'
lat_file = datpath + path_sep() + 'lat' + '.dat'
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 = intarr(dim[0], dim[1]) & count02 = intarr(dim[0], dim[1]) ; Count
count03 = intarr(dim[0], dim[1]) & count04 = intarr(dim[0], dim[1])
lat = fltarr(dim[0], dim[1]) ; latitude
lon = fltarr(dim[0], dim[1]) ; longitude
SZA = fltarr(dim[0], dim[1]) ; solor zenith angle
-
[2 단계] 입력 자료 읽기 및 물리적 변수로 환산
-
openr를 통해 각 채널에 대한 디지털값 (count01-04)
-
lat, lon : 위/경도 정보
-
채널별 디지털값의 경우 물리적 상수 (rad01-04)를 통해 복사 휘도로 환산
-
;******************************************************************
; 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, lat_file, /get_lun & readf, lun, lat & free_lun, lun
openr, lun, lon_file, /get_lun & readf, lun, lon & free_lun, lun
;******************************************************************
; 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
-
[3 단계] 각도 정보 계산
-
NASA 웹 사이트에서 IDL 소스 코드를 참조
-
RGB 영상을 게산하기 위한 각도 정보 계산
-
SZA : 태양 천정각
-
SAA : 태양 방위각
-
VZA : 위성 천정각
-
VAA : 위성 방위각
-
RAA : 상대 방위각
-
SA : 산란각
-
GA : 태양광 반사각
-
-
where를 통해 일정 각도 이상인 경우 결측값 처리
-
;******************************************************************
; 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
;====================================
; Angle Calc
;=====================================
SZA = solzen
SAA = az_not
VZA = theta
VAA = azimuth
RAA = fltarr(dim[0], dim[1])
SA = fltarr(dim[0], dim[1])
GA = fltarr(dim[0], dim[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
;******************************************************************
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
-
[4 단계] 채널별 반사율로 환산
-
각 채널에 따른 복사 휘도를 통해 반사율로 환산
-
where를 통해 일정 범위값 (최대값, 최소값) 이외의 경우 결측값 처리
-
;******************************************************************
; 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
-
[5 단계] RGB 합성 영상
-
scalemodis 및 채널별 반사율을 통해 RGB 합성
-
필수 전달 인자 및 fnMakePsPlot를 통해 RGB 합성 영상 이미지 저장
-
필수 전달 인자에 대한 세부 내용은 다음과 같다.
-
ps_name : 이미지 저장 파일명
-
main_name : 이미지 제목
-
color_name : 컬러바
-
-
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> True RGB
;******************************************************************
ps_name = "Img01_" + time
main_name = "Himawari8_True_RGB : " + time
color_name = ""
scaledBand = scalemodis(ref03, ref02, ref01)
fnMakePsPlot, lon, lat, scaledBand, dim, "false", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
-
[6 단계] RGB 합성 영상 (레일라이 산란 보정)
-
quick_Ray_crr_single를 통해 채널별 반사율에서 레일라이 산란 보정 수행
-
이 경우 where를 통해 일정 범위값 (최대값, 최소값) 이외의 결측값 처리
-
앞서 동일하게 fnMakePsPlot를 통해 RGB 합성 영상 (레일라이 산란 보정) 이미지 저장
-
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> RGB Rayleigh Correction
;******************************************************************
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
ps_name = "Img02_" + time
main_name = "Himawari8_RGB_Rayleigh_Correction : " + time
color_name = ""
scaledBand = scalemodis(ref03_ray, ref02_ray, ref01_ray)
fnMakePsPlot, lon, lat, scaledBand, dim, "false", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
-
[7 단계] RGB 합성 영상 (하이브리드)
-
레일라이 산란에 의해 보정된 반사율 (채널 02, 04)을 이용하여 하이브리드 채널 02 계산
-
scalemodis 및 fnMakePsPlot를 통해 RGB 합성 영상 (하이브리드) 이미지 저장
-
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> Hybrid RGB
;******************************************************************
factor = 0.16
hydro_ref02_ray = ((1-factor)*ref02_ray) + (factor*ref04_ray)
scaledBand = scalemodis(ref03_ray, hydro_ref02_ray, ref01_ray)
ps_name = "Img03_" + time
main_name = "Himawari8_Hybrid_RGB : " + time
color_name = ""
fnMakePsPlot, lon, lat, scaledBand, dim, "false", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
-
[8 단계] RGB 합성 영상 (하이브리드) + AOD 중첩
-
fnMakePsPlot에서 전달 인자 ("true")를 통해 RGB 합성 영상 (하이브리드)뿐만 아니라 천리안위성 1A호 (COMS/MI)의 에어로졸 광학두께 (AOD) 중첩 이미지 영상 저장
-
;******************************************************************
; Hybrid RGB + COMS AOD
;******************************************************************
factor = 0.16
hydro_ref02_ray = ((1-factor)*ref02_ray) + (factor*ref04_ray)
scaledBand = scalemodis(ref03_ray, hydro_ref02_ray, ref01_ray)
ps_name = "Img04_" + time
main_name = "Himawari8_Hybrid_RGB + COMS_AOD : " + time
color_name = "Aerosol Optical Depth"
fnMakePsPlot, lon, lat, scaledBand, dim, "true", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
-
사용자 편의성 함수 정의
-
전달 인자
-
2차원 위도 (lon), 2차원 경도 (lat), 3차원 값 (scaleBand), 배열 정보 (dim), COMS 영상 중첩 여부 (isComsPlot), COMS 중심 경도 (coms_lon), COMS 중심 위성 (coms_lat), COMS 값 (coms_val), COMS 배열 정보 (coms_dim), 이미지 저장 파일명 (ps_name), 제목명 (main_name), 컬러바명 (color_name)
-
-
기능
-
plot를 통해 위/경도에 따른 값 매핑
-
포스트 스크립트 (PS)를 PNG로 변환
-
-
;====================================================================
; Subroutine : Function Make Postscript to Png Graph
;====================================================================
pro fnMakePsPlot, lon, lat, scaledBand, dim, isComsPlot, coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
set_plot, 'ps'
device, filename = ps_name + ".ps", decomposed = 1, bits = 24, /color, xsize = 12, ysize = 13, /inches, /schoolbook, /bold
n_color = 256
start_color = 0 & end_color = 255 & colorn = end_color - start_color + 1
!p.font=0 & !p.charsize=2.0 & !p.charthick=1.6 & !p.multi=[0,1,1]
red = fltarr(dim[0], dim[1])
green = fltarr(dim[0], dim[1])
blue = fltarr(dim[0], dim[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 i = 0L, dim[0] - 1 do begin
for j = 0L, dim[1] - 1 do begin
nLon = lon[i,j]
nLat = lat[i,j]
nVal = red[i,j] + (green[i,j]*(2L^8)) + (blue[i,j]*(2L^16))
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 (nVal le 0 or finite(nVal, /nan) eq 1) then continue
plots, nLon, nLat, psym = 8, symsize = 1, color = nVal
endfor
endfor
if (isComsPlot eq "true") then begin
device, decomposed = 0
zmin = 0.0 & zmax = 1.0
cgloadct, 33, bottom=0, ncolor=n_color
for i = 0L, coms_dim[0] - 1 do begin
for j = 0L, coms_dim[1] - 1 do begin
nComsLon = coms_lon[i,j]
nComLat = coms_lat[i,j]
nComsVal = coms_val[i,j]
if (lonmin gt nComsLon or nComsLon gt lonmax or finite(nComsLon, /nan) eq 1) then continue
if (latmin gt nComLat or nComLat gt latmax or finite(nComLat, /nan) eq 1) then continue
if (nComsVal le 0 or finite(nComsVal, /nan) eq 1) then continue
plots, nComsLon, nComLat, psym = 8, symsize = 2.5, color = BYTSCL(nComsVal, zmin, zmax)
endfor
endfor
cgcolorbar, range=[zmin,zmax], bottom=0, ncolor=n_color, position=[0.15,0.09,0.85,0.13], $
yminor=0, ytickinterval=0, xminor=5, xstyle=1+8, ystyle=4, $
xthick=5, ythick=5, ticklen=0.29999, charsize=1.9
cgtext, 0.5, 0.02, color_name, charsize=1.9, charthick=1.0, alignment=0.5, color='black', font=0, /normal
endif
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, main_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, /hide, com
; file_delete, ps_name + ".ps"
return
end
-
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
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 ; ---------------------------------------------------------------------------------------
[전체]
-
RGB_Composite_Image_and_AOD_Using_Himawari8_Data.pro
PRO RGB_Composite_Image_and_AOD_Using_Himawari8_Data
dim = [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_'
lon_file = datpath + path_sep() + 'lon' + '.dat'
lat_file = datpath + path_sep() + 'lat' + '.dat'
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 = intarr(dim[0], dim[1]) & count02 = intarr(dim[0], dim[1]) ; Count
count03 = intarr(dim[0], dim[1]) & count04 = intarr(dim[0], dim[1])
lat = fltarr(dim[0], dim[1]) ; latitude
lon = fltarr(dim[0], dim[1]) ; longitude
SZA = fltarr(dim[0], dim[1]) ; solor zenith angle
;******************************************************************
; 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, lat_file, /get_lun & readf, lun, lat & free_lun, lun
openr, lun, lon_file, /get_lun & readf, lun, lon & free_lun, lun
;******************************************************************
; 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
;******************************************************************
; 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
;====================================
; Angle Calc
;=====================================
SZA = solzen
SAA = az_not
VZA = theta
VAA = azimuth
RAA = fltarr(dim[0], dim[1])
SA = fltarr(dim[0], dim[1])
GA = fltarr(dim[0], dim[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
;====================================
; COMS MI ADO File read
;=====================================
hdf5_file = datpath + path_sep() + 'coms_mi_le2_aod_cn_201402010000.h5'
; h5_browser(hdf5_file)
data_name = '/Product/Aerosol_Optical_Depth'
file_id = h5f_open(hdf5_file)
data_id = h5d_open(file_id, data_name)
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])
coms_val = data
coms_dim = size(coms_val, /dimensions)
hdf4_file = datpath + path_sep() + 'coms_cn_geos_lonlat.hdf'
; hdf_browser(hdf4_file)
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]
coms_lat = 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]
coms_lon = data
hdf_sd_endaccess, sds_id
hdf_sd_end, sd_id
;******************************************************************
; 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
;******************************************************************
; 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
;============================
; Statistical analysis
;============================
; print, min(lon, /nan), mean(lon, /nan), max(lon, /nan)
; print, min(lat, /nan), mean(lat, /nan), max(lat, /nan)
; print, min(var, /nan), mean(var, /nan), max(var, /nan)
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> True RGB
;******************************************************************
ps_name = "Img01_" + time
main_name = "Himawari8_True_RGB : " + time
color_name = ""
scaledBand = scalemodis(ref03, ref02, ref01)
fnMakePsPlot, lon, lat, scaledBand, dim, "false", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> RGB Rayleigh Correction
;******************************************************************
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
ps_name = "Img02_" + time
main_name = "Himawari8_RGB_Rayleigh_Correction : " + time
color_name = ""
scaledBand = scalemodis(ref03_ray, ref02_ray, ref01_ray)
fnMakePsPlot, lon, lat, scaledBand, dim, "false", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
;******************************************************************
; Himawari-8 Reflectance 1, 2, 3 -> Hybrid RGB
;******************************************************************
factor = 0.16
hydro_ref02_ray = ((1-factor)*ref02_ray) + (factor*ref04_ray)
scaledBand = scalemodis(ref03_ray, hydro_ref02_ray, ref01_ray)
ps_name = "Img03_" + time
main_name = "Himawari8_Hybrid_RGB : " + time
color_name = ""
fnMakePsPlot, lon, lat, scaledBand, dim, "false", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
;******************************************************************
; Hybrid RGB + COMS AOD
;******************************************************************
factor = 0.16
hydro_ref02_ray = ((1-factor)*ref02_ray) + (factor*ref04_ray)
scaledBand = scalemodis(ref03_ray, hydro_ref02_ray, ref01_ray)
ps_name = "Img04_" + time
main_name = "Himawari8_Hybrid_RGB + COMS_AOD : " + time
color_name = "Aerosol Optical Depth"
fnMakePsPlot, lon, lat, scaledBand, dim, "true", coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
end
;====================================================================
; Subroutine : Function Make Postscript to Png Graph
;====================================================================
pro fnMakePsPlot, lon, lat, scaledBand, dim, isComsPlot, coms_lon, coms_lat, coms_val, coms_dim, ps_name, main_name, color_name
set_plot, 'ps'
device, filename = ps_name + ".ps", decomposed = 1, bits = 24, /color, xsize = 12, ysize = 13, /inches, /schoolbook, /bold
n_color = 256
start_color = 0 & end_color = 255 & colorn = end_color - start_color + 1
!p.font=0 & !p.charsize=2.0 & !p.charthick=1.6 & !p.multi=[0,1,1]
red = fltarr(dim[0], dim[1])
green = fltarr(dim[0], dim[1])
blue = fltarr(dim[0], dim[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 i = 0L, dim[0] - 1 do begin
for j = 0L, dim[1] - 1 do begin
nLon = lon[i,j]
nLat = lat[i,j]
nVal = red[i,j] + (green[i,j]*(2L^8)) + (blue[i,j]*(2L^16))
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 (nVal le 0 or finite(nVal, /nan) eq 1) then continue
plots, nLon, nLat, psym = 8, symsize = 1, color = nVal
endfor
endfor
if (isComsPlot eq "true") then begin
device, decomposed = 0
zmin = 0.0 & zmax = 1.0
cgloadct, 33, bottom=0, ncolor=n_color
for i = 0L, coms_dim[0] - 1 do begin
for j = 0L, coms_dim[1] - 1 do begin
nComsLon = coms_lon[i,j]
nComLat = coms_lat[i,j]
nComsVal = coms_val[i,j]
if (lonmin gt nComsLon or nComsLon gt lonmax or finite(nComsLon, /nan) eq 1) then continue
if (latmin gt nComLat or nComLat gt latmax or finite(nComLat, /nan) eq 1) then continue
if (nComsVal le 0 or finite(nComsVal, /nan) eq 1) then continue
plots, nComsLon, nComLat, psym = 8, symsize = 2.5, color = BYTSCL(nComsVal, zmin, zmax)
endfor
endfor
cgcolorbar, range=[zmin,zmax], bottom=0, ncolor=n_color, position=[0.15,0.09,0.85,0.13], $
yminor=0, ytickinterval=0, xminor=5, xstyle=1+8, ystyle=4, $
xthick=5, ythick=5, ticklen=0.29999, charsize=1.9
cgtext, 0.5, 0.02, color_name, charsize=1.9, charthick=1.0, alignment=0.5, color='black', font=0, /normal
endif
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, main_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, /hide, com
; file_delete, ps_name + ".ps"
return
end
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
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 ; ---------------------------------------------------------------------------------------
참고 문헌
[논문]
- 없음
[보고서]
- 없음
[URL]
- 없음
문의사항
[기상학/프로그래밍 언어]
- sangho.lee.1990@gmail.com
[해양학/천문학/빅데이터]
- saimang0804@gmail.com
본 블로그는 파트너스 활동을 통해 일정액의 수수료를 제공받을 수 있음
'프로그래밍 언어 > IDL' 카테고리의 다른 글
| [IDL] 아이디엘 대기과학 전공자를 위한 IDL (Interactive Data Language) 소개 (0) | 2020.03.24 |
|---|---|
| [IDL] 아이디엘 이미지 처리 (Image Processing) 교육 연수 (0) | 2020.03.23 |
| [IDL] 아이디엘 히마와리 8호 (Himawari-8/AHI) 기상 위성 자료를 이용한 RGB 합성영상 가시화 (0) | 2019.08.30 |
| [IDL] 아이디엘 통계 결과를 테일러 다이어그램 (Taylor Diagram) 가시화 (0) | 2019.07.25 |
| [IDL] 아이디엘 Terra/CERES 극궤도 기상 위성 자료를 이용하여 지표면 특성 (Landcover) 및 운량 (Cloud Fraction)에 따른 가시화 (0) | 2019.07.25 |
최근댓글