添加地理数据处理脚本,使用 rasterio 和 geopandas 读取和转换 CRS,并计算统计结果
This commit is contained in:
49
main.py
Normal file
49
main.py
Normal file
@@ -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')
|
||||
Reference in New Issue
Block a user