计算异质性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参数的步骤，如果大家有更好的方法，希望及时告诉我。
