Python利用GDAL模块实现读取栅格数据并对指定数据加以筛选掩膜

孟子桓

本文介绍基于Python语言中

1. gdal

复制代码

模块，对遥感影像数据进行栅格读取与计算，同时基于QA波段对像元加以筛选、掩膜的操作。
本文所要实现的需求具体为：现有自行计算的全球叶面积指数（LAI）

1. .tif

复制代码

格式栅格产品（下称“自有产品”），为了验证其精确度，需要与已有学者提出的成熟产品——GLASS全球LAI

1. .hdf

复制代码

格式栅格产品（下称“GLASS产品”）进行做差对比；其中，自有产品除了LAI波段外，还有一个质量评估波段（QA），即自有产品在后期使用时，还需结合QA波段进行筛选、掩膜等处理。其中，二者均为基于MODIS

1. 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
   
6. 
   
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/"

复制代码

其中，

1. rt_file_path

复制代码

为自有产品的存放路径，

1. gl_file_path

复制代码

为GLASS产品的存放路径，

1. out_file_path

复制代码

为最终二者栅格做完差值处理后结果的存放路径。

1.2 栅格图像文件名读取与配对

接下来，需要将全部待处理的栅格图像用

1. os.listdir()

复制代码

进行获取，并用

1. 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

复制代码

其中，由于本文需求是对两种产品做差，因此首先需要结合二者的

1. hv

复制代码

分幅编号，将同一分幅编号的两景遥感影像放在一起；因此，依据自有产品文件名的特征，选择

1. .split()

复制代码

进行字符串分割，并随后截取获得遥感影像的

1. hv

复制代码

分幅编号。

1.3 输出文件名称准备

前述1.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 栅格文件数据与信息读取

接下来，利用

1. gdal

复制代码

模块对

1. .tif

复制代码

与

1. .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()

复制代码

首先，以上述代码的第一段为例进行讲解。其中，

1. gdal.Open()

复制代码

读取栅格图像；

1. .RasterCount

复制代码

获取栅格图像波段数量；

1. .ReadAsArray()

复制代码

将栅格图像各波段的信息读取为

1. Array

复制代码

格式，当波段数量大于

1. 1

复制代码

时，其共有三维，第一维为波段的个数；

1. rt_raster_array[0]

复制代码

表示取

1. Array

复制代码

中的第一个波段，在本文中也就是自有产品的LAI波段；

1. rt_qa_array=rt_raster_array[1]

复制代码

则表示取出第二个波段，在本文中也就是自有产品的QA波段；

1. .GetRasterBand(1)

复制代码

表示获取栅格图像中的第一个波段（注意，这里序号不是从

1. 0

复制代码

开始而是从

1. 1

复制代码

开始）；

1. np.where(rt_lai_array>30000,np.nan,rt_lai_array)

复制代码

表示利用

1. np.where()

复制代码

函数对

1. Array

复制代码

中第一个波段中像素

1. >30000

复制代码

加以选取，并将其设置为

1. nan

复制代码

，其他值不变。这一步骤是消除图像中填充值、

1. Nodata

复制代码

值的方法。最后一句

1. *0.001

复制代码

是将图层原有的缩放系数复原。
其次，上述代码第三段为获取栅格行、列数与投影变换信息。

1.5 差值计算与QA波段筛选

接下来，首先对自有产品与GLASS产品加以做差操作，随后需要对四种算法分别加以提取。

1.             lai_dif=rt_lai_array_fin-gl_lai_array_fin
   
2.             lai_dif=lai_dif*1000
   
3.             
   
4.             rt_qa_array_bin=copy.copy(rt_qa_array)
   
5.             rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
   
6.             for i in range(rt_qa_array_row):
   
7.                 for j in range(rt_qa_array_col):
   
8.                     rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
   
9.                     
   
10.             # DRT_pixel_pos=np.where((rt_qa_array_bin>=100) & (rt_qa_array_bin==11))
    
11.             # eco_pixel_pos=np.where((rt_qa_array_bin<100) & (rt_qa_array_bin==111))
    
12.             # wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
    
13.             # tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
    
14.             
    
15.             # colormap=plt.cm.Greens
    
