基于输入矢量边界自动分块处理栅格影像算法（不产生中间文件）

TOC](基于输入矢量边界自动分块处理栅格影像)

情景分析

在进行大区域栅格影像处理时，会存在数据量太大，导致进行数据分析处理时，存在内存不足的问题，针对于此现象可以通过规则机械分块（指定规则的行列数进行分块处理）或者通过导入自定义的矢量边界进行分块处理，本算法主要用于侵蚀沟提取中，由于侵蚀沟分布存在明显的流域分异特性，因此采用的提取的小流域矢量边界进行分块处理。接下来主要介绍基于自定义矢量边界进行影像分块处理的方法。

数据准备

自定义矢量边界数据

小流域划分矢量边界，包含多个要素的矢量数据。按每一个要素对象进行分块数据处理。

要处理的栅格影像数据

小流域划分矢量边界，包含多个要素的矢量数据。按每一个要素对象进行分块数据处理。

解决思路

基于GDAL进行影像数据时，其实是对数组进行计算。在对单波段影像进行处理时，就是对二维数组进行处理，在进行大区域影像处理时，有些算法比较复杂，导致出现内存不足的问题，此时，可先进行一个分块处理，通过设置掩膜数组的方式对读取的影像数组进行一个分区处理，具体步骤如下：
第一步：循环读取矢量数据中的每个要素
第二步：读取要素矢量边界，建立外接多边形，并计算读取影像的范围（左上角位置、行数和列数），并读取外界多边形内的影像
第三步：根据每个要素的矢量边界，建立掩膜数组对第二步读取的数组进行掩膜
第四步：对掩膜后的数组进行处理分析。

下图中红色边框为黄色要素的外接矩形

代码

下面展示一些 内联代码片。

// 导入包
import numpy as np
import numpy.ma as ma
from rasterio import features
from shapely.geometry import Polygon
from affine import Affine
from osgeo import gdal
//用于读写shp数据，是gdal包的一部分
from osgeo import ogr 
//读取shp数据
shp_test = ogr.Open(r"I:\hhhh\watershed40000.shp")
// 读取栅格数据
raster = gdal.Open(r'E:\gaofen\meanhill.tif') 
// 读取栅格数据的列数
cols = raster.RasterXSize
// 读取栅格数据的行数
rows = raster.RasterYSize
// 将栅格数据读取为数组
meanhill_array = raster.ReadAsArray(0,0,cols,rows)
//获取栅格影像的地理坐标信息（仿射变换参数）
transformgdal = raster.GetGeoTransform() 
//遥感图像的水平空间分辨率
pixelWidth = transformgdal[1] 
//遥感图像的垂直空间分辨率
pixelHeight = transformgdal[5]
//读取矢量的数据层
layer = shp_test.GetLayer() 
//得到矢量图层中有多少个要素
n = layer.GetFeatureCount() 
for i in range(2147,2148):   
//读取第i个要素数据层，以FID=2147对象为例
    feat=layer.GetFeature(i)  
    //读取第i个要素的几何形状
    geom = feat.GetGeometryRef()   
    //外接矩形边界
    polygonextent = geom.GetEnvelope()  
    //计算外接矩形最左边在第几列
    column0=int(math.floor((polygonextent[0]-transformgdal[0])/abs(pixelWidth)))
    //计算外接矩形最右边在第几列
    column1 = int(math.ceil((polygonextent[1]-transformgdal[0]) / pixelWidth))  
    //计算外接矩形在栅格影像中占多少列
    column=column1-column0 
    //计算外接矩形最上边在第几行
    row0=int(math.floor((transformgdal[3]-polygonextent[3]) / abs(pixelHeight)))   
	//计算外接矩形最下边在第几行
    row1 = int(math.ceil((transformgdal[3] - polygonextent[2]) / abs(pixelHeight)))
    //计算外接矩形在栅格影像中占多少行
    row=row1-row0
    //将外接矩形边界内的影像读取为数组
    meanhill_arrayi = raster.ReadAsArray(column0, row0, column, row)
    //将ogr.Geometry转换为shapely.Geometry.Polygon
    geom = [feat.GetGeometryRef()]
    geom_dict = eval(geom[0].ExportToJson())
    geom1 = [Polygon(geom_dict['coordinates'][0]).__geo_interface__]
    //读取外接矩形的最小x，y坐标
    xmin=polygonextent[0]
    ymax=polygonextent[3]
    //rasterio不支持gdal形式的transform，转换为Affine
    transform = Affine.from_gdal(xmin, pixelWidth, 0, ymax, 0, pixelHeight)
    //建立掩膜数组
    datamask=features.geometry_mask(geom1,[row,column], transform)
    //利用掩膜数组，对外接矩形内的数组进行掩膜
    clip_mask = ma.MaskedArray(meanhill_arrayi, datamask)
    //计算掩膜之后数组的平均值
    mean_clip_mask = np.nanmean(clip_mask) 
    //输出掩膜之后数组的行列数
    print(clip_mask.shape) 
    //输出掩膜之后数组的平均值
    print(mean_clip_mask)
    mean_meanhill_arrayi = np.nanmean(meanhill_arrayi) #计算读取外接矩形内数组的平均值
    //输出外接矩形内数组的行列数
    print(meanhill_arrayi.shape) 
    //输出外接矩形内数组的平均值
    print(mean_meanhill_arrayi)
