这是一个用于计算生长季始期和生长季末期的IDL程序，根据的算法是《近10 年青藏高原高寒草地物候时空变化特征分析》，需要指出的是，这里没有对NDVI数据进行滤波处理。
; Author: zhb
; E-mail: zhanghongbo@itpcas.ac.cn
; Description: compute the threshold of SOG and EOG.
; Date: 2013.10.10
PRO Compute_SOG
; 填写输入参数
; FileName 是原始文件，即包含有所有年份NDVI数据的文件，波段数应是23的倍数
FileName = ‘H:\yu for test data\SGVO_SGcalibrate_snow_namco_423.tif’ ; 完整的路径名
WorkDir = ‘H:\yu for test data\’ ; 工作目录/输出目录，别忘了路径的最后要有一个斜杠
ComputeMean = 1 ; 0 表示不需要计算均值， 1 表示相反
; 开始工作
print, ‘Start using ENVI routines, firstly get file info’
compile_opt IDL2 ; force idl to recognize vars only by square brackets
ENVI, /RESTORE_BASE_SAVE_FILES ; restore some needed savs
ENVI_BATCH_INIT ; start a new ENVI session with no GUI
ENVI_OPEN_FILE, FileName, R_FID = fid, /NO_REALIZE ; get fid
ENVI_FILE_QUERY, fid, NS = cols, NL = rows, DATA_IGNORE_VALUE = nodata, $
NB = bands, DATA_TYPE = type, FNAME = fname, DIMS = dims ; get some needed info
print, “Current data file: “, fname
print, FORMAT = ‘(“FID: “, I1, 2X, “rows: “, I4, 2X, “cols: “, I4, 2X, “bands: “, I3, 2X, “type: “, I2)’ $
, fid, rows, cols, bands, type
print, FORMAT = ‘(“Ignored Data: “, F)’, nodata
print, “DIMS: “, dims
mapinfo = ENVI_GET_MAP_INFO(FID = fid)
if ComputeMean then begin
; if ComputeMean is not set, we needn’t get its data
print, ‘Create a suitable data array’
;create data array
case type of
1: data = bytarr(cols, rows, bands) ; Bytes
2: data = intarr(cols, rows, bands) ; 16-bits integer
3: data = lonarr(cols, rows, bands) ; 32-bits integer
4: data = fltarr(cols, rows, bands) ; float
5: data = dblarr(cols, rows, bands) ; double
12: data = uintarr(cols, rows, bands) ; unsigned 16-bits integer
13: data = ulonarr(cols, rows, bands) ; unsigned 32-bits integer
endcase
print, ‘Read all data into this array’
; read all data
for i = 0, bands-1 do begin
iband = ENVI_GET_DATA(FID = fid, DIMS = dims, POS = i)
data[*,*,i] = iband
endfor
yrecs = 23 ; records per year
; get num of years
ynum = bands / yrecs
print, ‘Start computing each image”s mean value among all years if ”ComputeMean” is set’
for i = 0, yrecs-1 do begin
; create ave buffer and num array
num = make_array(cols, rows, /int, VALUE = 0)
sum = make_array(cols, rows, /float, VALUE = -9999)
for j = 0, ynum-1 do begin
tmp = data[*,*,j*yrecs + i]
sub = where(tmp ne nodata and sum ne -9999, /NULL)
if sub ne !NULL then sum[sub] += tmp[sub] ; add to total
delvar, sub
sub = where(tmp ne nodata and sum eq -9999, /NULL)
if sub ne !NULL then sum[sub] = tmp[sub] ; initialize unexlored situs
delvar, sub
num[where(tmp ne nodata, /NULL)]++
delvar, tmp
endfor
ave = make_array(cols, rows, /float, VALUE = nodata)
sub = where(num ne 0, /NULL)
ave[sub] = sum[sub] / num[sub]
outname = string(FORMAT = ‘(A, A,I02)’, WorkDir, ‘num’, i)
ENVI_WRITE_ENVI_FILE, ave, OUT_NAME = outname, $
MAP_INFO = mapinfo, /NO_OPEN
print, outname, ‘ has been output.’
delvar, num, sum, ave, sub
endfor
delvar, data
endif
; close file connection
ENVI_FILE_MNG, ID = fid, /REMOVE
print, ‘Start computing SOG-date and EOG-date’
; compute what we really want
diff0_p = make_array(cols, rows, /float, VALUE = -9999) ; positive difference, receiving last round’s values
diff0_n = make_array(cols, rows, /float, VALUE = -9999) ; negative difference, receiving last round’s values
diff1_p = make_array(cols, rows, /float, VALUE = -9999) ; positive difference, receiving current round’s values
diff1_n = make_array(cols, rows, /float, VALUE = -9999) ; negative difference, receiving current round’s values
flag_p = make_array(cols, rows, /integer, VALUE = -1) ; yrec mark for the max positive diff
flag_n = make_array(cols, rows, /integer, VALUE = -1) ; yrec mark for the max negative diff ( absolute value )
inname = string(FORMAT = ‘(A, A,I02)’, WorkDir, ‘num’, 0)
ENVI_OPEN_FILE, inname, R_FID = fid, /NO_REALIZE
; now, we don’t check the info of input file so you have to ensure it
data0 = ENVI_GET_DATA(FID = fid, DIMS = dims, POS = 0)
ENVI_FILE_MNG, ID = fid, /REMOVE
inname = string(FORMAT = ‘(A, A,I02)’, WorkDir, ‘num’, 1)
ENVI_OPEN_FILE, inname, R_FID = fid, /NO_REALIZE
data1 = ENVI_GET_DATA(FID = fid, DIMS = dims, POS = 0)
ENVI_FILE_MNG, ID = fid, /REMOVE
; here, we suppose the nodata value equals to the original one, so you have to ensure it
sub = where(data0 ne nodata and data1 ne nodata and data1 gt data0, /NULL)
if sub ne !NULL then diff0_p[sub] = data1[sub] – data0[sub]
flag_p[sub] = 1 ; 1 means SOG ocurred at the 1st 16days, currently
sub = where(data0 ne nodata and data1 ne nodata and data1 lt data0, /NULL)
if sub ne !NULL then diff0_n[sub] = data0[sub] – data1[sub]
flag_n[sub] = 1 ; 1 means EOG ocurred at the 1st 16days, currently
for i = 2, yrecs-1 do begin
data0 = data1
inname = string(FORMAT = ‘(A,A,I02)’, WorkDir, ‘num’, i)
ENVI_OPEN_FILE, inname, R_FID = fid, /NO_REALIZE
data1 = ENVI_GET_DATA(FID = fid, DIMS = dims, POS = 0)
ENVI_FILE_MNG, ID = fid, /REMOVE
; compute current round’s differece
sub = where(data0 ne nodata and data1 ne nodata and data1 gt data0, /NULL)
if sub ne !NULL then diff1_p[sub] = data1[sub] – data0[sub]
delvar, sub
sub = where(data0 ne nodata and data1 ne nodata and data1 lt data0, /NULL)
if sub ne !NULL then diff1_n[sub] = data0[sub] – data1[sub]
delvar, sub
; compute SOG
sub = where(diff1_p gt diff0_p, /NULL)
if sub ne !NULL then diff0_p[sub] = diff1_p[sub]
flag_p[sub] = i
delvar, sub
; compute EOG
sub = where(diff1_n gt diff0_n, /NULL)
if sub ne !NULL then diff0_n[sub] = diff1_n[sub]
flag_n[sub] = i
delvar, sub
endfor
delvar, data0, data1, diff0_p, diff0_n, diff1_p, diff1_n
print, ‘Start computing unexcepted cases’
; compute false_LOG
len = make_array(cols, rows, /integer, VALUE = -99)
sub = where(flag_n ne -1 and flag_p ne -1, /NULL)
if sub ne !NULL then len[sub] = flag_n[sub] – flag_p[sub]
delvar, sub
; create exception data
ex = make_array(cols, rows, /integer, VALUE = -99)
sub = where(len le 0, /NULL)
if sub ne !NULL then ex[sub] = len[sub]
delvar, sub
; 这里我写的有点繁琐，牺牲了运行速度，但是大大节约了系统内存，虽然其实意义不大
print, ‘Start computing the threshold of SOG for each grid’
sog = make_array(cols, rows, /float, VALUE = -1)
; find needed input file number
valid = flag_p[where(flag_p ne -1)]
nums = valid[uniq(valid, sort(valid))]
print, “needed num: “, nums
delvar, valid
for i = 0, N_ELEMENTS(nums)-1 do begin
num = nums[i]
inname = string(FORMAT = ‘(A,A,I02)’, WorkDir, ‘num’, num)
ENVI_OPEN_FILE, inname, R_FID = fid, /NO_REALIZE
data = ENVI_GET_DATA(FID = fid, DIMS = dims, POS = 0)
ENVI_FILE_MNG, ID = fid, /REMOVE
sub = where(flag_p eq num, /NULL)
if sub ne !NULL then sog[sub] = data[sub]
delvar, data, sub
endfor
print, ‘Start computing the threshold of EOG for each grid’
eog = make_array(cols, rows, /float, VALUE = -1)
; find needed input file number
valid = flag_n[where(flag_n ne -1)]
nums = valid[uniq(valid, sort(valid))]
print, “needed num: “, nums
delvar, valid
for i = 0, N_ELEMENTS(nums)-1 do begin
num = nums[i]
inname = string(FORMAT = ‘(A,A,I02)’, WorkDir, ‘num’, num)
ENVI_OPEN_FILE, inname, R_FID = fid, /NO_REALIZE
data = ENVI_GET_DATA(FID = fid, DIMS = dims, POS = 0)
ENVI_FILE_MNG, ID = fid, /REMOVE
sub = where(flag_n eq num, /NULL)
if sub ne !NULL then eog[sub] = data[sub]
delvar, data, sub
endfor
; let’s output
ENVI_WRITE_ENVI_FILE, flag_p, OUT_NAME = WorkDir + ‘SOG_flag’, $
MAP_INFO = mapinfo, /NO_OPEN
print, WorkDir + ‘SOG_flag’ + ‘ has been output’
ENVI_WRITE_ENVI_FILE, flag_n, OUT_NAME = WorkDir + ‘EOG_flag’, $
MAP_INFO = mapinfo, /NO_OPEN
print, WorkDir + ‘EOG_flag’ + ‘ has been output’
ENVI_WRITE_ENVI_FILE, sog, OUT_NAME = WorkDir + ‘SOG_threshold’, $
MAP_INFO = mapinfo, /NO_OPEN
print, WorkDir + ‘SOG_threshold’ + ‘has been output’
ENVI_WRITE_ENVI_FILE, eog, OUT_NAME = WorkDir + ‘EOG_threshold’, $
MAP_INFO = mapinfo, /NO_OPEN
print, WorkDir + ‘EOG_threshold’ + ‘has been output’
ENVI_WRITE_ENVI_FILE, ex, OUT_NAME = WorkDir + ‘EXCEPTION’, $
MAP_INFO = mapinfo, /NO_OPEN
print, WorkDir + ‘EXCEPTION’ + ‘ has been output’
print, ‘Task is finished’
ENVI_BATCH_EXIT ; terminate current ENVI session
END