hdfdatas_452.pro

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

;
;
;
;非自动生成功能部分，后添加功能
;
;
;
;
;----------------------------------
;关于HDF及其科学数据集的相关函数
;-----------------------
;打开HDF
pro OpenHDF_452, inFileName ,otHdfID, CREATE = CREATE, WRITE = WRITE, READ = READ
	IF File_Test(inFileName) EQ 0 THEN BEGIN ;如果文件不存在的话，即当输出产品HDF时调用
		otHdfID = HDF_OPEN(inFileName, /CREATE) ;根据文件名创建一个新的HDF文件
	ENDIF ELSE BEGIN
		otHdfID = HDF_OPEN(inFileName, CREATE = CREATE, WRITE = WRITE, READ = READ) ;根据输入的1KM文件的文件名打开其对应的HDF文件
	ENDELSE
end
;-----------------------
;开始SD,inType默认为只读方式
pro StartSD_452, inFileName, otSDInterfaceID, inType = inType ;函数的第三个为关键字，inType为打开科学数据集的方式，是只读或者是可读可写之类
	IF N_ELEMENTS(inType) EQ 0 THEN BEGIN
		otSDInterfaceID = HDF_SD_START(inFileName, /READ) ;如果没有指定inType的话则以只读方式打开科学数据集
		RETURN
	ENDIF

	IF (STRMATCH(inType, 'RDWR', /FOLD_CASE) EQ 1) THEN BEGIN ;如果指定了inType，则判读是以什么方式开始科学数据集
		otSDInterfaceID = HDF_SD_START(inFileName, /RDWR) ;以可读可写方式开始科学数据集
	ENDIF ELSE IF (STRMATCH(inType, 'CREATE', /FOLD_CASE) EQ 1) THEN BEGIN
		otSDInterfaceID = HDF_SD_START(inFileName, /CREATE) ;以重新创建一个新的HDF文件的方式开始科学数据集
	ENDIF ELSE BEGIN
		otSDInterfaceID = HDF_SD_START(inFileName, /READ) ;以只读方式开始科学数据集
	ENDELSE
end
;-----------------------
;根据属性名得到属性值
pro GetAttributeValue_452, inSDIcOrSDId, inAttrName, otAttrValues
	IF inSDIcOrSDId EQ -1 THEN RETURN

	;根据属性的名称得到属性的ID号
	attrIndex = HDF_SD_ATTRFIND(inSDIcOrSDId, inAttrName) & IF attrIndex EQ -1 THEN RETURN

	;读取属性信息
	HDF_SD_ATTRINFO, inSDIcOrSDId, attrIndex, DATA = Data, NAME = name,TYPE = type, COUNT = count;返回指定属性当中的值都返回到data当中，值的名称返回到name当中，值的类型返回到type当中，值的个数返回到count当中

	;返回属性的值
	otAttrValues = {DATA:Data, NAME:name, TYPE:type, COUNT:count}
end
;------------------------
;选择SD科学数据集
pro SelectSD_452, inSDInterfaceID, inSDName, otSDIndex, otSDId
	otSDIndex = HDF_SD_NAMETOINDEX(inSDInterfaceID, inSDName);Returns the specified SD dataset index number
	IF N_ELEMENTS(otSDIndex) EQ 0 THEN RETURN
	IF otSDIndex EQ -1 THEN RETURN
	otSDId = HDF_SD_SELECT(inSDInterfaceID, otSDIndex);Returns the specified SD dataset's identifier
end
;------------------------
;读取SD一层的数据(DN值)
PRO GetSDOneLayerData_452, inSDId, inLayerNO, otDatas
	HDF_SD_GETINFO, inSDId, DIMS = otDimens;返回'EV_1KM_RefSB'的维数，其维数是以颠倒的顺序返回，除非指定了NOREVERSE关键字
	IF N_ELEMENTS(otDimens) LE 2 THEN BEGIN ;otDimens最多为3唯，且其顺序是颠倒的所以otDimens[0]为列数，otDimens[1]为行数，otDimens[2]为层数
		HDF_SD_GETDATA, inSDId, otDatas, START = [0, 0], STRIDE = [0, 0];只有一层即layer为1
	ENDIF ELSE BEGIN
		dimLayer = otDimens[2];layer大于1
		HDF_SD_GETDATA, inSDId, otDatas, START = [0, 0, inLayerNO], STRIDE = [0, 0, dimLayer];则全部layer里的数据都读
	ENDELSE
END
;------------------------
;不选择SD
pro UnSelectSD_452, inSDID
	HDF_SD_ENDACCESS, inSDID ;不选择科学数据集，在不选择科学数据集之后才能关闭科学数据集
