计算异质性H值（运用arcgis和Python进行区域分析）

最近需要对ecognition分割结果进行统计分析，以此来进一步判断其分割结果中的欠分割和过分割对象，在看了一篇论文后，发现了可以用一个参数H来判断每个切割对象的异质性，由于此方法需要用到arcgis和Python来配合，因此记录下。

公式大概如下：

从中可以看出，如果需要计算出参数H，我们需要先计算出每个对象的归一化方差和归一化的莫兰指数。

在计算必须的参数前，我们需要准备的数据包括：

1.原始遥感图像

2.运用ecognition进行切割后产生的标签文件和矢量文件（shp文件）。

下面开始进行操作：

1.打开arcgis，导入矢量文件和标签文件，以及原始遥感图像。

2.调用分区统计工具，输入参数同上，计算出每个区域对应的标准差，输出格式为tif格式（也就是标签图那样的形式）。

3.将上一步中生成的标准差的图像文件路径输入到Python代码中，进行处理。（这里一定要记住生成的文件是带有空间参考系的，需要用GDAL库进行读取保存操作）。

这里我在网上找了一个人家编号的读取和保存函数，可以进行调用：

from osgeo import gdal


class IMAGE:
    # 读图像文件
    def read_img(self, filename):
        dataset = gdal.Open(filename)  # 打开文件
        im_width = dataset.RasterXSize  # 栅格矩阵的列数
        im_height = dataset.RasterYSize  # 栅格矩阵的行数
        im_bands = dataset.RasterCount  # 波段数
        im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率
        im_proj = dataset.GetProjection()  # 地图投影信息，字符串表示
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)

        del dataset

        return im_width, im_height, im_bands, im_proj, im_geotrans, im_data

    # 写GeoTiff文件
    def write_img(self, filename, im_proj, im_geotrans, im_data):

        # 判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        # 判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        driver = gdal.GetDriverByName("GTiff")
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影

        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

        del dataset

Python将标准差转为平差同时对图像进行归一化处理代码如下：

from service import IMAGE
import copy
import cv2
import numpy as np

std_var_path = 'data/stardvar.tif'
var_path = 'data/var.tif'

rw = IMAGE()
im_width, im_height, im_bands, im_proj, im_geotrans, im_data = rw.read_img(std_var_path)

var = copy.deepcopy(im_data)
rows, cols = im_data.shape

for row in range(rows):
    for col in range(cols):
        var[row][col] = im_data[row][col]**2

var = cv2.normalize(var,var,0,1,norm_type=cv2.NORM_MINMAX)            # 方差归一化


rw.write_img(var_path,im_proj,im_geotrans,var)

print('transform success!')

4.上面步骤计算出了归一化的方差数据，剩下的就是计算出归一化后的莫兰指数参数，moran指数参数是描述空间相关性的参数，进行计算的时候每个区域都要有一个指标才可以进行计算，本次记录中，我计算了每个区域的灰度的均值作为这个指标参数。

具体方法为：运用工具箱中的以“表格显示分区统计“工具，以shp文件为要素区域文件，原始遥感影像为赋值图像（软件会自动转为灰度图进行处理），选择字段为FID，保证唯一性。

5.上述操作最终会生成一个表格形式的文件，如下图所示：

我们需要将运用到这个表格中的mean参数，但是对于计算moran指数的工具来说，输入只能是shp矢量文件，因此我们需要将mean这一字段放到矢量文件中，可以在矢量文件字段表格中将两个表格进行连接操作。

以FID为连接字段（因为每个对象对应唯一的FID，方便进行确认），连接设置过程以及结果如下：

6.有了以上的带有均值字段的矢量文件，我们便可以进行局部莫兰指数的计算啦（如果是全局莫兰指数见我以前写的一篇吧），打开工具箱中的聚类和异常值分析工具，输入参数如下（记得一定要选择标准化！！！）：

最终会输出一个记录莫兰指数的tif文件，和记录归一化方差的tif图一样。

7.完成以上操作，我们所需要的数据便都准备好啦，下面需要的就是开始计算H参数，这一步骤我在Python中执行，同样，根据代码自己修改路径就行啦。

from service import IMAGE
import copy
import cv2
import numpy as np

# 通过归一化莫兰指数和归一化对象内方差计算H参数
# 后期提高效率可以通过标签图来统一修改图片
moran_path = 'data/moranout.tif'
var_path = 'data/var.tif'
H_path = 'data/H.tif'

rw = IMAGE()
moran_width, moran_height, moran_bands, moran_proj, moran_geotrans, moran_data = rw.read_img(moran_path)

var_width, var_height, var_bands, var_proj, var_geotrans, var_data = rw.read_img(var_path)

'''rows, cols = moran_data.shape
H_data = np.zeros_like(var_data)

for row in range(rows):
    for col in range(cols):
        nV = var_data[row][col]
        nLMI = moran_data[row][col]
        H_data[row][col] = (nV-nLMI)/(nV+nLMI)
        print(H_data[row][col],nV,nLMI)'''

H_data = (var_data-moran_data)/(var_data+moran_data)


rw.write_img(H_path,var_proj,var_geotrans,H_data)

print('H calculate completed!')

最终会生成一个记录H参数数值的tif标记图，达成目的。

以上便是我计算H参数的步骤，如果大家有更好的方法，希望及时告诉我。