本文介紹基于Python中的gdal模塊，批量讀取大量多波段遙感影像文件，分別對(duì)各波段數(shù)據(jù)加以數(shù)值處理，并將所得處理后數(shù)據(jù)保存為新的遙感影像文件的方法。

  首先，看一下本文的具體需求。我們現(xiàn)有一個(gè)文件夾，其中含有大量.tif格式的遙感影像文件；其中，這些遙感影像文件均含有4個(gè)波段，每1個(gè)波段都表示其各自的反射率數(shù)值。而對(duì)于這些遙感影像文件，有的文件其各波段數(shù)值已經(jīng)處于0至1的區(qū)間內(nèi)（也就是反射率數(shù)據(jù)的正常數(shù)值區(qū)間），而有的文件其各波段數(shù)值則是還沒(méi)有乘上縮放系數(shù)的（在本文中，縮放系數(shù)是0.0001）。

  例如，如下圖所示，即為文件夾中某一景遙感影像?？梢钥吹狡涓鞑ǘ螖?shù)值都是大于1的，這是因?yàn)槠鋽?shù)值都是還沒(méi)有乘上縮放系數(shù)的，即是真實(shí)的反射率數(shù)值的10000倍。

  我們希望實(shí)現(xiàn)的是，對(duì)于這些遙感影像中，還沒(méi)有乘上縮放系數(shù)0.0001的遙感影像，將其像元值乘上這個(gè)縮放系數(shù)；而對(duì)于已經(jīng)縮放過(guò)（也就是像元數(shù)值已經(jīng)落在0至1區(qū)間內(nèi)）的遙感影像，則不加以任何處理。最后，將經(jīng)過(guò)上述操作后的所有圖像（無(wú)論是否執(zhí)行縮放）均保存至指定的輸出結(jié)果文件夾中。

  本文所需代碼如下。

# -*- coding: utf-8 -*-
"""
Created on Thu Apr 18 12:37:22 2024
@author: fkxxgis
"""
import os
from osgeo import gdal
original_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_Original"
output_folder = r"E:\04_Reconstruction\99_MODIS\new_data\GF_Rec"
for filename in os.listdir(original_folder):
    if filename.endswith('.tif'):
        dataset = gdal.Open(os.path.join(original_folder, filename), gdal.GA_ReadOnly)
        width = dataset.RasterXSize
        height = dataset.RasterYSize
        band_count = dataset.RasterCount
        driver = gdal.GetDriverByName('GTiff')
        output_dataset = driver.Create(os.path.join(output_folder, "New_" + filename), width, height, band_count, gdal.GDT_Float32)
        for band_index in range(1, band_count + 1):
            band = dataset.GetRasterBand(band_index)
            data = band.ReadAsArray()
            if band_index == 1:
                data = data.astype(float)
                data[data > 1] /= 10000
            elif band_index == 2:
                data = data.astype(float)
                data[data > 1] /= 10000
            elif band_index == 3:
                data = data.astype(float)
                data[data > 1] /= 10000
            elif band_index == 4:
                data = data.astype(float)
                data[data > 1] /= 10000
            output_band = output_dataset.GetRasterBand(band_index)
            output_band.WriteArray(data)
            output_band.FlushCache()
        output_dataset.SetGeoTransform(dataset.GetGeoTransform())
        output_dataset.SetProjection(dataset.GetProjection())
        dataset = None
        output_dataset = None

  首先，我們使用os.listdir()函數(shù)遍歷原始數(shù)據(jù)文件夾中的所有文件，并使用if語(yǔ)句篩選出以.tif結(jié)尾的文件；隨后，使用gdal.Open()函數(shù)打開原始影像數(shù)據(jù)集，并指定只讀模式；接下來(lái)，使用dataset.RasterXSize和dataset.RasterYSize獲取影像數(shù)據(jù)集的寬度和高度。

  隨后，使用dataset.RasterCount獲取波段數(shù)量，并使用gdal.GetDriverByName()創(chuàng)建輸出數(shù)據(jù)集的驅(qū)動(dòng)程序?qū)ο?；緊接著，通過(guò)Create()方法創(chuàng)建輸出數(shù)據(jù)集，并指定輸出文件的路徑、寬度、高度、波段數(shù)量和數(shù)據(jù)類型（gdal.GDT_Float32表示浮點(diǎn)型）。

  接下來(lái)，就可以開始使用循環(huán)，對(duì)每個(gè)文件的每個(gè)波段進(jìn)行處理。首先，使用dataset.GetRasterBand()方法獲取當(dāng)前波段對(duì)象，然后使用band.ReadAsArray()將波段數(shù)據(jù)讀取為數(shù)組；根據(jù)波段索引的不同，對(duì)波段數(shù)據(jù)進(jìn)行處理。在本文中，對(duì)4個(gè)波段進(jìn)行的其實(shí)是相同的處理，即將大于1的像素值除以10000。

  其次，使用output_dataset.GetRasterBand()方法獲取輸出數(shù)據(jù)集中的當(dāng)前波段對(duì)象，并使用output_band.WriteArray()方法將處理后的數(shù)據(jù)寫入輸出數(shù)據(jù)集。

  再次，使用dataset.GetGeoTransform()和dataset.GetProjection()分別獲取原始數(shù)據(jù)集的地理轉(zhuǎn)換和投影信息，并使用output_dataset.SetGeoTransform()和output_dataset.SetProjection()設(shè)置輸出數(shù)據(jù)集的地理轉(zhuǎn)換和投影信息。

  最后一步，關(guān)閉數(shù)據(jù)集對(duì)象。至此，代碼就完成了對(duì)每個(gè).tif文件的處理，并將處理后的數(shù)據(jù)保存到輸出文件夾中。

  此時(shí)，打開本文開頭展示的那1景遙感影像，可以看到其像素?cái)?shù)值已經(jīng)是乘上縮放系數(shù)之后的了，也就是落在了0至1的區(qū)間內(nèi)；如下圖所示。

到此這篇關(guān)于Python代碼將大量遙感數(shù)據(jù)的值縮放指定倍數(shù)的方法的文章就介紹到這了,更多相關(guān)Python遙感數(shù)據(jù)的值縮放指定倍數(shù)內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

最新評(píng)論