計算異質性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參數的步驟，如果大家有更好的方法，希望及時告訴我。