mean_meanhill_array = np.nanmean(meanhill_array) #计算整个影像数组的平均值
//输出整个影像数组的行列数
print(meanhill_array.shape) 
//输出整个影像数组的平均值
print(mean_meanhill_array)
---输出结果---
(340, 219)    #掩膜数组的行列数
148.77304     #掩膜数组的平均值
(340, 219)    #外接矩形的行列数
151.44310    #外接矩形的平均值
D:\Software\Conda\envs\python37\lib\site-packages\numpy\core\fromnumeric.py:86: RuntimeWarning: overflow encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)     #计算整个影像数组平均值出现警告
(22597, 26500) #整个影像数组的行列数
-inf   #由于整个影像数组太大，平均值计算失败

为了检验其计算是否正确，分别在arcmap中查看整个影像、外接矩形影像以及掩膜影像的统计值，与计算结果进行对比一致

 整个影像行列数

)

   外接矩形影像行列数

  外接矩形影像统计平均值

掩膜影像行列数

掩膜影像统计平均值

适用场景分析

（1）在进行大区域影像数据处理时，数据处理速度慢需要分块处理
（2）在进行影像数据处理时，需要分区设定阈值或者分区设定算法
（3）除外之外，还可以拓展其应用，思路类似。现在许多机器学习算法对影像进行处理时，一般需要划定512*512大小的分块影像，可以通过构建渔网矢量进行循环读取，也可通过写一个递增数列，依次改变读取影像的初始左上角坐标。
（4）本方法中通过直接定位读取外接边界，可以大大减少数据处理运算时间，若直接通过对整个影像进行直接掩膜需要2817.463159561157毫秒，通过定位直接读取外接边界再掩膜仅需要30.91740608215332毫秒，速度提高91倍。

参考资料

[1] https://zhuanlan.zhihu.com/p/131341925
[2] https://www.osgeo.cn/python-gdal-utah-tutorial/ch02.html
[3] https://www.jianshu.com/p/f459d637cc47
[4] https://www.osgeo.cn/pygis/shapely-intro.html
[5] https://vimsky.com/zh-tw/examples/detail/python-method-affine.Affine.from_gdal.html
[6]https://blog.csdn.net/Prince999999/article/details/105650448
[7] https://programtalk.com/python-examples/rasterio.features.geometry_mask/
[8]https://cloud.tencent.com/developer/article/1618306
[9] https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html
[10] https://numpy.org.cn/reference/arrays/maskedarray.html#numpy-ma%E6%A8%A1%E5%9D%97
[11]https://blog.csdn.net/weixin_41712499/article/details/85208928