最近需要對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引數的步驟，如果大家有更好的方法，希望及時告訴