ndvimulti_803.bak

利用IDL读取MODIS影像L1B和MOD03数据(HDF格式),计算遥感植被指数NDVI

;
;
;
;非自动生成功能部分，后添加功能
;
;
;
;
;----------------------------------
;归一化植被指数NDVI产品
;-----------------------
;多文件的归一化植被指数NDVI计算
pro NDVIMulti_803, infilenames,  inSavePath = inSavePath, otSaveNames
    ;归一化植被指数NDVI计算
    iFileCount = N_ELEMENTS(inFileNames) & print,ifilecount & IF iFileCount EQ 0 THEN RETURN
	otSaveNames = STRARR(iFileCount) ;根据输入的MODIS02的文件数量，创建一个同样元素数量的数组，用于存放产品保存名称

    ;设置保存路径及保存文件名
	for ii = 0 , iFileCount - 1 do begin
		if N_ELEMENTS(inSavePath) eq 0 then begin ;如果使用默认保存路径
			strSaveHDFName = StrMid( inFileNames[ii], 0, StrLen(inFileNames[ii]) - 4)+'_NDVI.hdf' ;在hdf文件后面加NDVI保存
		endif else begin  ;使用指定路径进行保存
			strSaveHDFName = inSavePath + File_BaseName( inFileNames[ii], '.hdf', /FOLD_CASE) + '_NDVI.hdf'  ;File_BaseName能将文件名前面的路径去掉
		endelse

		otSaveNames[ii] = strSaveHDFName

		;处理每一组数据
		HDFNDVI_803, inFileNames[ii],  strSaveHDFName[0];计算各个文件的归一化植被指数NDVI
	endfor

