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