Python GDAL讀取柵格資料並基於質量評估波段QA對指定資料加以篩選掩膜

瘋狂學習GIS發表於2023-02-23

本文介紹基於Python語言中gdal模組，對遙感影像資料進行柵格讀取與計算，同時基於QA波段對像元加以篩選、掩膜的操作。

本文所要實現的需求具體為：現有自行計算的全球葉面積指數（LAI）.tif格式柵格產品（下稱“自有產品”），為了驗證其精確度，需要與已有學者提出的成熟產品——GLASS全球LAI.hdf格式柵格產品（下稱“GLASS產品”）進行做差對比；其中，自有產品除了LAI波段外，還有一個質量評估波段（QA），即自有產品在後期使用時，還需結合QA波段進行篩選、掩膜等處理。其中，二者均為基於MODIS hv分幅的產品。

本文分為兩部分，第一部分為程式碼的詳細分段講解，第二部分為完整程式碼。

1 程式碼分段講解

1.1 模組與路徑準備

首先，需要對用到的模組與存放柵格影像的各類路徑加以準備。

import os
import copy
import numpy as np
import pylab as plt
from osgeo import gdal

# rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Rc_Lai_A2018161_h12v03.tif"
# gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h12v03.2020323.hdf"
# out_file_path="G:/Postgraduate/LAI_Glass_RTlab/test.tif"
rt_file_path="I:/LAI_RTLab/A2018161/"
gl_file_path="I:/LAI_Glass/2018161/"
out_file_path="I:/LAI_Dif/"

其中，rt_file_path為自有產品的存放路徑，gl_file_path為GLASS產品的存放路徑，out_file_path為最終二者柵格做完差值處理後結果的存放路徑。

1.2 柵格影像檔名讀取與配對

接下來，需要將全部待處理的柵格影像用os.listdir()進行獲取，並用for迴圈進行迴圈批次處理操作的準備。

rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
    file_name_split=rt_file.split("_")
    rt_hv=file_name_split[3][:-4]
    
    gl_file_list=os.listdir(gl_file_path)
    for gl_file in gl_file_list:
        if rt_hv in gl_file:
            rt_file_tif_path=rt_file_path+rt_file
            gl_file_tif_path=gl_file_path+gl_file

其中，由於本文需求是對兩種產品做差，因此首先需要結合二者的hv分幅編號，將同一分幅編號的兩景遙感影像放在一起；因此，依據自有產品檔名的特徵，選擇.split()進行字串分割，並隨後擷取獲得遙感影像的hv分幅編號。

1.3 輸出檔名稱準備

前述1.1部分已經配置好了輸出檔案存放的路徑，但是還沒有進行輸出檔案檔名的配置；因此這裡我們需要配置好每一個做差後的遙感影像的檔案存放路徑與名稱。其中，我們就直接以遙感影像的hv編號作為輸出結果檔名。

            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=os.path.join(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=os.path.join(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=os.path.join(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=os.path.join(tim_out_file_path,rt_hv+".tif")

這一部分程式碼分為了四個部分，是因為自有產品的LAI是分別依據四種演算法得到的，在做差時需要每一種演算法分別和GLASS產品進行相減，因此配置了四個輸出路徑資料夾。

1.4 柵格檔案資料與資訊讀取

接下來，利用gdal模組對.tif與.hdf等兩種柵格影像加以讀取。

            rt_raster=gdal.Open(rt_file_path+rt_file)
            rt_band_num=rt_raster.RasterCount
            rt_raster_array=rt_raster.ReadAsArray()
            rt_lai_array=rt_raster_array[0]
            rt_qa_array=rt_raster_array[1]
            rt_lai_band=rt_raster.GetRasterBand(1)
            # rt_lai_nodata=rt_lai_band.GetNoDataValue()
            # rt_lai_nodata=32767
            # rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
            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
            
            gl_raster=gdal.Open(gl_file_path+gl_file)
            gl_band_num=gl_raster.RasterCount
            gl_raster_array=gl_raster.ReadAsArray()
            gl_lai_array=gl_raster_array
            gl_lai_band=gl_raster.GetRasterBand(1)
            gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
            gl_lai_array_fin=gl_lai_array_mask*0.01
            
            row=rt_raster.RasterYSize
            col=rt_raster.RasterXSize
            geotransform=rt_raster.GetGeoTransform()
            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產品加以做差操作，隨後需要對四種演算法分別加以提取。

            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[i][j]="{:012b}".format(rt_qa_array_bin[i][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<100) & (rt_qa_array_bin==111))
            # wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
            # tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
            
            # colormap=plt.cm.Greens
            # plt.figure(1)
            # # plt.subplot(2,4,1)
            # plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
            # plt.title("RT_LAI")
            # plt.colorbar()
            # plt.figure(2)
            # # plt.subplot(2,4,2)
            # plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
            # plt.title("GLASS_LAI")
            # plt.colorbar()
            # plt.figure(3)
            # dif_colormap=plt.cm.get_cmap("Spectral")
            # plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
            # plt.title("Difference_LAI (RT-GLASS)")
            # plt.colorbar()
            
            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<100) | (rt_qa_array_bin==111),
                                       np.nan,lai_dif)
            wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
                                       np.nan,lai_dif)
            tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
                                       np.nan,lai_dif)
            
            # plt.figure(4)
            # plt.imshow(DRT_lai_dif_array)
            # plt.colorbar()
            # plt.figure(5)
            # plt.imshow(eco_lai_dif_array)
            # plt.colorbar()
            # plt.figure(6)
            # plt.imshow(wat_lai_dif_array)
            # plt.colorbar()
            # plt.figure(7)
            # plt.imshow(tim_lai_dif_array)
            # plt.colorbar()

