반응형

     정보

    • 업무명    : 아이디엘 히마와리 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

     

    • 물리적 값으로 산출하기 위한 디지털 값에 대한 변환 계수

    calibration_coeff.txt
    0.00MB

     

    • 히마와리위성 8호에 대한 위/경도 자료

    GeoData.zip
    3.85MB

     

    • 천리안위성 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

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

     

    본 블로그는 파트너스 활동을 통해 일정액의 수수료를 제공받을 수 있음
    반응형
    • 네이버 블러그 공유하기
    • 네이버 밴드에 공유하기
    • 페이스북 공유하기
    • 카카오스토리 공유하기