前言
在处理栅格数据时,尽管用同一个矢量文件裁剪栅格数据,不同数据来源的栅格行列数也会出现不一致的情况。如果忽略或解决不好,会导致后续数据处理出现意想不到的误差或错误,尤其是利用编程实现数据处理时。因此,应当首先对栅格行列不一致的数据进行匹配处理,以降低出现BUG的风险。
本文利用Python算法实现栅格数据行列匹配,通过案例说明算法的有效性。
1、思路
利用同一个矢量文件裁剪的栅格数据,尽管会出现栅格行列数不一致的情况,但行列数差别不大,而且通过观察可以看出两个数据行列数的具体差别。
可以选择一个基准数据,将其余数据进行对比与匹配。
添加或删减数据某侧指定的行数或列数,输出与基准数据空间位置相一致的匹配结果。
2、代码实现
这里遥感影像读取主要使用GDAL读取。
# -*- coding: utf-8 -*-
import sys
import numpy
def DataMacth(in_base_raster, in_mtc_raster, out_raster, rc=[0, 0, 0, 0]):
'''
:param in_base_raster:
:param in_mtc_raster:
:param out_raster:
:param rc: <up, down, left, right>
:return:
'''
# 读取基准栅格数据与待匹配栅格数据及其属性参数
r_in_mtc = ReadRaster(in_mtc_raster)
r_in_base = ReadRaster(in_base_raster)
mtcdata = r_in_mtc.data
nCols = r_in_base.nCols
nRows = r_in_base.nRows
geoTrans = r_in_base.geotrans
srs = r_in_base.srs
NDV = r_in_mtc.noDataValue
# 新建与基准栅格数据行列一致的空矩阵
data_new = numpy.zeros((nRows, nCols))
# 根据指定的行列数进行增添与删减处理
for m in range(nRows):
for n in range(nCols):
if (m < rc[0] or m > nRows - rc[1] - 1) or (n < rc[2] or n > nCols - rc[3] - 1):
data_new[m][n] = NDV
else:
data_new[m][n] = mtcdata[m - rc[0]][n - rc[2]]
# 输出栅格数据
WriteGTiffFile(out_raster, nRows, nCols, data_new, geoTrans, srs, NDV, gdal.GDT_Float32)
print("\tSave as: %s" % out_raster)
调用方法与案例:
上述算法函数中的参数之一rc是一个数组类型的参数,即指定的行数或列数,四个数组元素分别表示“up”, “down”, “left”, “right”四个方位
例如:[-1, 0, 0, 0]表示数据上方减去一行。
具体调用方式如下
if __name__ == "__main__":
rootdir = <input data direction>
in_base_raster = rootdir + os.sep + r"in_base_raster.tif"
in_mtc_raster = rootdir + os.sep + r"in_mtc_raster.tif"
out_raster = rootdir + os.sep + r"out_raster.tif"
rc=[-1, 1, 0, 0]
DataMacth(in_base_raster, in_mtc_raster, out_raster, rc=rc)
3、案例
匹配前,基准栅格100 x 125,待匹配栅格101 x 126,像元数值与空间位置均不匹配。
rc = [0, -1, 0, -1]
程序执行后,基准栅格100 x 125,待匹配栅格100 x 125,像元数值与空间位置香匹配。