其中，上述程式碼前兩句為差值計算與資料化整。將資料轉換為整數，可以減少結果資料圖層的資料量（因為不需要儲存小數了）。

隨後，開始依據QA波段進行資料篩選與掩膜。其實各類遙感影像（例如MODIS、Landsat等）的QA波段都是比較近似的：透過一串二進位制碼來表示遙感影像的質量、資訊等，其中不同的位元位往往都代表著一種特性。例如下圖所示為Landsat Collection 2 Level-2的QA波段含義。

在這裡，QA波段原本為十進位制（一般遙感影像為了節省空間，QA波段都是寫成十進位制的形式），因此需要將其轉換為二進位制；隨後透過獲取指定需要的二進位制資料位數（在本文中也就是能確定自有產品中這一畫素來自於哪一種演算法的二進位制位數），從而判斷這一畫素所得LAI是透過哪一種演算法得到的，從而將每種演算法對應的畫素分別放在一起處理。DRT_lai_dif_array等四個變數分別表示四種演算法中，除了自己這一種演算法得到的畫素之外的其他所有畫素；之所以選擇這種方式，是因為後期我們可以將其直接掩膜掉，那麼剩下的就是這種演算法自身的畫素了。

其中，上述程式碼註釋掉的plt相關內容可以實現繪製空間分佈圖，大家感興趣可以嘗試使用。

1.6 結果柵格檔案寫入與儲存

接下來，將我們完成上述差值計算與依據演算法進行篩選後的影像儲存。

            driver=gdal.GetDriverByName("Gtiff")
            out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_DRT_lai.SetGeoTransform(geotransform)
            out_DRT_lai.SetProjection(projection)
            out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
            out_DRT_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_eco_lai.SetGeoTransform(geotransform)
            out_eco_lai.SetProjection(projection)
            out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
            out_eco_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_wat_lai.SetGeoTransform(geotransform)
            out_wat_lai.SetProjection(projection)
            out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
            out_wat_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_tim_lai.SetGeoTransform(geotransform)
            out_tim_lai.SetProjection(projection)
            out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
            out_tim_lai=None
            
            print(rt_hv)

其中，.GetDriverByName("Gtiff")表示儲存為.tif格式的GeoTIFF檔案；driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)表示按照路徑、行列數、波段數與資料格式等建立一個新的柵格圖層，作為輸出圖層的框架；其後表示分別將地理投影轉換資訊與畫素具體數值分別賦予這一新建的柵格圖層；最後=None表示將其從記憶體空間中釋放，完成寫入與儲存工作。

2 完整程式碼

本文所需完整程式碼如下：

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 15 19:36:15 2021