end
;------------------------
;结束SD
pro EndSD_452, inSDInterfaceID
	HDF_SD_END, inSDInterfaceID ;关闭科学数据集
end
;------------------------
;关闭HDF
pro CloseHDF_452, inHdfID
	IF inHdfID EQ -1 THEN RETURN
	HDF_CLOSE, inHdfID ;关闭HDF
end
;------------------------
;把DN值转换为RR值
pro DNsToRRs_452, inDNDatas, inScale, inOffset, otRRDatas
	otRRDatas = inScale * (Temporary(inDNDatas) - inOffset) ;Temporary用于创建一个临时空间来保存inDNDatas，以节省内存空间
end
;--------------------------
;创建产品标签
pro CreateEnvRSFlag_452, inFileName, inTypeDesc, inTypeCode
	AddAttributeToFile_452, inFileName, 'ProductDescription', '陆表植被应用分系统反演产品@中国科学院遥感应用研究所国家航天局航天遥感论证中心', /STRING
	AddAttributeToFile_452, inFileName, 'TypeDescription', inTypeDesc, /STRING
	AddAttributeToFile_452, inFileName, 'TypeCode', inTypeCode, /INT
end
;--------------------------
;把属性添加到文件中
pro AddAttributeToFile_452, inFileName, inAttrName, inAttrValue, inSDName = inSDName,$
						BYTE = BYTE, DFNT_CHAR = DFNT_CHAR, SHORT= SHORT, STRING = STRING , $
						DFNT_FLOAT32 = DFNT_FLOAT32, DFNT_FLOAT64 = DFNT_FLOAT64, DFNT_INT8 = DFNT_INT8, $
		 				DFNT_INT16 = DFNT_INT16, DFNT_INT32 = DFNT_INT32, DFNT_UINT8 = DFNT_UINT8,$
		 				DFNT_UINT16 = DFNT_UINT16, DFNT_UINT32 = DFNT_UINT32, DOUBLE = DOUBLE, $
						FLOAT = FLOAT, INT = INT, LONG = LONG
	bFileExist = File_Test(inFileName)
	IF bFileExist EQ 1 THEN OpenHDF_452, inFileName, otHdfID, /WRITE  ;如果文件存在则覆盖文件
	IF bFileExist EQ 0 THEN OpenHDF_452, inFileName, otHdfID, /CREATE ;如果文件不存在则重新创建一个新的HDF文件
	StartSD_452, inFileName, otSDInterfaceID, inType = 'RDWR'

	IF N_ELEMENTS(inSDName) EQ 0 THEN BEGIN
		HDF_SD_ATTRSET, otSDInterfaceID, inAttrName, inAttrValue, $
						BYTE = BYTE, DFNT_CHAR = DFNT_CHAR, SHORT= SHORT, STRING = STRING , $
						DFNT_FLOAT32 = DFNT_FLOAT32, DFNT_FLOAT64 = DFNT_FLOAT64, DFNT_INT8 = DFNT_INT8, $
		 				DFNT_INT16 = DFNT_INT16, DFNT_INT32 = DFNT_INT32, DFNT_UINT8 = DFNT_UINT8,$
		 				DFNT_UINT16 = DFNT_UINT16, DFNT_UINT32 = DFNT_UINT32, DOUBLE = DOUBLE, $
						FLOAT = FLOAT, INT = INT, LONG = LONG
	ENDIF ELSE BEGIN
		SelectSD_452, otSDInterfaceID, inSDName, otSDIndex, otSDId
		HDF_SD_ATTRSET, otSDId, inAttrName, inAttrValue, $
						BYTE = BYTE, DFNT_CHAR = DFNT_CHAR, SHORT= SHORT, STRING = STRING , $
						DFNT_FLOAT32 = DFNT_FLOAT32, DFNT_FLOAT64 = DFNT_FLOAT64, DFNT_INT8 = DFNT_INT8, $
		 				DFNT_INT16 = DFNT_INT16, DFNT_INT32 = DFNT_INT32, DFNT_UINT8 = DFNT_UINT8,$
		 				DFNT_UINT16 = DFNT_UINT16, DFNT_UINT32 = DFNT_UINT32, DOUBLE = DOUBLE, $
						FLOAT = FLOAT, INT = INT, LONG = LONG
		UnSelectSD_452, otSDId
	ENDELSE

	EndSD_452, otSDInterfaceID
	CloseHDF_452, otHdfID
end
;------------------------------
pro GetCoordEnvelope_612, inFileName, otEnvelope
    OpenHDF_452, inFileName, otHdfID, /READ
	StartSD_452, inFileName, otSDInterfaceID
	GetCoordEnvelopeFromHDF_322, otSDInterfaceID, otEnvelope
	EndSD_452, otSDInterfaceID
	CloseHDF_452, otHdfID
