影像的轉(zhuǎn)移矩陣用于分析土地利用類型或植被指數(shù)在不同年份間的變化。該技術(shù)通過計(jì)算原始影像中像元值在下一幅影像中的分布，幫助我們理解和量化地表特征的動(dòng)態(tài)變化。

分區(qū)像元值數(shù)量統(tǒng)計(jì)是一種常見的數(shù)據(jù)分析方法，用于統(tǒng)計(jì)特定矢量要素范圍內(nèi)的特定值的像元數(shù)量或各值的像元數(shù)量。

轉(zhuǎn)移矩陣

我使用該代碼統(tǒng)計(jì)了2005年到2015年各土地利用類型的變化

import rasterio
import numpy as np
import pandas as pd
# 讀取柵格數(shù)據(jù)
file01=r'C:\Users\HP\Desktop\ras\CLCD_2005.tif'
file02=r'C:\Users\HP\Desktop\ras\CLCD_2015.tif'

out_excel = r"C:\Users\HP\Desktop\pixel_c2.xlsx"

with rasterio.open(file01) as src:
 ras01 = src.read(1)
with rasterio.open(file02) as src:
 ras02 = src.read(1)

# 計(jì)算轉(zhuǎn)移矩陣
def calculate_transition_matrix(raster_2005, raster_2015):
 # 獲取柵格值的唯一值
 unique_2005 = np.unique(raster_2005)
 unique_2015 = np.unique(raster_2015)

 # 創(chuàng)建一個(gè)字典來存儲(chǔ)轉(zhuǎn)移矩陣
 transition_dict = {}

 # 使用 NumPy 的條件判斷來獲取轉(zhuǎn)移情況
 for value_2005 in unique_2005:
 for value_2015 in unique_2015:
 #先大循環(huán)被轉(zhuǎn)移柵格值，找到對(duì)應(yīng)值的區(qū)域，在小循環(huán)轉(zhuǎn)移柵格對(duì)應(yīng)區(qū)域的各值數(shù)量
 count = np.sum((raster_2005 == value_2005) & (raster_2015 == value_2015))
 #建立被轉(zhuǎn)移和轉(zhuǎn)移柵格的值的對(duì)應(yīng)柵格數(shù)量
 transition_dict[(value_2005, value_2015)] = count#在字典中直接疊加

 # 建立一個(gè)空DataFrame以存儲(chǔ)數(shù)據(jù)
 transition_matrix = pd.DataFrame(0, index=unique_2005, columns=unique_2015)

 # 填充矩陣
 for (i, j), count in transition_dict.items():
 transition_matrix.loc[i, j] = count

 return transition_matrix

# 計(jì)算轉(zhuǎn)移矩陣
transition_matrix = calculate_transition_matrix(ras01, ras02)

# 為了防止混淆年份，重命名行列
transition_matrix.columns = [(file01.split('_')[1].split('.')[0]) + str(c) for c in transition_matrix.columns]
transition_matrix.index = [(file01.split('_')[1].split('.')[0]) + str(i) for i in transition_matrix.index]
#結(jié)果導(dǎo)出
transition_matrix.to_excel(out_excel)</pre>

分區(qū)像元值數(shù)量統(tǒng)計(jì)

該代碼我是想知道這四年的各矢量對(duì)應(yīng)的森林變化量，所以我先對(duì)這四年的柵格影像進(jìn)行疊加，在使用矢量提取各年森林的像元數(shù)量，輸出表格

import rasterio as ras
import geopandas as gpd
import os
import glob
import numpy as np
import pandas as pd
from rasterio.mask import mask
from shapely.geometry import mapping

folder=r"F:\image_download\hebei"
shp=r"C:\Users\HP\Desktop\ddd\dd\shp\test01\pol02.shp"
# 存儲(chǔ)影像數(shù)據(jù)和年份
all_ras = []
years = []

# 獲取tif文件以對(duì)影像進(jìn)行疊加
tif_files = glob.glob(os.path.join(folder, '*.tif'))
for tif in tif_files:
 with ras.open(tif) as src:
 ras_crs = src.crs
 all_ras.append(src.read(1))  # 讀取影像數(shù)據(jù)
 years.append(int(tif.split('_')[3]))  # 從文件名提取年份
print(years)
print(ras_crs)

# 讀取矢量數(shù)據(jù)
pol = gpd.read_file(shp)
pol.set_crs(ras_crs)
# 存儲(chǔ)每個(gè)要素提取的結(jié)果
output_data = []

# 循環(huán)每個(gè)面要素
for index, row in pol.iterrows():
 geometry = row['geometry']
 # 獲取面要素的屬性名稱，以便于后續(xù)數(shù)據(jù)處理的識(shí)別
 field_name = row['PID']

 # 創(chuàng)建一個(gè)字典，包含該面要素的PID和年份對(duì)應(yīng)的像元值數(shù)量
 result_dict = {'Field_Name': field_name}

 # 循環(huán)每個(gè)影像（每年一個(gè)波段）
 for i, year in enumerate(years):#這段代碼能讓你得到索引和對(duì)應(yīng)的年份，索引正好可以作為讀取波段的索引
 # 獲取該年份影像數(shù)據(jù)
 ras_data = all_ras[i]

 # 使用矢量面要素進(jìn)行掩模，得到該面要素內(nèi)的像元值
 # 注意：這里的mask返回的是掩模區(qū)域的像元值和它的仿射變換
 with ras.open(tif_files[i]) as src:
 out_image, out_transform = mask(src, [mapping(geometry)], crop=True)
#這里為何輸入[mapping(geometry)]而不是直接輸入geometry，因?yàn)間eometry是一個(gè)幾何對(duì)象，而這里需要輸入shp這種格式的所以用mapping函數(shù)轉(zhuǎn)換一下

 # 計(jì)算掩模區(qū)域內(nèi)的像元值數(shù)量
 out_image = out_image[0]  # 由于mask返回的是三維數(shù)組，取第一個(gè)波段
 unique, counts = np.unique(out_image, return_counts=True)#統(tǒng)計(jì)唯一值的像元數(shù)量
 value_count = dict(zip(unique, counts))

 # 將該年份的像元數(shù)量添加到結(jié)果字典中
 result_dict[year] = value_count.get(2, 0)#將每年值為2的像元數(shù)量存儲(chǔ)在result_dict中，鍵為年份year

 # 將每個(gè)要素的結(jié)果添加到輸出數(shù)據(jù)列表
 output_data.append(result_dict)

# 將結(jié)果保存為DataFrame
df = pd.DataFrame(output_data)

# 輸出為Excel文件
output_excel = r"C:\Users\HP\Desktop\dd\ddd\table\pixel_counts2.xlsx"
df.to_excel(output_excel, index=False)</pre>