@author: fkxxgis
"""

import os
import copy
import numpy as np
import pylab as plt
from osgeo import gdal

# rt_file_path="G:/Postgraduate/LAI_Glass_RTlab/Rc_Lai_A2018161_h12v03.tif"
# gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h12v03.2020323.hdf"
# out_file_path="G:/Postgraduate/LAI_Glass_RTlab/test.tif"
rt_file_path="I:/LAI_RTLab/A2018161/"
gl_file_path="I:/LAI_Glass/2018161/"
out_file_path="I:/LAI_Dif/"

rt_file_list=os.listdir(rt_file_path)
for rt_file in rt_file_list:
    file_name_split=rt_file.split("_")
    rt_hv=file_name_split[3][:-4]
    
    gl_file_list=os.listdir(gl_file_path)
    for gl_file in gl_file_list:
        if rt_hv in gl_file:
            rt_file_tif_path=rt_file_path+rt_file
            gl_file_tif_path=gl_file_path+gl_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=os.path.join(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=os.path.join(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=os.path.join(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=os.path.join(tim_out_file_path,rt_hv+".tif")

            rt_raster=gdal.Open(rt_file_path+rt_file)
            rt_band_num=rt_raster.RasterCount
            rt_raster_array=rt_raster.ReadAsArray()
            rt_lai_array=rt_raster_array[0]
            rt_qa_array=rt_raster_array[1]
            rt_lai_band=rt_raster.GetRasterBand(1)
            # rt_lai_nodata=rt_lai_band.GetNoDataValue()
            # rt_lai_nodata=32767
            # rt_lai_mask=np.ma.masked_equal(rt_lai_array,rt_lai_nodata)
            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
            
            gl_raster=gdal.Open(gl_file_path+gl_file)
            gl_band_num=gl_raster.RasterCount
            gl_raster_array=gl_raster.ReadAsArray()
            gl_lai_array=gl_raster_array
            gl_lai_band=gl_raster.GetRasterBand(1)
            gl_lai_array_mask=np.where(gl_lai_array>1000,np.nan,gl_lai_array)
            gl_lai_array_fin=gl_lai_array_mask*0.01
            
            row=rt_raster.RasterYSize
            col=rt_raster.RasterXSize
            geotransform=rt_raster.GetGeoTransform()
            projection=rt_raster.GetProjection()
            
            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[i][j]="{:012b}".format(rt_qa_array_bin[i][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<100) & (rt_qa_array_bin==111))
            # wat_pixel_pos=np.where((rt_qa_array_bin<1000) & (rt_qa_array_bin==1011))
            # tim_pixel_pos=np.where((rt_qa_array_bin<1100) & (rt_qa_array_bin==1111))
            
            # colormap=plt.cm.Greens
            # plt.figure(1)
            # # plt.subplot(2,4,1)
            # plt.imshow(rt_lai_array_fin,cmap=colormap,interpolation='none')
            # plt.title("RT_LAI")
            # plt.colorbar()
            # plt.figure(2)
            # # plt.subplot(2,4,2)
            # plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation='none')
            # plt.title("GLASS_LAI")
            # plt.colorbar()
            # plt.figure(3)
            # dif_colormap=plt.cm.get_cmap("Spectral")
            # plt.imshow(lai_dif,cmap=dif_colormap,interpolation='none')
            # plt.title("Difference_LAI (RT-GLASS)")
            # plt.colorbar()
            
            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<100) | (rt_qa_array_bin==111),
                                       np.nan,lai_dif)
            wat_lai_dif_array=np.where((rt_qa_array_bin<1000) | (rt_qa_array_bin==1011),
                                       np.nan,lai_dif)
            tim_lai_dif_array=np.where((rt_qa_array_bin<1100) | (rt_qa_array_bin==1111),
                                       np.nan,lai_dif)
            
            # plt.figure(4)
            # plt.imshow(DRT_lai_dif_array)
            # plt.colorbar()
            # plt.figure(5)
            # plt.imshow(eco_lai_dif_array)
            # plt.colorbar()
            # plt.figure(6)
            # plt.imshow(wat_lai_dif_array)
            # plt.colorbar()
            # plt.figure(7)
            # plt.imshow(tim_lai_dif_array)
            # plt.colorbar()
            
            driver=gdal.GetDriverByName("Gtiff")
            out_DRT_lai=driver.Create(DRT_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_DRT_lai.SetGeoTransform(geotransform)
            out_DRT_lai.SetProjection(projection)
            out_DRT_lai.GetRasterBand(1).WriteArray(DRT_lai_dif_array)
            out_DRT_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_eco_lai=driver.Create(eco_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_eco_lai.SetGeoTransform(geotransform)
            out_eco_lai.SetProjection(projection)
            out_eco_lai.GetRasterBand(1).WriteArray(eco_lai_dif_array)
            out_eco_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_wat_lai=driver.Create(wat_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_wat_lai.SetGeoTransform(geotransform)
            out_wat_lai.SetProjection(projection)
            out_wat_lai.GetRasterBand(1).WriteArray(wat_lai_dif_array)
            out_wat_lai=None
            
            driver=gdal.GetDriverByName("Gtiff")
            out_tim_lai=driver.Create(tim_out_file_tif_path,row,col,1,gdal.GDT_Float32)
            out_tim_lai.SetGeoTransform(geotransform)
            out_tim_lai.SetProjection(projection)
            out_tim_lai.GetRasterBand(1).WriteArray(tim_lai_dif_array)
            out_tim_lai=None
            
            print(rt_hv)

至此，大功告成。

Python基於Excel資料加以反距離加權空間插值並掩膜圖層
2024-04-10
PythonExcel
Python批次讀取HDF多波段柵格資料並繪製像元直方圖
2023-03-01
Python直方圖
試驗資料的篩選和質量視覺化
2021-02-22
視覺化
Python 柵格資料處理教程（一）
2024-08-13
Python
Python 柵格資料處理教程（二）
2024-08-14
Python
億信華辰:怎樣去斷定一份資料的質量高低？資料質量如何評估？
2022-01-18
讀資料質量管理：資料可靠性與資料質量問題解決之道01資料質量
2024-11-12
C# 讀取四波段遙感影像生成植被覆蓋度柵格(TIF)
2024-11-15
C#
基於python的大資料分析-pandas資料讀取（程式碼實戰）
2019-08-29
Python大資料
讀資料質量管理：資料可靠性與資料質量問題解決之道02資料湖倉
2024-11-13
讀資料質量管理：資料可靠性與資料質量問題解決之道03資料目錄
2024-11-14
讀資料質量管理：資料可靠性與資料質量問題解決之道06資料測試
2024-11-17
讀資料質量管理：資料可靠性與資料質量問題解決之道10資料平臺
2024-11-21
在DataSet資料集中 DataView篩選資料
2024-06-14
View
GeoServer釋出PostGIS資料庫中的柵格資料
2018-08-02
Server資料庫
資料治理--資料質量
2024-06-04
Excel上傳並讀取資料
2021-09-09
Excel
基於MaxCompute的數倉資料質量管理
2019-04-18
讀資料質量管理：資料可靠性與資料質量問題解決之道05資料標準化
2024-11-16
讀資料質量管理：資料可靠性與資料質量問題解決之道09資料可靠性
2024-11-20
Python依據遙感影像的分幅篩選出對應的柵格檔案
2024-06-01
Python
可用於資料庫對比評估的FURPS+模型
2023-03-03
資料庫模型
python爬取股票資料並存到資料庫
2021-03-29
Python資料庫
python讀取MySQL資料
2021-01-02
PythonMySql
KVC原理與資料篩選
2022-11-29
R：SNP資料篩選.R
2024-06-19
讀資料質量管理：資料可靠性與資料質量問題解決之道04收集與清洗
2024-11-15
用Python計算柵格資料的真實面積
2024-11-11
Python
MongoDB資料庫順序讀效能評估測試
2016-08-11
MongoDB資料庫
讀取mysq資料庫l資料，並使用dataview顯示
2013-12-04
資料庫View
讀資料質量管理：資料可靠性與資料質量問題解決之道07異常檢測
2024-11-18
用好工業大資料的基礎是資料質量
2023-10-09
大資料
資料資產價值評估常用方法及對比
2020-12-21
資料治理：資料質量管理策略！
2023-12-14
資料治理之資料質量管理
2024-04-25
Python爬取CSDN部落格資料
2019-01-03
Python
Python讀取YAML配置資料
2024-05-22
PythonYAML
python讀取串列埠資料
2024-08-27
Python串列埠