栅格数据行列匹配处理

Row-column matching for raster data

在处理栅格数据时，尽管用同一个矢量文件裁剪栅格数据，不同数据来源的栅格行列数也会出现不一致的情况。如果忽略或解决不好，会导致后续数据处理出现意想不到的误差或错误，尤其是利用编程实现数据处理时。因此，应当首先对栅格行列不一致的数据进行匹配处理，以降低出现BUG的风险。

本文利用Python算法实现栅格数据行列匹配，通过案例说明算法的有效性，并附算法源码。

（实现方法其实很简单~）

思路

利用同一个矢量文件裁剪的栅格数据，尽管会出现栅格行列数不一致的情况，但行列数差别不大，而且通过观察可以看出两个数据行列数的具体差别。

可以选择一个基准数据，将其余数据进行对比与匹配。

添加或删减数据某侧指定的行数或列数，输出与基准数据空间位置相一致的匹配结果。

算法实现

其中栅格数据的读取和输出程序可参考：栅格图像补偿简单处理方法，算法依赖库：GDAL in Python

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# -*- 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]表示数据上方减去一行。

具体调用方式如下

案例

匹配前，基准栅格100 x 125，待匹配栅格101 x 126，像元数值与空间位置均不匹配。

rc = [0, -1, 0, -1]

程序执行后，基准栅格100 x 125，待匹配栅格100 x 125，像元数值与空间位置香匹配。

Fighting, GISer!