歡迎關注博主的微信公眾號：“智慧遙感”。

該公眾號將為您奉上Python地學分析、爬蟲、資料分析、Web開發、機器學習、深度學習等熱門原始碼。

本人的GitHub程式碼資料主頁（持續更新中，多給Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

" 當單幅遙感影像較大時，也就是解析度較高或者像元數量較多時，如果批量處理這些影像，計算機記憶體可能不夠，程式容易報錯：記憶體溢位。這時需要對影像進行分塊讀取與處理，也是本節所要分享的重點。"

按塊讀取遙感影像

在上一節中，程式一次讀取並儲存了整個波段的資料。但是，如果單幅影象尺寸較大（行列數較大）的話，我們可以將其分解為塊來讀取。可能是因為你只需要影象中的某一塊，或者你的本本沒有足夠的記憶體來同時儲存所有資料。下面我們就來一起看看如何訪問影象的子集而不是整幅影象。

ReadAsArray函式有幾個可選引數，儘管它們根據你使用的是資料集還是波段而有所不同。這是波段版本的引數：

band.ReadAsArray([xoff], [yoff], [win_xsize], [win_ysize], [buf_xsize],[buf_ysize], [buf_obj])
 

• xoff是開始讀取的列。預設值為0。
• yoff是開始讀取的行。預設值為0。
• win_xsize是要讀取的列數。預設是全部讀取它們。
• win_ysize是要讀取的行數。預設是全部讀取它們。
• buf_xsize是輸出陣列中的列數。預設設定是使用win_xsize值。如果此值與win_xsize不同，則將重新取樣資料。
• buf_ysize是輸出陣列中的行數。預設設定是使用win_ysize值。如果此值與win_ysize不同，則將重新取樣資料。
• buf_obj是一個NumPy陣列，用於放入資料而不是建立新陣列。如果需要，將重新取樣資料以適合此陣列。值也將轉換為此陣列的資料型別。

xoff和yoff引數分別指開始讀取的列和行偏移量。預設設定是從第一行和第一列開始讀取。win_xsize和win_ysize引數是指要讀取的行數和列數，預設值是全部讀取它們。buf_xsize和buf_ysize引數允許你指定輸出陣列的大小。如果這些值與win_xsize和win_ysize值不同，則資料將在讀取時重新取樣，以匹配輸出陣列大小。buf_obj引數是一個NumPy陣列，資料將被儲存，而不是正在建立的新陣列。將更改像元資料型別以匹配此陣列的資料型別。如果提供的buf_xsize和buf_ysize值與此陣列的維度不匹配，則會出現錯誤。

例如，要讀取從圖1中所示的第6000行和第1400列開始的六列和三行，你可以執行以下操作：

data = band.ReadAsArray(1400, 6000, 6, 3)

圖1 使用ReadAsArray（1400,6000,6,3）從第6000行和第1400列開始讀取六列和三行。

如果你需要像元值是浮點型而不是位元組型，你可以在讀完之後使用NumPy轉換它們，如下所示：

data = band.ReadAsArray(1400, 6000, 6, 3).astype(float)

或者你可以讓GDAL在讀取資料時為你進行轉換。要使用此方法，先建立一個浮點陣列，然後將其作為buf_obj引數傳遞給ReadAsArray。確保建立的陣列與正在讀取的資料具有相同的尺寸。

import numpy as np

data = np.empty((3, 6), dtype=float)
band.ReadAsArray(1400, 6000, 6, 3, buf_obj=data)

NumPy空函式建立一個尚未使用任何值初始化的陣列，因此它包含垃圾，直到你以某種方式填充它。函式的第一個引數是一個元組，其中包含要建立的陣列的尺寸。如果它是二維陣列，則元組包含行數，然後包含列數。dtype引數是可選的，它指定陣列將儲存的資料型別。如果未提供，則陣列將保留浮點數。

要將資料陣列寫入其他資料集中的特定位置，請將偏移量傳遞給WriteArray。它會寫出你傳遞給函式的陣列中的所有資料，從你提供的偏移量開始寫。

band2.WriteArray(data, 1400, 6000)

閱讀部分資料集時要記住一件重要的事情，那就是你必須確保讀取的子資料尺寸不要大於父資料尺寸，否則程式會報錯。例如，如果影象有100行，並且你要求它以偏移75開始讀取並讀取30行，那麼它將超過影象的末尾並將報錯。

你如何使用此資訊處理不適合RAM的大型資料集？一種方法是一次處理一個塊。請記住，柵格將資料儲存在磁碟中。因為塊中的資料一起儲存在磁碟上，所以處理這些塊中的影象是有效的。

以下程式碼顯示瞭如何將數字高程模型從米轉換為英尺，一次一個塊。這是一個小資料集，因此你可能不需要將其拆分，但是你將以相同的方式處理大型資料集。這也向你展示了在柵格中處理NoData值的示例，這些像元被認為具有空值。像元必須具有值，但是特定值可以指定為NoData，因此會被忽略。

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/12/6 17:19'

import os
import numpy as np
from osgeo import gdal

os.chdir(r'D:\osgeopy-data\Washington\dem')

in_ds = gdal.Open('gt30w140n90.tif')
in_band = in_ds.GetRasterBand(1)
xsize = in_band.XSize
ysize = in_band.YSize
block_xsize, block_ysize = in_band.GetBlockSize()
nodata = in_band.GetNoDataValue()

out_ds = in_ds.GetDriver().Create(
    'dem_feet.tif', xsize, ysize, 1, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)

for x in range(0, xsize, block_ysize):
    if x + block_xsize < xsize:
        cols = block_xsize
    else:
        cols = xsize - x
    for y in range(0, ysize, block_ysize):
        if y + block_ysize < ysize:
            rows = block_ysize
        else:
            rows = ysize - y
        data = in_band.ReadAsArray(x, y, cols, rows)
        # 將米轉換為英尺
        data = np.where(data == nodata, nodata, data * 3.28084)
        out_band.WriteArray(data, x, y)

out_band.FlushCache()
out_band.SetNoDataValue(nodata)
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
del out_ds

程式一開始，開啟資料集並獲取有關波段的資訊，包括塊的大小和NoData值。建立輸出資料集後，開始在水平(x)方向上迴圈遍歷塊。從第0列開始，然後轉到最後一列，用索引xsize表示。執行每次迴圈時，你將x增加一個塊中的列數（範圍的第三個引數是要增加的量），所以你從一個塊的開頭跳到下一個塊的開頭。然後儲存要在名為cols的變數中讀取的列數。如果剩下要讀取的完整塊的列數，則將此變數設定為塊中的列數。但是如果沒有足夠的列，如圖1中x等於10的情況那樣，則使用剩餘列的數量（圖中為3）。你需要這樣做，因為如果你嘗試讀取比已有資料更多的行或列，你的程式將報錯。

計算要讀取的列數後，重複執行行數。如圖1所示，通過第二個迴圈的前兩次讀取五行，但第三次只剩下一行讀取。處理完所有行的前五列後，進入外迴圈的下一次迭代並處理接下來的五列，然後處理最後三列。

data = in_band.ReadAsArray(x, y, cols, rows)

處理完所有塊後，你可以計算統計資料並構建概檢視。要從統計計算中排除NoData像元，你必須在呼叫ComputeStatistics之前告訴波段哪個值代表NoData。你可能想要計算迴圈內的統計資訊，但是你希望統計資訊基於波段中的所有像元，因此你需要等到所有像元值都已計算完畢。