05-GDAL 空间分块计算-a基本思路

一、空间分块思路

空间分块计算首先是将一幅较大的影像数据分成多个小块，然后分别对这些小块进行计算。这里先列举自己设想的分块的方法，后面再进一步探索多进程分块处理的方法。

比如一幅影像，拆分成2*2块tiles计算NDVI的值。

假如图像是400行，640列，则每块所占的行列为：

假设图像的存储顺序为：

data[bn,rows,cols]
波段数、行数、列数

那我们要获取的每一个分块实际上为：

data[:,begin_rows:end_rows;begin_cols:end_cols]

获取的实际上是各个波段的，这一块的数据。

二、计算思路

主要思路是先从行开始计算，i控制行，每一行内的块的rows是不变的；变动的是每一块的列cols。

比如拆分成3*3=9块，i=0控制的第一行，j遍历0，1，2；来控制列的增加。然后i+1=1，j遍历0，1，2。

所以i控制外面的大循环： c_rows为计算的平均每块图像的行数
beginrow：i*c_rows
endrow:(i+1)*c_rows j控制内部的小循环： c_cols为计算的平均每块图像的列数
begincol：j*c_cols
endcol:(j+1)*c_clos

另外还考虑到一个问题：
影像数据可能有奇数行列，或者除不尽的。这里考虑的解决办法是前面分块取小、取整。最后一行或者最后一列的分块，取到末尾行和列。

三、实现脚本代码

拆分成2*2块tiles计算NDVI的值。

from osgeo import gdal
import numpy as np

# Read Image data
dataset = gdal.Open("G:\mypython3\gdal\data\can_tmr.img")
cols = dataset.RasterXSize
rows = dataset.RasterYSize
proj=dataset.GetProjection()
im_bands = dataset.RasterCount

datas = dataset.ReadAsArray(0,0,cols,rows)
dimen = np.array(datas).shape
print(dimen)

# Get redband and NIR band
reddata = np.array(datas[2,:,:])
reddata.astype(np.float32)
nirdata = np.array(datas[3,:,:])
nirdata.astype(np.float32)
print(cols,rows)
print(np.array(reddata).shape)

# Set the tiling number，4=2*2
ntiles = 4
unittiles = int(math.sqrt(ntiles))
# Get the number of cols and rows for each tile
c_cols=int(cols/unittiles)
c_rows=int(rows/unittiles)

# Create a new dataset to store ndvi output array
outputNDVIdata = np.zeros([rows,cols],float)
print("输出ndvi栅格大小：")
print(np.array(outputNDVIdata).shape)

# Use parameter i to iterate in rows; use parameter j to iterate in cols
for i in range(unittiles):
    if i == (unittiles - 1):
        # if this is the last row
        for j in range(unittiles):
            if j == (unittiles - 1):
                # if this is the last col
                begin_rows = (c_rows) * i
                end_rows = rows
                begin_cols = (c_cols) * j
                end_cols = cols

                ndvi = (nirdata[begin_rows:end_rows, begin_cols:end_cols] - (
                    reddata[begin_rows:end_rows, begin_cols:end_cols])) / (
                               reddata[begin_rows:end_rows, begin_cols:end_cols] + nirdata[begin_rows:end_rows,
                                                                                   begin_cols:end_cols])
                outputNDVIdata[begin_rows:end_rows, begin_cols:end_cols] = ndvi

                print('i:' + str(i) + 'j:' + str(j))
                print(ndvi)
            else:
                # if this is the last row,but not the last col
                begin_rows = (c_rows) * i
                end_rows = rows
                begin_cols = (c_cols) * j
                end_cols = (j + 1) * c_cols

                ndvi = (nirdata[begin_rows:end_rows, begin_cols:end_cols] - (
                    reddata[begin_rows:end_rows, begin_cols:end_cols])) / (
                               reddata[begin_rows:end_rows, begin_cols:end_cols] + nirdata[begin_rows:end_rows,
                                                                                   begin_cols:end_cols])
                outputNDVIdata[begin_rows:end_rows, begin_cols:end_cols] = ndvi

                print('i:' + str(i) + 'j:' + str(j))
                print(ndvi)
    else:
        for j in range(unittiles):
            if j == (unittiles - 1):
                # if this is the last col, but not the last row
                begin_rows = (c_rows) * i
                end_rows = (i + 1) * c_rows
                begin_cols = (c_cols) * j
                end_cols = cols

                ndvi = (nirdata[begin_rows:end_rows, begin_cols:end_cols] - (
                    reddata[begin_rows:end_rows, begin_cols:end_cols])) / (
                               reddata[begin_rows:end_rows, begin_cols:end_cols] + nirdata[begin_rows:end_rows,
                                                                                   begin_cols:end_cols])
                outputNDVIdata[begin_rows:end_rows, begin_cols:end_cols] = ndvi

                print('i:' + str(i) + 'j:' + str(j))
                print(ndvi)
            else:
                # if this is not the last col,  not the last row neither
                begin_rows = (c_rows) * i
                end_rows = (i + 1) * c_rows
                begin_cols = (c_cols) * j
                end_cols = (j + 1) * c_cols

                ndvi = (nirdata[begin_rows:end_rows, begin_cols:end_cols] - (
                    reddata[begin_rows:end_rows, begin_cols:end_cols])) / (
                               reddata[begin_rows:end_rows, begin_cols:end_cols] + nirdata[begin_rows:end_rows,
                                                                                   begin_cols:end_cols])
                outputNDVIdata[begin_rows:end_rows, begin_cols:end_cols] = ndvi

                print('i:' + str(i) + 'j:' + str(j))
                print(ndvi)



print("输出的NDVI数组：")
print(np.array(outputNDVIdata).shape)
print(outputNDVIdata)

# Create output raster
driver = gdal.GetDriverByName("GTiff")
outDs = driver.Create("G:\mypython3\gdal\data\\ndvi.tif",cols,rows,1,gdal.GDT_Float32)
#dataset.SetProjection(proj) #写入投影
outDs.SetProjection(proj)
outBand = outDs.GetRasterBand(1)
outBand.WriteArray(outputNDVIdata)

# Close raster file
outDs = None
del outDs, outBand