影像的转移矩阵用于分析土地利用类型或植被指数在不同年份间的变化。该技术通过计算原始影像中像元值在下一幅影像中的分布，帮助我们理解和量化地表特征的动态变化。
分区像元值数量统计是一种常见的数据分析方法，用于统计特定矢量要素范围内的特定值的像元数量或各值的像元数量。
转移矩阵
我使用该代码统计了2005年到2015年各土地利用类型的变化
import rasterio
import numpy as np
import pandas as pd
# 读取栅格数据
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)
# 计算转移矩阵
def calculate_transition_matrix(raster_2005, raster_2015):
# 获取栅格值的唯一值
unique_2005 = np.unique(raster_2005)
unique_2015 = np.unique(raster_2015)
# 创建一个字典来存储转移矩阵
transition_dict = {}
# 使用 NumPy 的条件判断来获取转移情况
for value_2005 in unique_2005:
for value_2015 in unique_2015:
#先大循环被转移栅格值，找到对应值的区域，在小循环转移栅格对应区域的各值数量
count = np.sum((raster_2005 == value_2005) & (raster_2015 == value_2015))
#建立被转移和转移栅格的值的对应栅格数量
transition_dict[(value_2005, value_2015)] = count#在字典中直接叠加
# 建立一个空DataFrame以存储数据
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
# 计算转移矩阵
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]
#结果导出
transition_matrix.to_excel(out_excel)</pre>
分区像元值数量统计
该代码我是想知道这四年的各矢量对应的森林变化量，所以我先对这四年的栅格影像进行叠加，在使用矢量提取各年森林的像元数量，输出表格
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"
# 存储影像数据和年份
all_ras = []
years = []
# 获取tif文件以对影像进行叠加
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)) # 读取影像数据
years.append(int(tif.split('_')[3])) # 从文件名提取年份
print(years)
print(ras_crs)
# 读取矢量数据
pol = gpd.read_file(shp)
pol.set_crs(ras_crs)
# 存储每个要素提取的结果
output_data = []
# 循环每个面要素
for index, row in pol.iterrows():
geometry = row['geometry']
# 获取面要素的属性名称，以便于后续数据处理的识别
field_name = row['PID']
# 创建一个字典，包含该面要素的PID和年份对应的像元值数量
result_dict = {'Field_Name': field_name}
# 循环每个影像（每年一个波段）
for i, year in enumerate(years):#这段代码能让你得到索引和对应的年份，索引正好可以作为读取波段的索引
# 获取该年份影像数据
ras_data = all_ras[i]
# 使用矢量面要素进行掩模，得到该面要素内的像元值
# 注意：这里的mask返回的是掩模区域的像元值和它的仿射变换
with ras.open(tif_files[i]) as src:
out_image, out_transform = mask(src, [mapping(geometry)], crop=True)
#这里为何输入[mapping(geometry)]而不是直接输入geometry，因为geometry是一个几何对象，而这里需要输入shp这种格式的所以用mapping函数转换一下
# 计算掩模区域内的像元值数量
out_image = out_image[0] # 由于mask返回的是三维数组，取第一个波段
unique, counts = np.unique(out_image, return_counts=True)#统计唯一值的像元数量
value_count = dict(zip(unique, counts))
# 将该年份的像元数量添加到结果字典中
result_dict[year] = value_count.get(2, 0)#将每年值为2的像元数量存储在result_dict中，键为年份year
# 将每个要素的结果添加到输出数据列表
output_data.append(result_dict)
# 将结果保存为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>
