정보
-
업무명 : 아이디엘 HDF 형식인 Aqua·Terra/MODIS 극궤도 기상 위성 자료를 이용한 에어로졸 광학두께 가시화
-
작성자 : 이상호
-
작성일 : 2020-03-25
-
설 명 :
-
수정이력 :
내용
[개요]
-
안녕하세요? 기상 연구 및 웹 개발을 담당하고 있는 해솔입니다.
-
대학원 석사 1학기에 배운 지구물리 원격탐사에 대한 실습 내용을 다루고자 합니다.
-
오늘 포스팅에서는 HDF 형식인 Aqua·Terra/MODIS 극궤도 기상 위성 자료를 이용한 에어로졸 광학두께 가시화를 소개해 드리고자 합니다.
-
추가로 지구물리 원격탐사에 대한 이론적 배경을 소개한 링크를 보내드립니다.
[강릉원주대 대기환경과학과] 2015년 2학기 전선 지구물리 원격탐사 소개 및 과제물
정보 업무명 : 2015년 2학기 전선 지구물리 원격탐사 소개 및 과제물 작성자 : 이상호 작성일 : 2019-12-20 설 명 : 수정이력 : 내용 [특징] 지구물리 원격탐사 수업에 대한 이해를 돕기위해 작성 [기능] 소개 주..
shlee1990.tistory.com
[특징]
-
HDF 극궤도 기상 위성 자료를 이용하여 정성적으로 가시화하기 위한 도구가 필요하며 이 프로그램은 이러한 목적을 달성하기 위해 고안된 소프트웨어
[기능]
-
Aqua/MODIS 및 Terra/MODIS 극궤도 기상 위성 자료에 대한 위/경도 및 에어로졸 광학두께 읽기
-
에어로졸 광학두께 가시화 및 이미지 저장
[활용 자료]
-
자료명 : MYD04_L2.A2012289.0455.051.2012289222139.hdf
MOD04_L2.A2003142.1730.005.2006352204928.hdf -
자료 종류 : 에어로졸 광학두께
-
확장자 : HDF
-
영역 : 전지구
-
기간 : 2012년 10월 15일 17:00 UTC
2003년 04월 22일 17:30 UTC -
시간 해상도 : 하루에 2번 관측
-
제공처 : NASA
MYD04_L2.A2012289.0455.051.2012289222139.hdf
1.57MB
MOD04_L2.A2003142.1730.005.2006352204928.hdf
1.30MB
[사용법]
-
입력 자료를 동일 디렉터리에 위치
-
소스 코드를 실행 (idl -e Visualization_Using_MODIS_Aerosol_Data)
-
가시화 결과를 확인
[사용 OS]
-
Linux
-
Windows 10
[사용 언어]
-
IDL v8.5
소스 코드
-
해당 작업을 하기 위한 컴파일 및 실행 코드는 IDL로 작성되었으며 가시화를 위한 라이브러리는 CIMSS/SSEC Library를 이용하였습니다.
-
소스 코드는 단계별로 수행하며 상세 설명은 다음과 같습니다.
-
1 단계는 주 프로그램은 작업 경로 설정, HDF 파일 읽기, 변수 전처리하여 메모리상에 저장하고 가시화를 위한 초기 설정합니다.
-
일반적으로 포스트 스크립트 (PS) 형식에서 PNG로 변환하나 오랜 시간이 소요되는 단점이 있습니다. 따라서 2 단계에서는 imagemap를 통해 영상 표출하여 이미지 형식으로 저장합니다.
-
[명세]
-
전역 변수 설정
-
결측값 (fillvnan), 최소값 (small), 그림 창 크기 (nxsize, nysize)
-
입력 자료 선택 옵션 (iop_aqua)에 따라 설정
-
즉 0인 경우 Terra와 Aqua 위성이고 1 및 2에서는 각각 Aqua 위성 및 Terra 위성을 의미
-
-
fillvnan=!values.f_nan
small=1.0e-6
nxsize=800 & nysize=800
iop_aqua=0
case iop_aqua of
0: begin
filterstr='M*.hdf' ; Terra and Aqua MODIS data
end
1: begin
filterstr='MYD04*.hdf' ; Aqua/MODIS Aerosol data only
end
2: begin
filterstr='MOD04*.hdf' ; Terra/MODIS Aerosol data only
end
else: begin
filterstr='*.hdf' ; any hdf file
end
endcase
-
MODIS (MOD04*) 파일 설정 및 파일 찾기
-
dialog_pickfile를 통해 별도의 검색 창에서 선택
-
pname='C:\SYSTEM\PROG\IDL\RSEE_idl_2_read_plot_MOD04_L2\'
pathout=pname
print, "pname : " + pname
pn_modin=dialog_pickfile(path=pname, filter=filterstr)
-
파일명 (pn_modis)에서 패턴 추출 및 이미지 저장 파일 (fn_img0) 설정
grnuleidstr=strmid(pn_modin,34,13,/reverse_offset)
satidstr=strmid(pn_modin,43,3,/reverse_offset)
iyrmodstr=strmid(pn_modin,33,4,/reverse_offset)
reads, iyrmodstr, iyrin, format='(i4)'
iop_deepblue=1
iop_img=1
fn_img0=pname+'Visualization_Using_MODIS_Aerosol_Data'+satidstr+'04_'+grnuleidstr+'_img0.png'
-
HDF 파일 읽기 및 변수 전처리
-
read_single_MOD04_sds를 통해 육/해상 에어로졸 광학두께 변수 선택
-
target_sds_name_in='Optical_Depth_Land_And_Ocean'
read_single_MOD04_sds, pn_modin, target_sds_name_in, ndim_op_aot550, op_aot550, ndim_latlon, xlat_mod, xlon_mod
-
관심 영역 (ROI)에서 에어로졸 광학두께 평균 계산
locx=where(xlat_mod ge 36. and xlat_mod lt 38. and xlon_mod ge 125. and xlon_mod lt 128., nlocx)
caot = op_aot550[locx]
avg = mean(caot, /nan)
print, 'avg=', avg
-
가시화를 위한 초기 설정
-
폰트 설정
-
loadct를 통해 컬러바 설정
-
!p.charsize=1.3
!p.charthick=2.5
!p.thick=2.5
!p.symsize=3.0
!x.thick=2.5
!x.charsize=1.3
!y.thick=2.5
!y.charsize=1.3
loadct, 39 ; load color table
-
에어로졸 광학두께를 위한 설정 및 가시화
-
imagemap를 통해 위도 (lon), 경도 (lat)에 따른 에어로졸 광학두께 (aod) 맵핑
-
saveimage를 통해 png 형식으로 이미지 저장
-
;*************************************
window, 0, xsize=nxsize, ysize=nysize
device, decomposed=0
;*************************************
mag=5
nxx_mod=ndim_latlon[0] & nyy_mod=ndim_latlon[1]
minx=0.0 & maxx=2.0
plotlimit=[min(xlat_mod),min(xlon_mod),max(xlat_mod),max(xlon_mod)]
ndivbar=5
nclev=255
delclev = (maxx - minx)/nclev
lat=congrid(xlat_mod, mag*nxx_mod, mag*nyy_mod, /interp)
lon=congrid(xlon_mod, mag*nxx_mod, mag*nyy_mod, /interp)
aod=congrid(op_aot550, mag*nxx_mod, mag*nyy_mod,/minus_one)
locfill=where(aod ge maxx-delclev,nlocfill)
if(nlocfill gt 0) then aod[locfill]=maxx-delclev
map_set, limit=plotlimit, /noborder, /cylindrical, xmargin=2, ymargin=5, /isotropic
imagemap, aod, lat, lon, range=[minx,maxx], limit=plotlimit, missing=0 $
, /noborder, /current, title=titlestr, /isotropic, position=[.01,.2,.99,.95]
color = !d.table_size -1
map_continents, color=color, /countries, /coasts, /hires
map_grid, color=color, latdel=latlon_intvl, londel=latlon_intvl, charsize=1.2, charthick=2.0, /box_axes
colorbar, position=[0.2,0.05,0.8,0.08], range=[minx,maxx] $
, /horizontal, background=0, color=255, divisions=ndivbar $
, ticknames=ticknames, format='(f4.1)'
if(iop_img eq 1) then saveimage, fn_img0, /png
-
Aqua/MODIS를 이용한 에어로졸 광학두께
-
Terra/MODIS를 이용한 에어로졸 광학두께
-
가시화에 필요한 라이브러리
;##################################################################################
;##################################################################################
pro read_single_MOD04_sds, pn_modin, target_sds_name_in, ndim_target_sds, target_sds, ndim_latlon, xlat_mod, xlon_mod
;##################################################################################
;##################################################################################
fillvnan=!values.f_nan
small=1.0e-5
;
;=============================================================================
; open Level 2 MODIS aerosol data(MOD04_L2 or MYD04_L2) - hdf file
;-----------------------------------------------------------------------------
print, "pn_modin : ", pn_modin
fidmod=hdf_sd_start(pn_modin)
;-----------------------------------------------------------------------------
;***************************************
; Latitude
;***************************************
index = hdf_sd_nametoindex(fidmod, 'Latitude')
sdsid = hdf_sd_select(fidmod, index)
hdf_sd_getinfo, sdsid, name=name, ndims=ndims, natts=natts, dims=dims
ndim_latlon=dims
if( natts gt 0 ) then begin
idfill = hdf_sd_attrfind(sdsid,'_FillValue')
hdf_sd_attrinfo, sdsid, idfill, data=fillv_modis
idscl = hdf_sd_attrfind(sdsid,'scale_factor')
hdf_sd_attrinfo, sdsid, idscl, data=sf_modis
idscl = hdf_sd_attrfind(sdsid,'valid_range')
hdf_sd_attrinfo, sdsid, idscl, data=vrng_modis
endif
hdf_sd_getdata, sdsid, data
hdf_sd_endaccess, sdsid
xlat_mod=sf_modis[0]*data
bd=where(data eq fillv_modis[0] or data lt vrng_modis[0] or data gt vrng_modis[1],nbd)
if(nbd gt 0) then xlat_mod[bd]=fillvnan
;***************************************
; Longitude
;***************************************
index = hdf_sd_nametoindex(fidmod, 'Longitude')
sdsid = hdf_sd_select(fidmod, index)
hdf_sd_getinfo, sdsid, name=name, ndims=ndims, natts=natts, dims=dims
if( natts gt 0 ) then begin
idfill = hdf_sd_attrfind(sdsid,'_FillValue')
hdf_sd_attrinfo, sdsid, idfill, data=fillv_modis
idscl = hdf_sd_attrfind(sdsid,'scale_factor')
hdf_sd_attrinfo, sdsid, idscl, data=sf_modis
idscl = hdf_sd_attrfind(sdsid,'valid_range')
hdf_sd_attrinfo, sdsid, idscl, data=vrng_modis
endif
hdf_sd_getdata, sdsid, data
hdf_sd_endaccess, sdsid
xlon_mod=sf_modis[0]*data
bd=where(data eq fillv_modis[0] or data lt vrng_modis[0] or data gt vrng_modis[1],nbd)
if(nbd gt 0) then xlon_mod[bd]=fillvnan
;***************************************
; Target SDS
;***************************************
index = hdf_sd_nametoindex(fidmod, target_sds_name_in)
sdsid = hdf_sd_select(fidmod, index)
hdf_sd_getinfo, sdsid, name=name, ndims=ndims, natts=natts, dims=dims
ndim_target_sds=dims
if( natts gt 0 ) then begin
idfill = hdf_sd_attrfind(sdsid,'_FillValue')
hdf_sd_attrinfo, sdsid, idfill, data=fillv_modis
if( natts gt 3 ) then begin
idscl = hdf_sd_attrfind(sdsid,'scale_factor')
hdf_sd_attrinfo, sdsid, idscl, data=sf_modis
idscl = hdf_sd_attrfind(sdsid,'add_offset')
hdf_sd_attrinfo, sdsid, idscl, data=offset_modis
idscl = hdf_sd_attrfind(sdsid,'valid_range')
hdf_sd_attrinfo, sdsid, idscl, data=vrng_modis
endif
endif
hdf_sd_getdata, sdsid, data
hdf_sd_endaccess, sdsid
if( natts gt 3 ) then begin
target_sds=sf_modis[0]*data+offset_modis[0]
bd=where(data eq fillv_modis[0] or data lt vrng_modis[0] or data gt vrng_modis[1],nbd)
if(nbd gt 0) then target_sds[bd]=fillvnan
endif else begin
target_sds=data
bd=where(abs(data-fillv_modis[0]) le small,nbd)
if(nbd gt 0) then target_sds[bd]=fillvnan
endelse
; conserve memory
data=fillvnan
;=============================================================================
; closing the hdf file & the end of data file access
;-----------------------------------------------------------------------------
hdf_sd_end, fidmod
;=============================================================================
;##################################################################################
return
; The end of subroutine reading selected SDS from MOD08M3
end
;##################################################################################
;=========================================================================
; Compute Julian day for a given date (month, day)
function mmdd2jd, ileap, imo, ida
;=========================================================================
case ileap of
0: begin
day2add=[0,31,59,90,120,151,181,212,243,273,304,334]
end
1: begin
day2add=[0,31,60,91,121,152,182,213,244,274,305,335]
end
else: begin
print, 'Error in ILEAP value - out of range'
end
endcase
jday = day2add[imo-1]+ida
return, jday
end
;=========================================================================
PRO SAVEIMAGE, FILE, $
BMP=BMP, PNG=PNG, GIF=GIF, JPEG=JPEG, TIFF=TIFF, $
QUALITY=QUALITY, DITHER=DITHER, CUBE=CUBE, QUIET=QUIET
;- Check arguments
if (n_params() ne 1) then $
message, 'Usage: SAVEIMAGE, FILE'
if (n_elements(file) eq 0) then $
message, 'Argument FILE is undefined'
if (n_elements(quality) eq 0) then quality = 75
;- Get output file type
output = 'JPEG'
if keyword_set(bmp) then output = 'BMP'
if keyword_set(png) then output = 'PNG'
if keyword_set(gif) then output = 'GIF'
if keyword_set(jpeg) then output = 'JPEG'
if keyword_set(tiff) then output = 'TIFF'
;- Check if GIF output is available
version = float(!version.release)
if (version ge 5.4) and (output eq 'GIF') then $
message, 'GIF output is not available'
;- Get contents of graphics window,
;- and color table if needed
image = screenread(depth=depth)
if (depth le 8) then tvlct, r, g, b, /get
;- Write 8-bit output file
if (output eq 'BMP') or $
(output eq 'PNG') or $
(output eq 'GIF') then begin
;- If image depth is 24-bit, convert to 8-bit
if (depth gt 8) then begin
case keyword_set(cube) of
0 : image = color_quan(image, 1, r, g, b, $
colors=256, dither=keyword_set(dither))
1 : image = color_quan(image, 1, r, g, b, cube=6)
endcase
endif
;- Reverse PNG image order if required
if (output eq 'PNG') and (version le 5.3) then $
image = reverse(temporary(image), 2)
;- Save the image
case output of
'BMP' : write_bmp, file, image, r, g, b
'PNG' : write_png, file, image, r, g, b
'GIF' : write_gif, file, image, r, g, b
endcase
endif
;- Write 24-bit output file
if (output eq 'JPEG') or $
(output eq 'TIFF') then begin
;- Convert 8-bit image to 24-bit
if (depth le 8) then begin
info = size(image)
nx = info[1]
ny = info[2]
true = bytarr(3, nx, ny)
true[0, *, *] = r[image]
true[1, *, *] = g[image]
true[2, *, *] = b[image]
image = temporary(true)
endif
;- Reverse TIFF image order
if (output eq 'TIFF') then $
image = reverse(temporary(image), 3)
;- Write the image
case output of
'JPEG' : write_jpeg, file, image, $
true=1, quality=quality
'TIFF' : write_tiff, file, image, 1
endcase
endif
;- Report to user
if (keyword_set(quiet) eq 0) then $
print, file, output, $
format='("Created ",a," in ",a," format")'
END
FUNCTION SCREENREAD, X0, Y0, NX, NY, DEPTH=DEPTH
;- Check arguments
if (n_elements(x0) eq 0) then x0 = 0
if (n_elements(y0) eq 0) then y0 = 0
if (n_elements(nx) eq 0) then nx = !d.x_vsize - x0
if (n_elements(ny) eq 0) then ny = !d.y_vsize - y0
;- Check for TVRD capable device
tvrd_true = !d.flags and 128
if (tvrd_true eq 0) then message, $
'TVRD is not supported on this device: ' + !d.name
;- On devices which support windows, check for open window
win_true = !d.flags and 256
if (win_true gt 0) and (!d.window lt 0) then message, $
'No graphics window are open'
;- Get IDL version number
version = float(!version.release)
;- Get display depth
depth = 8
if (win_true gt 0) then begin
if (version ge 5.1) then begin
device, get_visual_depth=depth
endif else begin
if (!d.n_colors gt 256) then depth = 24
endelse
endif
;- Set decomposed color mode on 24-bit displays
if (depth gt 8) then begin
entry_decomposed = 0
if (version gt 5.1) then $
device, get_decomposed=entry_decomposed
device, decomposed=1
endif
;- Get the contents of the window
if (depth gt 8) then true = 1 else true = 0
image = tvrd(x0, y0, nx, ny, order=0, true=true)
;- Restore decomposed color mode on 24-bit displays
if (depth gt 8) then device, decomposed=entry_decomposed
;- Return result to caller
return, image
END
pro imagemap, image, lat, lon, newimage = newimage, range = range, $
limit = limit, position = position, isotropic = isotropic, title = title, $
xoffset = xoffset, yoffset = yoffset, xsize = xsize, ysize = ysize, $
missing = missing, noborder = noborder, noerase = noerase, lowres = lowres, $
current = current
;+
;PURPOSE:
; Display an image which has latitude and longitude defined for
; each pixel on a map projection. If a map projection is not currently
; defined, then a Mercator map projection is created which corresponds to
; the lat/lon limits of the image.
;
;CALLING SEQUENCE:
; IMAGEMAP, IMAGE, LAT, LON
;
;INPUT:
; IMAGE Array (2D) or vector (1D) of image values
; LAT Array or vector of latitudes corresponding to image values
; (degrees, -90.0 = S, +90.0 = N)
; LON Array or vector of longitudes corresponding to image values
; (degrees, -180.0 = W, +180.0 = E)
;
;OPTIONAL KEYWORDS:
; NEWIMAGE Named variable in which resampled image array is returned.
; Note that this image is always scaled to a BYTE array.
; RANGE Range of image values used for brightness scaling, [MIN,MAX]
; (default is [MIN(IMAGE),MAX(IMAGE)])
; LIMIT Limits of map display, [LATMIN,LONMIN,LATMAX,LONMAX]
; (default is [MIN(LAT),MIN(LON),MAX(LAT),MAX(LON)])
; POSITION Normalized coordinates for map display window [X1,Y1,X2,Y2]
; (default is to let MAP_SET determine the window size)
; This is useful when used in conjunction with the
; ESRG BOXPOS procedure. For example,
; IMAGEMAP, IMAGE, LAT, LON, POS = BOXPOS( /RMARG )
; will leave room at the right for a COLOR_KEY colorbar.
; ISOTROPIC If set, creates an isotropic map projection (default=non-isotropic).
; TITLE String variable containing image title (default=no title).
; XOFFSET Named variable in which the lower left device X coordinate
; of the displayed image is returned.
; YOFFSET Named variable in which the lower left device Y coordinate
; of the displayed image is returned.
; XSIZE Named variable in which the width of the displayed image
; is returned (used by devices which have scalable pixels
; such as Postscript).
; YSIZE Named variable in which the height of the displayed image
; is returned (used by devices which have scalable pixels
; such as Postscript).
; MISSING Byte value to use for missing (unfilled) portions of image
; (default is zero).
; NOBORDER If set, do not draw border around image (default=draw border).
; NOERASE If set, do not erase window before creating image (default=erase).
; LOWRES If set, draw image in low resolution mode (default=high resolution).
; CURRENT If set, use the current map projection.
;
;OUTPUT:
; The resampled image is displayed in the current graphics window
; in map coordinates. Continental outlines and lat/lon grids may be
; overlaid with MAP_CONTINENTS AND MAP_GRID.
;
;CREATED:
; Liam Gumley, CIMSS/SSEC, 26-JUL-1996
; liam.gumley@ssec.wisc.edu
; http://cimss.ssec.wisc.edu/~gumley/index.html
;
;REVISED:
; Liam Gumley, CIMSS/SSEC, 17-SEP-1996
; Modified to work with Postscript output.
; Added XOFFSET, YOFFSET, XSIZE, YSIZE keywords.
; Mercator map projection is now created only if no map projection exists.
;
; Liam Gumley, CIMSS/SSEC, 15-OCT-1996
; Added MISSING keyword to set missing values in image.
; Added NOBORDER keyword.
;
; Liam Gumley, CIMSS/SSEC, 25-NOV-1996
; Now uses MISSING keyword properly.
; Added NOERASE keyword.
; Added LOWRES keyword (useful for low resolution images, e.g. HIRS, AMSU).
;
; Liam Gumley, CIMSS/SSEC, 26-MAR-1999
; Now uses NOERASE keyword properly.
;
; Liam Gumley, CIMSS/SSEC, 13-AUG-1999
; Added CURRENT keyword.
;
;NOTES:
; (1) Hermann Mannstein (h.mannstein@dlr.de) suggested this IDL method.
; (2) This has been tested on MAS, AVHRR, GOES, and simulated MODIS data.
; It will not work well on low resolution data like HIRS or MSU.
; (3) You might run into problems with data over the poles - I've only
; tried mid-latitude imagery.
; (4) This procedure was designed for display purposes *only*.
; If you use the resampled data for any other purpose, you do so at
; your own risk.
; (5) The example takes about 15.6 seconds to execute
; on a SGI Power Indigo2 (R8000/75MHz) with 256 MB RAM.
;
;EXAMPLE:
;
;; This example is adopted from the ESRG library 'REGRID' procedure.
;; First, create latitude, longitude, and image arrays at 250x200 size.
;;
;c = complex(2,2) + cmplxgen(250,200,/center)/100
;c = c + c^2
;lon = float(c) - 100.0
;lat = 20 + imaginary(c)
;image = sqrt(abs(sin(lon*!pi)*sin(lat*!pi)))^0.3
;;
;; Resize arrays to simulate 1km resolution imagery.
;;
;lat = congrid(lat,1000,800,interp=1)
;lon = congrid(lon,1000,800,interp=1)
;image = congrid(image,1000,800,interp=1)
;;
;; Resample data to Mercator projection, and overlay coastline and grid
;;
;t0 = systime(1.0)
;imagemap,image,lat,lon
;print,'Elapsed time (sec) = ',systime(1.0)-t0
;map_continents
;map_grid
;-
;- check limit keyword
if keyword_set( limit ) then begin
if n_elements( limit ) ne 4 then $
message, 'LIMIT must be a 4 element vector of the form [LATMIN,LONMIN,LATMAX,LONMAX]'
latmin = limit( 0 )
lonmin = limit( 1 )
latmax = limit( 2 )
lonmax = limit( 3 )
endif else begin
latmin = min( lat, max = latmax )
lonmin = min( lon, max = lonmax )
endelse
;- check keywords
if not keyword_set( isotropic ) then isotropic = 0
if not keyword_set( noerase ) then noerase = 0
if n_elements( title ) eq 0 then title = ' '
if n_elements( missing ) eq 0 then missing = 0B
missing = byte( ( missing > 0 ) < ( !d.table_size - 1 ) )
;- create Mercator map projection if necessary after checking position keyword
if not keyword_set( current ) then begin
if n_elements( position ) gt 0 then begin
if n_elements( position ) ne 4 then $
message, 'POSITION must be a 4 element vector of the form [X1,Y1,X2,Y2]'
map_set, /mercator, limit = [ latmin, lonmin, latmax, lonmax ], $
title = title, isotropic = isotropic, position = position, /noborder, $
noerase = noerase
endif else begin
map_set, /mercator, limit = [ latmin, lonmin, latmax, lonmax ], $
title = title, isotropic = isotropic, /noborder, $
noerase = noerase
endelse
endif
;- compute scaling range for byte image after checking range keyword
if keyword_set( range ) then begin
if n_elements( range ) ne 2 then $
message, 'RANGE must be a 2 element vector of the form [MIN,MAX]'
imin = range( 0 )
imax = range( 1 )
endif else begin
imin = min( image, max = imax )
endelse
;- set number of samples and lines for warped image
ns = !d.x_size
nl = !d.y_size
if ( !d.flags and 1 ) then begin
ns = 640L
nl = long( float( ns ) * float( !d.y_size ) / float( !d.x_size ) )
endif
;- create resampled byte image
p = convert_coord( lon, lat, /data, /to_normal )
newimage = replicate( 0B, ns, nl )
newimage( p( 0, * ) * ( ns - 1 ), p( 1, * ) * ( nl - 1 ) ) = $
bytscl( image, min = imin, max = imax, top = !d.table_size - 2 ) + 1B
;- extract portion of image which fits within map boundaries
x = !x.window * ns
y = !y.window * nl
newimage = temporary( newimage( x(0):x(1), y(0):y(1) ) )
;- compute image offset and size (device coordinates)
p = convert_coord( [ x(0), x(1) ] / float( ns ), [ y(0), y(1) ] / float( nl ), $
/normal, /to_device )
xoffset = p(0,0)
yoffset = p(1,0)
xsize = p(0,1) - p(0,0)
ysize = p(1,1) - p(1,0)
;- fill holes in resampled image
nxr = 2
nyr = 2
if keyword_set( lowres ) then begin
nxr = 5
nyr = 5
if ( !d.flags and 1 ) then begin
nxr = 7
nyr = 7
endif
endif
fill = dilate( newimage, replicate( 1, nxr, nyr ), /gray )
loc = where( ( fill ge 1b ) and ( newimage eq 0 ), count )
if count ge 1 then newimage( loc ) = fill( loc )
fill = 0
;- fill remaining undefined areas of image with the missing value
loc = where( newimage eq 0B, count )
if ( count ge 1 ) and ( missing gt 0B) then newimage( loc ) = missing
;- display resampled image
tv, newimage, xoffset, yoffset, xsize = xsize, ysize = ysize
;- draw map borders
if not keyword_set( noborder ) then begin
plots, [ p(0,0), p(0,1) ], [ p(1,0), p(1,0) ], /device
plots, [ p(0,1), p(0,1) ], [ p(1,0), p(1,1) ], /device
plots, [ p(0,1), p(0,0) ], [ p(1,1), p(1,1) ], /device
plots, [ p(0,0), p(0,0) ], [ p(1,1), p(1,0) ], /device
endif
end
;---
;+
; NAME:
; COLORBAR
;
; PURPOSE:
;
; The purpose of this routine is to add a color bar to the current
; graphics window.
;
; AUTHOR:
;
; FANNING SOFTWARE CONSULTING
; David Fanning, Ph.D.
; 1645 Sheely Drive
; Fort Collins, CO 80526 USA
; Phone: 970-221-0438
; E-mail: davidf@dfanning.com
; Coyote's Guide to IDL Programming: http://www.dfanning.com/
;
; CATEGORY:
;
; Graphics, Widgets.
;
; CALLING SEQUENCE:
;
; COLORBAR
;
; INPUTS:
;
; None.
;
; KEYWORD PARAMETERS:
;
; BOTTOM: The lowest color index of the colors to be loaded in
; the bar.
;
; CHARSIZE: The character size of the color bar annotations. Default is 1.0.
;
; COLOR: The color index of the bar outline and characters. Default
; is !P.Color..
;
; DIVISIONS: The number of divisions to divide the bar into. There will
; be (divisions + 1) annotations. The default is 6.
;
; FONT: Sets the font of the annotation. Hershey: -1, Hardware:0, True-Type: 1.
;
; FORMAT: The format of the bar annotations. Default is '(I5)'.
;
; INVERTCOLORS: Setting this keyword inverts the colors in the color bar.
;
; MAXRANGE: The maximum data value for the bar annotation. Default is
; NCOLORS.
;
; MINRANGE: The minimum data value for the bar annotation. Default is 0.
;
; MINOR: The number of minor tick divisions. Default is 2.
;
; NCOLORS: This is the number of colors in the color bar.
;
; POSITION: A four-element array of normalized coordinates in the same
; form as the POSITION keyword on a plot. Default is
; [0.88, 0.10, 0.95, 0.90] for a vertical bar and
; [0.10, 0.88, 0.90, 0.95] for a horizontal bar.
;;
; RANGE: A two-element vector of the form [min, max]. Provides an
; alternative way of setting the MINRANGE and MAXRANGE keywords.
;
; RIGHT: This puts the labels on the right-hand side of a vertical
; color bar. It applies only to vertical color bars.
;
; TICKNAMES: A string array of names or values for the tick marks.
;
; TITLE: This is title for the color bar. The default is to have
; no title.
;
; TOP: This puts the labels on top of the bar rather than under it.
; The keyword only applies if a horizontal color bar is rendered.
;
; VERTICAL: Setting this keyword give a vertical color bar. The default
; is a horizontal color bar.
;
; COMMON BLOCKS:
;
; None.
;
; SIDE EFFECTS:
;
; Color bar is drawn in the current graphics window.
;
; RESTRICTIONS:
;
; The number of colors available on the display device (not the
; PostScript device) is used unless the NCOLORS keyword is used.
;
; EXAMPLE:
;
; To display a horizontal color bar above a contour plot, type:
;
; LOADCT, 5, NCOLORS=100
; CONTOUR, DIST(31,41), POSITION=[0.15, 0.15, 0.95, 0.75], $
; C_COLORS=INDGEN(25)*4, NLEVELS=25
; COLORBAR, NCOLORS=100, POSITION=[0.15, 0.85, 0.95, 0.90]
;
; MODIFICATION HISTORY:
;
; Written by: David W. Fanning, 10 JUNE 96.
; 10/27/96: Added the ability to send output to PostScript. DWF
; 11/4/96: Substantially rewritten to go to screen or PostScript
; file without having to know much about the PostScript device
; or even what the current graphics device is. DWF
; 1/27/97: Added the RIGHT and TOP keywords. Also modified the
; way the TITLE keyword works. DWF
; 7/15/97: Fixed a problem some machines have with plots that have
; no valid data range in them. DWF
; 12/5/98: Fixed a problem in how the colorbar image is created that
; seemed to tickle a bug in some versions of IDL. DWF.
; 1/12/99: Fixed a problem caused by RSI fixing a bug in IDL 5.2. Sigh... DWF.
; 3/30/99: Modified a few of the defaults. DWF.
; 3/30/99: Used NORMAL rather than DEVICE coords for positioning bar. DWF.
; 3/30/99: Added the RANGE keyword. DWF.
; 3/30/99: Added FONT keyword. DWF
; 5/6/99: Many modifications to defaults. DWF.
; 5/6/99: Removed PSCOLOR keyword. DWF.
; 5/6/99: Improved error handling on position coordinates. DWF.
; 5/6/99. Added MINOR keyword. DWF.
; 5/6/99: Set Device, Decomposed=0 if necessary. DWF.
; 2/9/99: Fixed a problem caused by setting BOTTOM keyword, but not NCOLORS. DWF.
; 8/17/99. Fixed a problem with ambiguous MIN and MINOR keywords. DWF
; 8/25/99. I think I *finally* got the BOTTOM/NCOLORS thing sorted out. :-( DWF.
; 10/10/99. Modified the program so that current plot and map coordinates are
; saved and restored after the colorbar is drawn. DWF.
; 3/18/00. Moved a block of code to prevent a problem with color decomposition. DWF.
; 4/28/00. Made !P.Font default value for FONT keyword. DWF.
; 9/26/00. Made the code more general for scalable pixel devices. DWF.
; 1/16/01. Added INVERTCOLORS keyword. DWF.
; 5/11/04. Added TICKNAME keyword. DWF.
;-
;
;###########################################################################
;
; LICENSE
;
; This software is OSI Certified Open Source Software.
; OSI Certified is a certification mark of the Open Source Initiative.
;
; Copyright � 2000-2004 Fanning Software Consulting.
;
; This software is provided "as-is", without any express or
; implied warranty. In no event will the authors be held liable
; for any damages arising from the use of this software.
;
; Permission is granted to anyone to use this software for any
; purpose, including commercial applications, and to alter it and
; redistribute it freely, subject to the following restrictions:
;
; 1. The origin of this software must not be misrepresented; you must
; not claim you wrote the original software. If you use this software
; in a product, an acknowledgment in the product documentation
; would be appreciated, but is not required.
;
; 2. Altered source versions must be plainly marked as such, and must
; not be misrepresented as being the original software.
;
; 3. This notice may not be removed or altered from any source distribution.
;
; For more information on Open Source Software, visit the Open Source
; web site: http://www.opensource.org.
;
;###########################################################################
PRO COLORBAR, BOTTOM=bottom, CHARSIZE=charsize, COLOR=color, DIVISIONS=divisions, $
FORMAT=format, POSITION=position, MAXRANGE=maxrange, MINRANGE=minrange, NCOLORS=ncolors, $
TITLE=title, VERTICAL=vertical, TOP=top, RIGHT=right, MINOR=minor, $
RANGE=range, FONT=font, TICKLEN=ticklen, _EXTRA=extra, INVERTCOLORS=invertcolors, $
TICKNAMES=ticknames
; Return to caller on error.
On_Error, 2
; Save the current plot state.
bang_p = !P
bang_x = !X
bang_Y = !Y
bang_Z = !Z
bang_Map = !Map
; Are scalable pixels available on the device?
IF (!D.Flags AND 1) NE 0 THEN scalablePixels = 1 ELSE scalablePixels = 0
; Which release of IDL is this?
thisRelease = Float(!Version.Release)
; Check and define keywords.
IF N_ELEMENTS(ncolors) EQ 0 THEN BEGIN
; Most display devices to not use the 256 colors available to
; the PostScript device. This presents a problem when writing
; general-purpose programs that can be output to the display or
; to the PostScript device. This problem is especially bothersome
; if you don't specify the number of colors you are using in the
; program. One way to work around this problem is to make the
; default number of colors the same for the display device and for
; the PostScript device. Then, the colors you see in PostScript are
; identical to the colors you see on your display. Here is one way to
; do it.
IF scalablePixels THEN BEGIN
oldDevice = !D.NAME
; What kind of computer are we using? SET_PLOT to appropriate
; display device.
thisOS = !VERSION.OS_FAMILY
thisOS = STRMID(thisOS, 0, 3)
thisOS = STRUPCASE(thisOS)
CASE thisOS of
'MAC': SET_PLOT, thisOS
'WIN': SET_PLOT, thisOS
ELSE: SET_PLOT, 'X'
ENDCASE
; Here is how many colors we should use.
ncolors = !D.TABLE_SIZE
SET_PLOT, oldDevice
ENDIF ELSE ncolors = !D.TABLE_SIZE
ENDIF
IF N_ELEMENTS(bottom) EQ 0 THEN bottom = 0B
IF N_ELEMENTS(charsize) EQ 0 THEN charsize = 1.0
IF N_ELEMENTS(format) EQ 0 THEN format = '(I5)'
IF N_ELEMENTS(color) EQ 0 THEN color = !P.Color
IF N_ELEMENTS(minrange) EQ 0 THEN minrange = 0
IF N_ELEMENTS(maxrange) EQ 0 THEN maxrange = ncolors
IF N_ELEMENTS(ticklen) EQ 0 THEN ticklen = 0.2
IF N_ELEMENTS(minor) EQ 0 THEN minor = 2
IF N_ELEMENTS(range) NE 0 THEN BEGIN
minrange = range[0]
maxrange = range[1]
ENDIF
IF N_ELEMENTS(divisions) EQ 0 THEN divisions = 6
IF N_ELEMENTS(font) EQ 0 THEN font = !P.Font
IF N_ELEMENTS(title) EQ 0 THEN title = ''
; You can't have a format set *and* use ticknames.
IF N_ELEMENTS(ticknames) NE 0 THEN format = ""
IF KEYWORD_SET(vertical) THEN BEGIN
bar = REPLICATE(1B,20) # BINDGEN(ncolors)
IF Keyword_Set(invertcolors) THEN bar = Reverse(bar, 2)
IF N_ELEMENTS(position) EQ 0 THEN BEGIN
position = [0.88, 0.1, 0.95, 0.9]
ENDIF ELSE BEGIN
IF position[2]-position[0] GT position[3]-position[1] THEN BEGIN
position = [position[1], position[0], position[3], position[2]]
ENDIF
IF position[0] GE position[2] THEN Message, "Position coordinates can't be reconciled."
IF position[1] GE position[3] THEN Message, "Position coordinates can't be reconciled."
ENDELSE
ENDIF ELSE BEGIN
bar = BINDGEN(ncolors) # REPLICATE(1B, 20)
IF Keyword_Set(invertcolors) THEN bar = Reverse(bar, 1)
IF N_ELEMENTS(position) EQ 0 THEN BEGIN
position = [0.1, 0.88, 0.9, 0.95]
ENDIF ELSE BEGIN
IF position[3]-position[1] GT position[2]-position[0] THEN BEGIN
position = [position[1], position[0], position[3], position[2]]
ENDIF
IF position[0] GE position[2] THEN Message, "Position coordinates can't be reconciled."
IF position[1] GE position[3] THEN Message, "Position coordinates can't be reconciled."
ENDELSE
ENDELSE
; Scale the color bar.
bar = BYTSCL(bar, TOP=(ncolors-1 < (255-bottom))) + bottom
; Get starting locations in NORMAL coordinates.
xstart = position(0)
ystart = position(1)
; Get the size of the bar in NORMAL coordinates.
xsize = (position(2) - position(0))
ysize = (position(3) - position(1))
; Display the color bar in the window. Sizing is
; different for PostScript and regular display.
IF scalablePixels THEN BEGIN
TV, bar, xstart, ystart, XSIZE=xsize, YSIZE=ysize, /Normal
ENDIF ELSE BEGIN
bar = CONGRID(bar, CEIL(xsize*!D.X_VSize), CEIL(ysize*!D.Y_VSize), /INTERP)
; Decomposed color off if device supports it.
CASE StrUpCase(!D.NAME) OF
'X': BEGIN
IF thisRelease GE 5.2 THEN Device, Get_Decomposed=thisDecomposed
Device, Decomposed=0
ENDCASE
'WIN': BEGIN
IF thisRelease GE 5.2 THEN Device, Get_Decomposed=thisDecomposed
Device, Decomposed=0
ENDCASE
'MAC': BEGIN
IF thisRelease GE 5.2 THEN Device, Get_Decomposed=thisDecomposed
Device, Decomposed=0
ENDCASE
ELSE:
ENDCASE
TV, bar, xstart, ystart, /Normal
; Restore Decomposed state if necessary.
CASE StrUpCase(!D.NAME) OF
'X': BEGIN
IF thisRelease GE 5.2 THEN Device, Decomposed=thisDecomposed
ENDCASE
'WIN': BEGIN
IF thisRelease GE 5.2 THEN Device, Decomposed=thisDecomposed
ENDCASE
'MAC': BEGIN
IF thisRelease GE 5.2 THEN Device, Decomposed=thisDecomposed
ENDCASE
ELSE:
ENDCASE
ENDELSE
; Annotate the color bar.
IF KEYWORD_SET(vertical) THEN BEGIN
IF KEYWORD_SET(right) THEN BEGIN
PLOT, [minrange,maxrange], [minrange,maxrange], /NODATA, XTICKS=1, $
YTICKS=divisions, XSTYLE=1, YSTYLE=9, $
POSITION=position, COLOR=color, CHARSIZE=charsize, /NOERASE, $
YTICKFORMAT='(A1)', XTICKFORMAT='(A1)', YTICKLEN=ticklen , $
YRANGE=[minrange, maxrange], FONT=font, _EXTRA=extra, YMINOR=minor
AXIS, YAXIS=1, YRANGE=[minrange, maxrange], YTICKFORMAT=format, YTICKS=divisions, $
YTICKLEN=ticklen, YSTYLE=1, COLOR=color, CHARSIZE=charsize, $
FONT=font, YTITLE=title, _EXTRA=extra, YMINOR=minor, YTICKNAME=ticknames
ENDIF ELSE BEGIN
PLOT, [minrange,maxrange], [minrange,maxrange], /NODATA, XTICKS=1, $
YTICKS=divisions, XSTYLE=1, YSTYLE=9, YMINOR=minor, $
POSITION=position, COLOR=color, CHARSIZE=charsize, /NOERASE, $
YTICKFORMAT=format, XTICKFORMAT='(A1)', YTICKLEN=ticklen , $
YRANGE=[minrange, maxrange], FONT=font, YTITLE=title, _EXTRA=extra, $
YTICKNAME=ticknames
AXIS, YAXIS=1, YRANGE=[minrange, maxrange], YTICKFORMAT='(A1)', YTICKS=divisions, $
YTICKLEN=ticklen, YSTYLE=1, COLOR=color, CHARSIZE=charsize, $
FONT=font, _EXTRA=extra, YMINOR=minor
ENDELSE
ENDIF ELSE BEGIN
IF KEYWORD_SET(top) THEN BEGIN
PLOT, [minrange,maxrange], [minrange,maxrange], /NODATA, XTICKS=divisions, $
YTICKS=1, XSTYLE=9, YSTYLE=1, $
POSITION=position, COLOR=color, CHARSIZE=charsize, /NOERASE, $
YTICKFORMAT='(A1)', XTICKFORMAT='(A1)', XTICKLEN=ticklen, $
XRANGE=[minrange, maxrange], FONT=font, _EXTRA=extra, XMINOR=minor
AXIS, XTICKS=divisions, XSTYLE=1, COLOR=color, CHARSIZE=charsize, $
XTICKFORMAT=format, XTICKLEN=ticklen, XRANGE=[minrange, maxrange], XAXIS=1, $
FONT=font, XTITLE=title, _EXTRA=extra, XCHARSIZE=charsize, XMINOR=minor, $
XTICKNAME=ticknames
ENDIF ELSE BEGIN
PLOT, [minrange,maxrange], [minrange,maxrange], /NODATA, XTICKS=divisions, $
YTICKS=1, XSTYLE=1, YSTYLE=1, TITLE=title, $
POSITION=position, COLOR=color, CHARSIZE=charsize, /NOERASE, $
YTICKFORMAT='(A1)', XTICKFORMAT=format, XTICKLEN=ticklen, $
XRANGE=[minrange, maxrange], FONT=font, XMinor=minor, _EXTRA=extra, $
XTICKNAME=ticknames
ENDELSE
ENDELSE
; Restore the previous plot and map system variables.
!P = bang_p
!X = bang_x
!Y = bang_y
!Z = bang_z
!Map = bang_map
END
[전체]
참고 문헌
[논문]
- 없음
[보고서]
- 없음
[URL]
- 없음
문의사항
[기상학/프로그래밍 언어]
- sangho.lee.1990@gmail.com
[해양학/천문학/빅데이터]
- saimang0804@gmail.com
최근댓글