16.             # plt.figure(1)
    
17.             # # plt.subplot(2,4,1)
    
18.             # plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
    
19.             # plt.title("RT_LAI")
    
20.             # plt.colorbar()
    
21.             # plt.figure(2)
    
22.             # # plt.subplot(2,4,2)
    
23.             # plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
    
24.             # plt.title("GLASS_LAI")
    
25.             # plt.colorbar()
    
26.             # plt.figure(3)
    
27.             # dif_colormap=plt.cm.get_cmap("Spectral")
    
28.             # plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
    
29.             # plt.title("Difference_LAI (RT-GLASS)")
    
30.             # plt.colorbar()
    
31.             
    
32.             DRT_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
    
33.                                        np.nan,lai_dif)
    
34.             eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
    
35.                                        np.nan,lai_dif)
    
36.             wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
    
37.                                        np.nan,lai_dif)
    
38.             tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
    
39.                                        np.nan,lai_dif)
    
40.             
    
41.             # plt.figure(4)
    
42.             # plt.imshow(DRT_lai_dif_array)
    
43.             # plt.colorbar()
    
44.             # plt.figure(5)
    
45.             # plt.imshow(eco_lai_dif_array)
    
46.             # plt.colorbar()
    
47.             # plt.figure(6)
    
48.             # plt.imshow(wat_lai_dif_array)
    
49.             # plt.colorbar()
    
50.             # plt.figure(7)
    
51.             # plt.imshow(tim_lai_dif_array)
    
52.             # plt.colorbar()

复制代码

其中，上述代码前两句为差值计算与数据化整。将数据转换为整数，可以减少结果数据图层的数据量（因为不需要存储小数了）。
随后，开始依据QA波段进行数据筛选与掩膜。其实各类遥感影像（例如MODIS、Landsat等）的QA波段都是比较近似的：通过一串二进制码来表示遥感影像的质量、信息等，其中不同的比特位往往都代表着一种特性。例如下图所示为Landsat Collection 2 Level-2的QA波段含义。

在这里，QA波段原本为十进制（一般遥感影像为了节省空间，QA波段都是写成十进制的形式），因此需要将其转换为二进制；随后通过获取指定需要的二进制数据位数（在本文中也就是能确定自有产品中这一像素来自于哪一种算法的二进制位数），从而判断这一像素所得LAI是通过哪一种算法得到的，从而将每种算法对应的像素分别放在一起处理。

1. DRT_lai_dif_array

复制代码

等四个变量分别表示四种算法中，除了自己这一种算法得到的像素之外的其他所有像素；之所以选择这种方式，是因为后期我们可以将其直接掩膜掉，那么剩下的就是这种算法自身的像素了。
其中，上述代码注释掉的

1. plt

复制代码

相关内容可以实现绘制空间分布图，大家感兴趣可以尝试使用。

1.6 结果栅格文件写入与保存

接下来，将我们完成上述差值计算与依据算法进行筛选后的图像保存。

1.             driver=gdal.GetDriverByName("Gtiff")
   
2.             out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
   
3.             out_DRT_lai.SetGeoTransform(geotransform)
   
4.             out_DRT_lai.SetProjection(projection)
   
5.             out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
   
6.             out_DRT_lai=None
   
7.             
   
8.             driver=gdal.GetDriverByName("Gtiff")
   
9.             out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
   
10.             out_eco_lai.SetGeoTransform(geotransform)
    
11.             out_eco_lai.SetProjection(projection)
    
12.             out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
    
13.             out_eco_lai=None
    
14.             
    
15.             driver=gdal.GetDriverByName("Gtiff")
    
16.             out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
    
17.             out_wat_lai.SetGeoTransform(geotransform)
    
18.             out_wat_lai.SetProjection(projection)
    
19.             out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
    
20.             out_wat_lai=None
    
21.             
    
22.             driver=gdal.GetDriverByName("Gtiff")
    
23.             out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
    
24.             out_tim_lai.SetGeoTransform(geotransform)
    
25.             out_tim_lai.SetProjection(projection)
    
26.             out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
    
27.             out_tim_lai=None
    
28.             
    
29.             print(rt_hv)

复制代码

其中，

