|
;+
; :Description:
; 自动统计置信区间（如5%）对应的DN值
; 可用于百分比线性拉伸、植被覆盖度NDVImin和NDVImax的计算等
; 注1：用于调用了ENVIRasterStatistics，此代码最低支持ENVI52/IDL84版本
; 注2：本程序实现初衷是为了计算植被覆盖度，如果仅仅是为了图像拉伸，
; 可以调用一系列虚拟栅格或ENVITask，如 ENVILinearPercentStretchRaster
;
; :Author: duhj@geoscene.cn
;
; :Date: 2016-8-10 11:48:30
;-
;拉伸函数，输入raster只能是单波段或三波段栅格
;percent输入百分比，默认为0.02
;返回值为拉伸后数组
FUNCTION LinearPercentStretch, raster, percent=percent
COMPILE_OPT idl2
;percent默认为0.02
IF ~N_ELEMENTS(percent) THEN percent=0.02
;获取波段数，如果不为1或3，则返回空
nb = raster.NBANDS
IF nb NE 1 AND nb NE 3 THEN RETURN, !NULL
;初始化字节型数组，用于存储拉伸后结果
result = BYTARR(raster.NCOLUMNS, raster.NROWS, nb)
;循环获取每个波段拉伸结果，存储在result中
FOR i=0,nb-1 DO BEGIN
subRaster = ENVISubsetRaster(raster, BANDS=[i])
;根据置信区间统计 minValue和maxValue，
;可用于植被覆盖度计算（对应NDVImin和NDVImax）
;当统计结果不理想时，说明异常值较多，
;可以设置HISTOGRAM_NBINS=1000，或更大的值
;这里使用了ENVI接口进行直方图统计，
;也可以换成IDL的HISTOGRAM函数，用法查看帮助
stats = ENVIRasterStatistics(subRaster, /HISTOGRAMS, $
HISTOGRAM_NBINS=256)
;统计直方图，获取累计像元个数
hist = stats.HISTOGRAMS[0]
totalCounts = TOTAL(hist['counts'], /CUMULATIVE)
;获取累积百分比
totalPercents = totalCounts/totalCounts[-1]
;获取置信区间的索引
tmp = MIN(totalPercents-percent, beginIndex, /ABSOLUTE)
tmp = MIN(totalPercents-(1.0-percent), endIndex, /ABSOLUTE)
;根据索引计算置信区间最小和最大值
minValue = hist['min']+beginIndex*hist['binsize']
maxValue = hist['min']+endIndex*hist['binsize']
data = subRaster.GetData()
data = BYTSCL(data, MIN=minValue, MAX=maxValue)
result[*,*,i] = data
ENDFOR
;返回结果
RETURN, result
END
;主程序
PRO test_Confidence
COMPILE_OPT idl2
e=envi(/headless)
;使用ENVI自带数据
file = FILEPATH('qb_boulder_msi', ROOT_DIR=e.ROOT_DIR, $
SUBDIRECTORY=['data'])
raster = e.OpenRaster(file)
;获取RGB真彩色栅格
rgbRaster = ENVISubsetRaster(raster, BANDS=[2,1,0])
;获取2%线性拉伸结果
result = LinearPercentStretch(rgbRaster, percent=0.02)
;显示使用BYTSCL拉伸效果（即线性拉伸）
i = IMAGE(BYTSCL(rgbRaster.GetData()), LAYOUT=[2,1,1], $
DIMENSIONS=[800,400],/ORDER, TITLE='线性拉伸结果')
;获取百分比线性拉伸结果（如2%线性拉伸）
i = IMAGE(result, LAYOUT=[2,1,2] ,/ORDER, $
TITLE='2%线性拉伸结果', /CURRENT)
END
|
