From f9edc3ea134cd11a8ee110d789c220bf2e44998d Mon Sep 17 00:00:00 2001
From: missum <missum.cn@gmail.com>
Date: Sun, 9 Nov 2025 03:58:05 +0800
Subject: [PATCH] =?UTF-8?q?=E6=B7=BB=E5=8A=A0=E5=9C=B0=E7=90=86=E6=95=B0?=
 =?UTF-8?q?=E6=8D=AE=E5=A4=84=E7=90=86=E8=84=9A=E6=9C=AC=EF=BC=8C=E4=BD=BF?=
 =?UTF-8?q?=E7=94=A8=20rasterio=20=E5=92=8C=20geopandas=20=E8=AF=BB?=
 =?UTF-8?q?=E5=8F=96=E5=92=8C=E8=BD=AC=E6=8D=A2=20CRS=EF=BC=8C=E5=B9=B6?=
 =?UTF-8?q?=E8=AE=A1=E7=AE=97=E7=BB=9F=E8=AE=A1=E7=BB=93=E6=9E=9C?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 main.py | 49 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 main.py

diff --git a/main.py b/main.py
new file mode 100644
index 0000000..fd0ad0a
--- /dev/null
+++ b/main.py
@@ -0,0 +1,49 @@
+import os
+import sys
+import geopandas as gpd
+import rasterio
+from exactextract import exact_extract
+
+try:
+    proj_lib_path = os.path.join(sys.prefix, 'Lib', 'site-packages', 'rasterio', 'proj_data')
+    os.environ['PROJ_LIB'] = proj_lib_path
+except Exception as e:
+    print(f"Warning: Could not automatically set PROJ_LIB. Please set it manually. Error: {e}")
+
+# 定义文件路径
+raster_path = "D:/工作/三普成果编制/出图数据/容县/栅格0925/AB/AB.tif"
+vector_path = "D:/测试文件夹/容县耕园林草.shp"
+
+# 1. 使用 rasterio 读取栅格的坐标参考系统 (CRS)
+with rasterio.open(raster_path) as src:
+    raster_crs = src.crs
+
+# 2. 使用 geopandas 读取矢量文件
+gdf = gpd.read_file(vector_path)
+
+# 3. 检查并转换矢量数据的 CRS
+print(f"原始矢量CRS: {gdf.crs}")
+print(f"目标栅格CRS: {raster_crs}")
+
+if gdf.crs != raster_crs:
+    print("CRS不匹配，正在转换矢量数据的CRS...")
+    # 使用 .to_crs() 方法进行转换
+    gdf = gdf.to_crs(raster_crs)
+    print("转换完成。")
+
+# 4. 现在，将已经对齐了坐标系的 GeoDataFrame 传递给 exact_extract
+#    注意：可以直接传递 GeoDataFrame 对象，而不仅仅是文件路径
+stats_to_calculate = ['mean', 'sum', 'count', 'min', 'max']
+results = exact_extract(raster_path, gdf, stats_to_calculate)
+
+# 5. 将结果合并回 GeoDataFrame
+#    exact_extract 在处理 GeoDataFrame 时，会保留原始的行顺序
+for stat in stats_to_calculate:
+    # 从结果列表中提取每个要素的'properties'字典中的统计值
+    gdf[stat] = [res['properties'][stat] for res in results]
+
+# 打印最终带有统计结果的 GeoDataFrame
+print("\n分区统计结果:")
+print(gdf.head())
+
+gdf.to_file("ddddl.shp", driver='ESRI Shapefile', encoding='utf-8')
\ No newline at end of file