end
;-----------------------
pro HDFNDVI_803, inFileName, inOtFileName

    ;初始化进度条
	strReportInit = ['输入文件：', inFileName, $
	'输出文件：', inOtFileName] ;处理过程当中显示的内容
	;Envi_Report_Init是envi里面的一个函数，用于显示工作进度的
	Envi_Report_Init, strReportInit, Title = '归一化植被指数NDVI计算进度', Base = reportBaseID, /INTERRUPT ;初始化Envi进程报告
	WIDGET_CONTROL, reportBaseID[3], Set_Value = '取消'
	iDen =4 & iNum = 0
	;---------------------------
	;读取MOD02数据
	Envi_Report_Stat, reportBaseID, ++iNum, iDen, CANCEL = cancelValue
	IF cancelValue EQ 1 THEN BEGIN  ;当按下cancel是取消进度条
		Envi_Report_Init,' ', Base = reportBaseID, /FINISH & RETURN ;点击进度窗体上的取消按钮时，则将进度报告窗体清空再返回
	ENDIF
    GetCoordEnvelope_612, inFileName, Envelope ;读取MODIS02文件，获取东西南北四个角点的坐标
    ;---------------------------
    ;创建新产品的HDF文件
    ENVI_REPORT_STAT, reportBaseID, ++iNum, iDen, CANCEL = cancelValue
	IF cancelValue EQ 1 THEN BEGIN
		Envi_Report_Init, ' ', Base = reportBaseID, /FINISH & RETURN
	ENDIF
	OpenHDF_452, inOtFileName, otHdfID, /CREATE, /WRITE  ;用可读可写的方式打开HDF
	CreateEnvRSFlag_452, inOtFileName, '归一化植被指数NDVI', 202;创建产品标签
	StartSD_452, inOtFileName, otSDInterfaceID, inType = 'RDWR'

    ;写入坐标信息

	HDF_SD_ATTRSET, otSDInterfaceID, 'MOD02 input file',  File_BaseName(inFileName) , /STRING
	HDF_SD_ATTRSET, otSDInterfaceID, 'NORTHBOUNDINGCOORDINATE', Envelope[0], /DFNT_FLOAT64
	HDF_SD_ATTRSET, otSDInterfaceID, 'SOUTHBOUNDINGCOORDINATE', Envelope[1], /DFNT_FLOAT64
	HDF_SD_ATTRSET, otSDInterfaceID, 'EASTBOUNDINGCOORDINATE', Envelope[2], /DFNT_FLOAT64
	HDF_SD_ATTRSET, otSDInterfaceID, 'WESTBOUNDINGCOORDINATE', Envelope[3], /DFNT_FLOAT64
    ;---------------------------
    ;判断文件类型
    ENVI_REPORT_STAT, reportBaseID, ++iNum, iDen, CANCEL = cancelValue
	IF cancelValue EQ 1 THEN BEGIN
		Envi_Report_Init, ' ', Base = reportBaseID, /FINISH & RETURN
	ENDIF

    ;判断数据是1km、500米、205米
    FileIsMOD021KM_322, inFileName, ot1KMFlag
    FileIsMOD02HKM_322, inFileName, otHKMFlag
    FileIsMOD02QKM_322, inFileName, otQKMFlag


    ;---------------------------
    ;根据MODIS02的不同文件计算归一化植被指数
    ENVI_REPORT_STAT, reportBaseID, ++iNum, iDen, CANCEL = cancelValue
	IF cancelValue EQ 1 THEN BEGIN
		Envi_Report_Init, ' ', Base = reportBaseID, /FINISH & RETURN
	ENDIF
    ;1km数据
    if ot1KMFlag eq 1 then begin
		;1波段为红外波段，2为近红外波段

			GetReflectance_610, inFileName, 'EV_250_Aggr1km_RefSB', 1, otRed, otValidData1
			GetReflectance_610, inFileName, 'EV_250_Aggr1km_RefSB', 2, otNearInfrared, otValidData2
			otNDVI = (otNearInfrared - otRed)/(otNearInfrared + otRed)
			index = where(otValidData1 eq 0)
			otNDVI[index] =  255
	        CreateDSAndAddData_324, otSDInterfaceID, 'NDVI', '归一化植被指数', otNDVI, /Float
		    SelectSD_452, otSDInterfaceID, 'NDVI', otSDIndex, otSDId
		    HDF_SD_ATTRSET, otSDId, '_FillValue', 255, /DFNT_UINT16
		    HDF_SD_ATTRSET, otSDId, 'ScaleValue', 100, /Float
		    HDF_SD_ATTRSET, otSDId, 'offset', 0, /FLOAT
		    UnSelectSD_452, otSDId

    endif

    ;500m数据
    if otHKMFlag eq 1 then begin
		;1波段为红外波段，2为近红外波段

			GetReflectance_610, inFileName, 'EV_250_Aggr500_RefSB', 1, otRed, otValidData1
			GetReflectance_610, inFileName, 'EV_250_Aggr500_RefSB', 2, otNearInfrared, otValidData2
			otNDVI = (otNearInfrared - otRed)/(otNearInfrared + otRed)
			index = where(otValidData1 eq 0)
			otNDVI[index] =  255
	        CreateDSAndAddData_324, otSDInterfaceID, 'NDVI', '归一化植被指数', otNDVI, /Float
		    SelectSD_452, otSDInterfaceID, 'NDVI', otSDIndex, otSDId
		    HDF_SD_ATTRSET, otSDId, '_FillValue', 255, /DFNT_UINT16
		    HDF_SD_ATTRSET, otSDId, 'ScaleValue', 100, /Float
		    HDF_SD_ATTRSET, otSDId, 'offset', 0, /FLOAT
		    UnSelectSD_452, otSDId
    endif

    ;250米数据
    if otQKMFlag eq 1 then begin
		;1波段为红外波段，2为近红外波段

			GetReflectance_610, inFileName, 'EV_250_RefSB', 1, otRed, otValidData1
			GetReflectance_610, inFileName, 'EV_250_RefSB', 2, otNearInfrared, otValidData2
			otNDVI = (otNearInfrared - otRed)/(otNearInfrared + otRed)
			index = where(otValidData1 eq 0)
			otNDVI[index] =  255
	        CreateDSAndAddData_324, otSDInterfaceID, 'NDVI', '归一化植被指数', otNDVI, /Float
		    SelectSD_452, otSDInterfaceID, 'NDVI', otSDIndex, otSDId
		    HDF_SD_ATTRSET, otSDId, '_FillValue', 255, /DFNT_UINT16
		    HDF_SD_ATTRSET, otSDId, 'ScaleValue', 100, /Float
		    HDF_SD_ATTRSET, otSDId, 'offset', 0, /FLOAT
		    UnSelectSD_452, otSDId

    endif

    EndSD_452, otSDInterfaceID
	CloseHDF_452, otHdfID

    ;finish reportBar
    Envi_Report_Init, ' ', Base = reportBaseID, /FINISH
end