1. .GetDriverByName("Gtiff")

复制代码

表示保存为

1. .tif

复制代码

格式的GeoTIFF文件；

1. driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)

复制代码

表示按照路径、行列数、波段数与数据格式等建立一个新的栅格图层，作为输出图层的框架；其后表示分别将地理投影转换信息与像素具体数值分别赋予这一新建的栅格图层；最后

1. =None

复制代码

表示将其从内存空间中释放，完成写入与保存工作。

2 完整代码

本文所需完整代码如下：

1. # -*- coding: utf-8 -*-"""Created on Thu Jul 15 19:36:15 2021@author: fkxxgis"""import os
   
2. import copy
   
3. import numpy as np
   
4. import pylab as plt
   
5. from osgeo import gdal
   
6. 
   
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_list=os.listdir(rt_file_path)
    
13. for rt_file in rt_file_list:
    
14.     file_name_split=rt_file.split("_")
    
15.     rt_hv=file_name_split[3][:-4]
    
16.    
    
17.     gl_file_list=os.listdir(gl_file_path)
    
18.     for gl_file in gl_file_list:
    
19.         if rt_hv in gl_file:
    
20.             rt_file_tif_path=rt_file_path+rt_file
    
21.             gl_file_tif_path=gl_file_path+gl_file                        DRT_out_file_path=out_file_path+"DRT/"
    
22.             if not os.path.exists(DRT_out_file_path):
    
23.                 os.makedirs(DRT_out_file_path)
    
24.             DRT_out_file_tif_path=os.path.join(DRT_out_file_path,rt_hv+".tif")
    
25.             
    
26.             eco_out_file_path=out_file_path+"eco/"
    
27.             if not os.path.exists(eco_out_file_path):
    
28.                 os.makedirs(eco_out_file_path)
    
29.             eco_out_file_tif_path=os.path.join(eco_out_file_path,rt_hv+".tif")
    
30.             
    
31.             wat_out_file_path=out_file_path+"wat/"
    
32.             if not os.path.exists(wat_out_file_path):
    
33.                 os.makedirs(wat_out_file_path)
    
34.             wat_out_file_tif_path=os.path.join(wat_out_file_path,rt_hv+".tif")
    
35.             
    
36.             tim_out_file_path=out_file_path+"tim/"
    
37.             if not os.path.exists(tim_out_file_path):
    
38.                 os.makedirs(tim_out_file_path)
    
39.             tim_out_file_tif_path=os.path.join(tim_out_file_path,rt_hv+".tif")            rt_raster=gdal.Open(rt_file_path+rt_file)
    
40.             rt_band_num=rt_raster.RasterCount
    
41.             rt_raster_array=rt_raster.ReadAsArray()
    
42.             rt_lai_array=rt_raster_array[0]
    
43.             rt_qa_array=rt_raster_array[1]
    
44.             rt_lai_band=rt_raster.GetRasterBand(1)
    
45.             # rt_lai_nodata=rt_lai_band.GetNoDataValue()
    
46.             # rt_lai_nodata=32767
    
47.             # rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
    
48.             rt_lai_array_mask=np.where(rt_lai_array>30000,np.nan,rt_lai_array)
    
49.             rt_lai_array_fin=rt_lai_array_mask*0.001
    
50.             
    
51.             gl_raster=gdal.Open(gl_file_path+gl_file)
    
52.             gl_band_num=gl_raster.RasterCount
    
53.             gl_raster_array=gl_raster.ReadAsArray()
    
54.             gl_lai_array=gl_raster_array
    
55.             gl_lai_band=gl_raster.GetRasterBand(1)
    
56.             gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
    
57.             gl_lai_array_fin=gl_lai_array_mask*0.01
    
58.             
    
59.             row=rt_raster.RasterYSize
    
60.             col=rt_raster.RasterXSize
    
61.             geotransform=rt_raster.GetGeoTransform()
    
62.             projection=rt_raster.GetProjection()                        lai_dif=rt_lai_array_fin-gl_lai_array_fin
    
63.             lai_dif=lai_dif*1000
    
64.             
    
65.             rt_qa_array_bin=copy.copy(rt_qa_array)
    