end
;------------------------------
;获取外接矩形经纬度范围
PRO GetCoordEnvelopeFromHDF_322, inSDInterfaceID, otEnvelope
	strArchive = 'ArchiveMetadata.0'						;----李家国修改12.08，获取投影后的标签范围
	GetAttributeValue_452, inSDInterfaceID, 'OldArchiveMetadata.0', otAttrValues
	IF N_ELEMENTS(otAttrValues) NE 0 THEN strArchive = 'OldArchiveMetadata.0'

	AttrName = ['NORTHBOUNDINGCOORDINATE', 'SOUTHBOUNDINGCOORDINATE', 'EASTBOUNDINGCOORDINATE', $
	'WESTBOUNDINGCOORDINATE', strArchive[0]]
	n = -1
	for i=0, 4 do begin
		attrIndex = HDF_SD_ATTRFIND(inSDInterfaceID, AttrName[i])
		if attrIndex ne -1 then n = n+i
	endfor
	otEnvelope = FltArr(4)
	if (n ne 3) && (n ne 5) then begin
		otEnvelope = [0.0, 0.0, 0.0, 0.0]
		RETURN
	endif
	if n eq 5 then begin
		for ii=0, 3 do begin
			GetAttributeValue_452, inSDInterfaceID, AttrName[ii], otAttrValues
			otEnvelope[ii] = otAttrValues.DATA
		endfor
	endif
	if n eq 3 then begin
		GetAttributeValue_452, inSDInterfaceID, AttrName[4], otAttrValues
		CoreMD = otAttrValues.DATA
		for jj=0, 3 do begin
			Pos1 = STRPOS(CoreMD, AttrName[jj])
			Pos2 = STRPOS(CoreMD, AttrName[jj], /REVERSE_SEARCH)
			Str1 = STRMID(CoreMD, Pos1+1, Pos2-Pos1-1)

			p1 = STRPOS(Str1, 'VALUE')
			p2 = STRPOS(Str1, 'END_OBJECT')
			Str2 = STRMID(Str1, p1+1, p2-p1-1)
			p3 = STRPOS(Str2, '=')
			Str3 = STRMID(Str2, p3+1)
			otEnvelope[jj] = Float(StrTrim(Str3, 2))
		endfor
	endif
END
;-------------------------------
;;得到Reflectance（反射率）值
pro GetReflectance_610, inFileName, inSDName, inBandName, otReflectanceData, otValidData
	;打开HDF，打开科学数据集（SD），选中SD
	OpenHDF_452, inFileName, otHdfID

	StartSD_452, inFileName, otSDInterfaceID

	SelectSD_452, otSDInterfaceID, inSDName, otSDIndex, otSDId


    ;;从科学数据集中获取数据
	;根据属性名得到属性值,band_names是某个数据集里面的一个属性文件
	GetAttributeValue_452, otSDId, 'band_names', otBandValus
	bandNames = otBandValus.data & bandNames = STRSPLIT(bandNames, ',', ESCAPE=' ', /EXTRACT)

	;判断输入的波段名称与数据集里面的波段数中处于第几个（如1公里的有：21~36（除了26））
	FOR ii = 0, N_ELEMENTS(bandNames) DO BEGIN
		IF inBandName EQ bandNames[ii] THEN BEGIN
			iBandIndex = ii
			BREAK
		ENDIF ELSE iBandIndex = -1
	ENDFOR
	IF iBandIndex EQ -1 THEN RETURN
	;根据波段索引值，读取SD一层的数据(DN值)，其中otLayerData储存了dn值的数据
	GetSDOneLayerData_452, otSDId, iBandIndex, otLayerData

	;判断DN值是否有效，并记录到矩阵otValidData
    otValidData = otLayerData le 32767

	;获取属性信息
	strAttrName = ['reflectance_scales', 'reflectance_offsets']
	GetAttributeValue_452, otSDId, strAttrName[0], otReflectanceScaleValues
	GetAttributeValue_452, otSDId, strAttrName[1], otReflectanceOffsetValues

	;把DN值转换为反射率的值
	scaleValue = otReflectanceScaleValues.Data[iBandIndex]
	offsetValue = otReflectanceOffsetValues.Data[iBandIndex]
	DNsToRRs_452, otLayerData, scaleValue, offsetValue, otReflectanceData

	;关闭SD，关闭HDF
	UnSelectSD_452, otSDId
	EndSD_452, otSDInterfaceID
	CloseHDF_452, otHdfID
end