社群
學習
工具庫
AI工具
休閒
搜尋
繁体中文
目錄
1 程式碼分段講解
1.1 模組與路徑準備
1.2 柵格圖像檔案名稱讀取與配對
1.3 輸出檔案名稱準備
1.4 柵格檔案資料與資訊讀取
1.5 差值計算與QA波段篩選
1.6 结果栅格文件写入与保存
2 完整代码
首頁 後端開發 Python教學 Python如何利用GDAL模組實作讀取柵格資料並對指定資料加以篩選遮罩

Python如何利用GDAL模組實作讀取柵格資料並對指定資料加以篩選遮罩

王林
May 15, 2023 am 11:16 AM
python gdal

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_h22v03.tif"
# gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h22v03.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=&#39;none&#39;)
            # plt.title("RT_LAI")
            # plt.colorbar()
            # plt.figure(2)
            # # plt.subplot(2,4,2)
            # plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation=&#39;none&#39;)
            # plt.title("GLASS_LAI")
            # plt.colorbar()
            # plt.figure(3)
            # dif_colormap=plt.cm.get_cmap("Spectral")
            # plt.imshow(lai_dif,cmap=dif_colormap,interpolation=&#39;none&#39;)
            # 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_h22v03.tif"
# gl_file_path="G:/Postgraduate/LAI_Glass_RTlab/GLASS01E01.V50.A2018161.h22v03.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=&#39;none&#39;)
            # plt.title("RT_LAI")
            # plt.colorbar()
            # plt.figure(2)
            # # plt.subplot(2,4,2)
            # plt.imshow(gl_lai_array_fin,cmap=colormap,interpolation=&#39;none&#39;)
            # plt.title("GLASS_LAI")
            # plt.colorbar()
            # plt.figure(3)
            # dif_colormap=plt.cm.get_cmap("Spectral")
            # plt.imshow(lai_dif,cmap=dif_colormap,interpolation=&#39;none&#39;)
            # 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如何利用GDAL模組實作讀取柵格資料並對指定資料加以篩選遮罩的詳細內容。更多資訊請關注PHP中文網其他相關文章！

本網站聲明
本文內容由網友自願投稿，版權歸原作者所有。本站不承擔相應的法律責任。如發現涉嫌抄襲或侵權的內容，請聯絡admin@php.cn

熱AI工具

Undresser.AI Undress

Undresser.AI Undress

人工智慧驅動的應用程序，用於創建逼真的裸體照片

AI Clothes Remover

AI Clothes Remover

用於從照片中去除衣服的線上人工智慧工具。

Undress AI Tool

Undress AI Tool

免費脫衣圖片

Clothoff.io

Clothoff.io

AI脫衣器

AI Hentai Generator

AI Hentai Generator

免費產生 AI 無盡。

顯示更多

熱門文章

R.E.P.O.能量晶體解釋及其做什麼（黃色晶體）
2 週前 By 尊渡假赌尊渡假赌尊渡假赌
倉庫：如何復興隊友
4 週前 By 尊渡假赌尊渡假赌尊渡假赌
Hello Kitty Island冒險：如何獲得巨型種子
3 週前 By 尊渡假赌尊渡假赌尊渡假赌
擊敗分裂小說需要多長時間？
3 週前 By DDD
R.E.P.O.保存文件位置：在哪里以及如何保護它？
3 週前 By DDD
顯示更多

熱工具

記事本++7.3.1

記事本++7.3.1

好用且免費的程式碼編輯器

SublimeText3漢化版

SublimeText3漢化版

中文版，非常好用

禪工作室 13.0.1

禪工作室 13.0.1

強大的PHP整合開發環境

Dreamweaver CS6

Dreamweaver CS6

視覺化網頁開發工具

SublimeText3 Mac版

SublimeText3 Mac版

神級程式碼編輯軟體(SublimeText3)

顯示更多

熱門話題

gmail信箱登陸入口在哪裡
7314
9
Java教學
1625
14
CakePHP 教程
1348
46
Laravel 教程
1260
25
PHP教程
1207
29
顯示更多
Related knowledge
模板化的優點和缺點有哪些？ 模板化的優點和缺點有哪些？ May 08, 2024 pm 03:51 PM

模板化的優點和缺點有哪些？

Google AI 為開發者發佈 Gemini 1.5 Pro 和 Gemma 2 Google AI 為開發者發佈 Gemini 1.5 Pro 和 Gemma 2 Jul 01, 2024 am 07:22 AM

Google AI 為開發者發佈 Gemini 1.5 Pro 和 Gemma 2

怎麼下載deepseek 小米 怎麼下載deepseek 小米 Feb 19, 2025 pm 05:27 PM

怎麼下載deepseek 小米

deepseek怎麼問他 deepseek怎麼問他 Feb 19, 2025 pm 04:42 PM

deepseek怎麼問他

evaluate函數怎麼保存 evaluate函數怎麼保存 May 07, 2024 am 01:09 AM

evaluate函數怎麼保存

NET40是什麼軟體 NET40是什麼軟體 May 10, 2024 am 01:12 AM

NET40是什麼軟體

deepseek該怎麼搜索 deepseek該怎麼搜索 Feb 19, 2025 pm 05:18 PM

deepseek該怎麼搜索

瀏覽器插件是什麼語言寫的 瀏覽器插件是什麼語言寫的 May 08, 2024 pm 09:36 PM

瀏覽器插件是什麼語言寫的

See all articles