Python GDAL读取栅格数据并基于质量评估波段QA对指定数据加以筛选掩膜

正义侠士

  本文介绍基于Python语言中gdal模块，对遥感影像数据进行栅格读取与计算，同时基于QA波段对像元加以筛选、掩膜的操作。
  本文所要实现的需求具体为：现有自行计算的全球叶面积指数（LAI）.tif格式栅格产品（下称“自有产品”），为了验证其精确度，需要与已有学者提出的成熟产品——GLASS全球LAI.hdf格式栅格产品（下称“GLASS产品”）进行做差对比；其中，自有产品除了LAI波段外，还有一个质量评估波段（QA），即自有产品在后期使用时，还需结合QA波段进行筛选、掩膜等处理。其中，二者均为基于MODIS hv分幅的产品。
  本文分为两部分，第一部分为代码的详细分段讲解，第二部分为完整代码。
1 代码分段讲解

1.1 模块与路径准备

  首先，需要对用到的模块与存放栅格图像的各类路径加以准备。

1. import os
2. import copy
3. import numpy as np
4. import pylab as plt
5. from osgeo import gdal
7. # rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Rc_Lai_A2018161_h12v03.tif"
8. # gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h12v03.2020323.hdf"
9. # out_file_path="G:/Postgraduate/LAI_Glass_RTlab/test.tif"
10. rt_file_path="I:/LAI_RTLab/A2018161/"
11. gl_file_path="I:/LAI_Glass/2018161/"
12. out_file_path="I:/LAI_Dif/"

复制代码

  其中，rt_file_path为自有产品的存放路径，gl_file_path为GLASS产品的存放路径，out_file_path为最终二者栅格做完差值处理后结果的存放路径。
1.2 栅格图像文件名读取与配对

  接下来，需要将全部待处理的栅格图像用os.listdir()进行获取，并用for循环进行循环批量处理操作的准备。

1. rt_file_list=os.listdir(rt_file_path)
2. for rt_file in rt_file_list:
3.     file_name_split=rt_file.split("_")
4.     rt_hv=file_name_split[3][:-4]
5.    
6.     gl_file_list=os.listdir(gl_file_path)
7.     for gl_file in gl_file_list:
8.         if rt_hv in gl_file:
9.             rt_file_tif_path=rt_file_path+rt_file
10.             gl_file_tif_path=gl_file_path+gl_file

复制代码

  其中，由于本文需求是对两种产品做差，因此首先需要结合二者的hv分幅编号，将同一分幅编号的两景遥感影像放在一起；因此，依据自有产品文件名的特征，选择.split()进行字符串分割，并随后截取获得遥感影像的hv分幅编号。
1.3 输出文件名称准备

  前述1.1部分已经配置好了输出文件存放的路径，但是还没有进行输出文件文件名的配置；因此这里我们需要配置好每一个做差后的遥感影像的文件存放路径与名称。其中，我们就直接以遥感影像的hv编号作为输出结果文件名。

1.             DRT_out_file_path=out_file_path+"DRT/"
2.             if not os.path.exists(DRT_out_file_path):
3.                 os.makedirs(DRT_out_file_path)
4.             DRT_out_file_tif_path=os.path.join(DRT_out_file_path,rt_hv+".tif")
5.             
6.             eco_out_file_path=out_file_path+"eco/"
7.             if not os.path.exists(eco_out_file_path):
8.                 os.makedirs(eco_out_file_path)
9.             eco_out_file_tif_path=os.path.join(eco_out_file_path,rt_hv+".tif")
10.             
11.             wat_out_file_path=out_file_path+"wat/"
12.             if not os.path.exists(wat_out_file_path):
13.                 os.makedirs(wat_out_file_path)
14.             wat_out_file_tif_path=os.path.join(wat_out_file_path,rt_hv+".tif")
15.             
16.             tim_out_file_path=out_file_path+"tim/"
17.             if not os.path.exists(tim_out_file_path):
18.                 os.makedirs(tim_out_file_path)
19.             tim_out_file_tif_path=os.path.join(tim_out_file_path,rt_hv+".tif")

复制代码

  这一部分代码分为了四个部分，是因为自有产品的LAI是分别依据四种算法得到的，在做差时需要每一种算法分别和GLASS产品进行相减，因此配置了四个输出路径文件夹。
1.4 栅格文件数据与信息读取

  接下来，利用gdal模块对.tif与.hdf等两种栅格图像加以读取。

1.             rt_raster=gdal.Open(rt_file_path+rt_file)
2.             rt_band_num=rt_raster.RasterCount
3.             rt_raster_array=rt_raster.ReadAsArray()
4.             rt_lai_array=rt_raster_array[0]
5.             rt_qa_array=rt_raster_array[1]
6.             rt_lai_band=rt_raster.GetRasterBand(1)
7.             # rt_lai_nodata=rt_lai_band.GetNoDataValue()
8.             # rt_lai_nodata=32767
9.             # rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
10.             rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
11.             rt_lai_array_fin=rt_lai_array_mask*0.001
12.             
13.             gl_raster=gdal.Open(gl_file_path+gl_file)
14.             gl_band_num=gl_raster.RasterCount
15.             gl_raster_array=gl_raster.ReadAsArray()
16.             gl_lai_array=gl_raster_array
17.             gl_lai_band=gl_raster.GetRasterBand(1)
18.             gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
19.             gl_lai_array_fin=gl_lai_array_mask*0.01
20.             
21.             row=rt_raster.RasterYSize
22.             col=rt_raster.RasterXSize
23.             geotransform=rt_raster.GetGeoTransform()
24.             projection=rt_raster.GetProjection()

复制代码

  首先，以上述代码的第一段为例进行讲解。其中，gdal.Open()读取栅格图像；.RasterCount获取栅格图像波段数量；.ReadAsArray()将栅格图像各波段的信息读取为Array格式，当波段数量大于1时，其共有三维，第一维为波段的个数；rt_raster_array[0]表示取Array中的第一个波段，在本文中也就是自有产品的LAI波段；rt_qa_array=rt_raster_array[1]则表示取出第二个波段，在本文中也就是自有产品的QA波段；.GetRasterBand(1)表示获取栅格图像中的第一个波段（注意，这里序号不是从0开始而是从1开始）；np.where(rt_lai_array>30000,np.nan,rt_lai_array)表示利用np.where()函数对Array中第一个波段中像素>30000加以选取，并将其设置为nan，其他值不变。这一步骤是消除图像中填充值、Nodata值的方法。最后一句*0.001是将图层原有的缩放系数复原。
  其次，上述代码第三段为获取栅格行、列数与投影变换信息。
1.5 差值计算与QA波段筛选

  接下来，首先对自有产品与GLASS产品加以做差操作，随后需要对四种算法分别加以提取。
[code]            lai_dif=rt_lai_array_fin-gl_lai_array_fin            lai_dif=lai_dif*1000                        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:]                                # DRT_pixel_pos=np.where((rt_qa_array_bin>=100) & (rt_qa_array_bin==11))            # eco_pixel_pos=np.where((rt_qa_array_bin