结合以上基本信息,我们已经对MCD15A3H产品的基本信息有了一定了解。接下来就可以进行栅格数据的读取与处理、筛选了。
在这里需要注意的是,之前的两篇文章Python GDAL读取栅格数据并基于质量评估波段QA对指定数据加以筛选掩膜以及Python批量读取HDF多波段栅格数据并绘制像元直方图已经对本次所要用到的大部分需求与代码加以实现并进行了详细讲解,这里就不再赘述。本文代码所实现功能与上述第一篇博客中的需求一致,唯一不同的是将GLASS产品更改为了MCD15A3H产品,且仅需对MCD15A3H产品的主算法像元加以做差计算(也就是筛选出MCD15A3H产品中第一个QC波段对应二进制数的第一位为0的像元,其它像元就不用参与差值计算了)。
具体代码如下:
[code]# -*- coding: utf-8 -*-"""Created on Sun Jul 25 14:57:45 2021@author: fkxxgis"""import osimport copyimport numpy as npfrom osgeo import gdalrt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/RT_LAI/"mcd15_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/mcd15A3H/"out_file_path="G:/Postgraduate/LAI_Glass_RTlab/Test_DRT/"rt_file_list=os.listdir(rt_file_path)for rt_file in rt_file_list: rt_file_split=rt_file.split("_") rt_hv=rt_file_split[3][:-4] mcd15_file_list=os.listdir(mcd15_file_path) for mcd15_file in mcd15_file_list: if rt_hv in mcd15_file: rt_file_tif_path=rt_file_path+rt_file mcd15_file_tif_path=mcd15_file_path+mcd15_file drt_out_file_path=out_file_path+"drt/" if not os.path.exists(drt_out_file_path): os.makedirs(drt_out_file_path) drt_out_file_tif_path=drt_out_file_path+rt_hv+".tif" eco_out_file_path=out_file_path+"eco/" if not os.path.exists(eco_out_file_path): os.makedirs(eco_out_file_path) eco_out_file_tif_path=eco_out_file_path+rt_hv+".tif" wat_out_file_path=out_file_path+"wat/" if not os.path.exists(wat_out_file_path): os.makedirs(wat_out_file_path) wat_out_file_tif_path=wat_out_file_path+rt_hv+".tif" tim_out_file_path=out_file_path+"tim/" if not os.path.exists(tim_out_file_path): os.makedirs(tim_out_file_path) tim_out_file_tif_path=tim_out_file_path+rt_hv+".tif" rt_raster=gdal.Open(rt_file_tif_path) rt_raster_array=rt_raster.ReadAsArray() rt_lai_array=rt_raster_array[0] rt_qa_array=rt_raster_array[1] rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array) rt_lai_array_fin=rt_lai_array_mask*0.001 mcd15_raster=gdal.Open(mcd15_file_tif_path) mcd15_sub_dataset=mcd15_raster.GetSubDatasets() # for sub_dataset in mcd15_sub_dataset: # print(sub_dataset[1]) # print(mcd15_sub_dataset[1][1]) # print(mcd15_sub_dataset[2][1]) mcd15_sub_lai=gdal.Open(mcd15_sub_dataset[1][0]) mcd15_sub_qc=gdal.Open(mcd15_sub_dataset[2][0]) mcd15_lai_array=mcd15_sub_lai.ReadAsArray() mcd15_qc_array=mcd15_sub_qc.ReadAsArray() mcd15_lai_array_mask=np.where(mcd15_lai_array>248,np.nan,mcd15_lai_array) mcd15_lai_array_fin=mcd15_lai_array_mask*0.1 rt_qa_array_bin=copy.copy(rt_qa_array) rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape for i in range(rt_qa_array_row): for j in range(rt_qa_array_col): rt_qa_array_bin[j]="{:012b}".format(rt_qa_array_bin[j])[-4:] mcd15_qc_array_bin=copy.copy(mcd15_qc_array) mcd15_qc_array_row,mcd15_qc_array_col=mcd15_qc_array.shape for i in range(mcd15_qc_array_row): for j in range(mcd15_qc_array_col): mcd15_qc_array_bin[j]="{:08b}".format(mcd15_qc_array[j])[-1:] mcd15_lai_main_array=np.where(mcd15_qc_array_bin==1,np.nan,mcd15_lai_array_fin) lai_dif=rt_lai_array_fin-mcd15_lai_main_array lai_dif=lai_dif*1000 drt_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11), np.nan,lai_dif) eco_lai_dif_array=np.where((rt_qa_array_bin