66.             rt_qa_array_row,rt_qa_array_col=rt_qa_array.shape
    
67.             for i in range(rt_qa_array_row):
    
68.                 for j in range(rt_qa_array_col):
    
69.                     rt_qa_array_bin[i][j]="{:012b}".format(rt_qa_array_bin[i][j])[-4:]
    
70.                     
    
71.             # DRT_pixel_pos=np.where((rt_qa_array_bin>=100) & (rt_qa_array_bin==11))
    
72.             # eco_pixel_pos=np.where((rt_qa_array_bin<100) & (rt_qa_array_bin==111))
    
73.             # wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
    
74.             # tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
    
75.             
    
76.             # colormap=plt.cm.Greens
    
77.             # plt.figure(1)
    
78.             # # plt.subplot(2,4,1)
    
79.             # plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
    
80.             # plt.title("RT_LAI")
    
81.             # plt.colorbar()
    
82.             # plt.figure(2)
    
83.             # # plt.subplot(2,4,2)
    
84.             # plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
    
85.             # plt.title("GLASS_LAI")
    
86.             # plt.colorbar()
    
87.             # plt.figure(3)
    
88.             # dif_colormap=plt.cm.get_cmap("Spectral")
    
89.             # plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
    
90.             # plt.title("Difference_LAI (RT-GLASS)")
    
91.             # plt.colorbar()
    
92.             
    
93.             DRT_lai_dif_array=np.where((rt_qa_array_bin>=100) | (rt_qa_array_bin==11),
    
94.                                        np.nan,lai_dif)
    
95.             eco_lai_dif_array=np.where((rt_qa_array_bin<100) | (rt_qa_array_bin==111),
    
96.                                        np.nan,lai_dif)
    
97.             wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
    
98.                                        np.nan,lai_dif)
    
99.             tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
    
100.                                        np.nan,lai_dif)
     
101.             
     
102.             # plt.figure(4)
     
103.             # plt.imshow(DRT_lai_dif_array)
     
104.             # plt.colorbar()
     
105.             # plt.figure(5)
     
106.             # plt.imshow(eco_lai_dif_array)
     
107.             # plt.colorbar()
     
108.             # plt.figure(6)
     
109.             # plt.imshow(wat_lai_dif_array)
     
110.             # plt.colorbar()
     
111.             # plt.figure(7)
     
112.             # plt.imshow(tim_lai_dif_array)
     
113.             # plt.colorbar()                        driver=gdal.GetDriverByName("Gtiff")
     
114.             out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
     
115.             out_DRT_lai.SetGeoTransform(geotransform)
     
116.             out_DRT_lai.SetProjection(projection)
     
117.             out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
     
118.             out_DRT_lai=None
     
119.             
     
120.             driver=gdal.GetDriverByName("Gtiff")
     
121.             out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
     
122.             out_eco_lai.SetGeoTransform(geotransform)
     
123.             out_eco_lai.SetProjection(projection)
     
124.             out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
     
125.             out_eco_lai=None
     
126.             
     
127.             driver=gdal.GetDriverByName("Gtiff")
     
128.             out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
     
129.             out_wat_lai.SetGeoTransform(geotransform)
     
130.             out_wat_lai.SetProjection(projection)
     
131.             out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
     
132.             out_wat_lai=None
     
133.             
     
134.             driver=gdal.GetDriverByName("Gtiff")
     
135.             out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
     
136.             out_tim_lai.SetGeoTransform(geotransform)
     
137.             out_tim_lai.SetProjection(projection)
     
138.             out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
     
139.             out_tim_lai=None
     
140.             
     
141.             print(rt_hv)

复制代码

到此这篇关于Python利用GDAL模块实现读取栅格数据并对指定数据加以筛选掩膜的文章就介绍到这了,更多相关Python GDAL读取栅格数据内容请搜索脚本之家以前的文章或继续浏览下面的相关文章希望大家以后多多支持脚本之家！

来源:https://www.jb51.net/article/276257.htm
免责声明：由于采集信息均来自互联网，如果侵犯了您的权益，请联系我们【E-Mail:cb@itdo.tech】 我们会及时删除侵权内容，谢谢合作！