批处理---以表格显示分区统计

在影像处理中，经常需要使用arcgis以表格显示分区统计工具以根据矢量对栅格进行空间统计，就是统计对应矢量要素内像元的最大值、均值以及最小值等。但是如果栅格数量比较多，使用软件就有点麻烦，如下代码循环对各栅格进行区域统计，并将输出结果统一保存到Excel表中以方便我们处理。此外我还添加一个以点矢量数据提取对应栅格值的数据，并将坐标点和对应的值保存到Excel表中。

rasterstats（https://pythonhosted.org/rasterstats/manual.html#zonal-statistics）是一个用于区域统计的包，能实现面、点矢量对栅格数据的提取。

以下两个库rasterio、geopandas分别是常用于处理栅格和矢量的库，以下链接为其使用指南：
rasterio user guide：https://rasterio.readthedocs.io/en/stable/
geopandas user guide：https://geopandas.org/en/stable/docs/user_guide.html

以表格显示分区统计

统计各矢量要素内像素的最大值、最小值、均值等，并将结果输出保存至表格

import geopandas as gpd
import rasterio
from rasterstats import zonal_stats
import pandas as pd
import os
import glob
#读取文件
vector_path = r"C:\Users\HP\Desktop\test011\shi.shp"
gdf = gpd.read_file(vector_path)

folder_path = r'E:\project_bfu\ddd\ddd\cc\zz\Night01'
output_excel = r"C:\Users\HP\Desktop\test011\sum03.xlsx"
tif_files = glob.glob(os.path.join(folder_path, '*.tif'))
all_results = []
# 循环读取每个.tif文件
for tif_file in tif_files:

    with rasterio.open(tif_file) as src:
        affine = src.transform#获取栅格图像的仿射变换矩阵
        array = src.read(1)  # 读取栅格的第一波段

    # 计算分区统计,affine确保栅格数据的像素能够正确映射到矢量区域的坐标系统,也就是保证两者坐标系一致
    stats = zonal_stats(vector_path, tif_file, stats=["mean", "min", "max", "std"], affine=affine)

    # 结果整理为DataFrame
    stats_df = pd.DataFrame(stats)
    result = gdf.join(stats_df)#统计结果与矢量属性字段相连

    # 生成最终结果，仅保留相关字段，并增加影像名称作为新的一列
    result_filtered = result[["shi_mc", "mean", "min", "max", "std"]].copy()
    result_filtered.insert(0, "TIF_Name", os.path.basename(tif_file))  # 在最前面插入影像名称

    # 追加到总列表
    all_results.append(result_filtered)
#注意这里all_results还是一个列表，但他的格式是[df1,df2,df3],列表里的每个元素都是一个 DataFrame，所以下文中你可以直接使用concat()函数合并

#将所有结果合并并输出到 Excel

final_df = pd.concat(all_results,ignore_index=True)
print(all_results)

#写入 Excel
final_df.to_excel(output_excel, index=False, sheet_name="NightLightStats")

点数据提取对应位置栅格的值（值提取至点）

在点要素类的指定位置提取栅格像元值，并将这些值输出保存至表格

import geopandas as gpd
from rasterstats import zonal_stats, point_query
import pandas as pd

#读取文件
vector_path = r"C:\Users\HP\Desktop\test011\dd.shp"
tif_file=r"C:\Users\HP\Desktop\test011\night2014.tif"

#读取文件的坐标并转为DataFrame格式
gdf = gpd.read_file(vector_path)
gdf['X'] = gdf.geometry.x  # 获取经度
gdf['Y'] = gdf.geometry.y   # 获取纬度
df = gdf[['X', 'Y']]
df = pd.DataFrame(df)

#读取对应坐标点的栅格值
pts = point_query(vector_path , tif_file)
value=pd.DataFrame(pts,columns=["value"])
sum=pd.concat([df,value],axis=1)#将坐标位置与栅格值连接输出
print(sum)