SuiteForge
103 lines
3.4 KiB
Python
103 lines
3.4 KiB
Python
#!/usr/bin/env python3
|
|
"""
|
|
Step 2: 叠加分析
|
|
判断各城市是否在省界内（spatial join），输出分析结果
|
|
"""
|
|
import sys
|
|
import json
|
|
import os
|
|
import glob
|
|
import geopandas as gpd
|
|
|
|
|
|
def main():
|
|
province_shp = os.environ.get("PARAM_PROVINCE_SHP", "")
|
|
cities_dir = os.environ.get("PARAM_CITIES_DIR", "/tmp/output/shp")
|
|
output_dir = os.environ.get("PARAM_OUTPUT_DIR", "/tmp/output/result")
|
|
|
|
if not province_shp or not os.path.isfile(province_shp):
|
|
print(json.dumps({"error": f"省界 SHP 文件不存在: {province_shp}", "success": False}))
|
|
sys.exit(1)
|
|
|
|
os.makedirs(output_dir, exist_ok=True)
|
|
|
|
# 读取省界
|
|
print(json.dumps({"msg": f"读取省界 SHP: {province_shp}"}))
|
|
province = gpd.read_file(province_shp)
|
|
|
|
if province.empty:
|
|
print(json.dumps({"error": "省界数据为空", "success": False}))
|
|
sys.exit(1)
|
|
|
|
# 确保坐标系一致
|
|
source_crs = province.crs or "EPSG:4326"
|
|
# 使用投影坐标系（墨卡托）计算中心点以避免 geographic CRS 警告
|
|
projected_crs = "EPSG:3857"
|
|
province_proj = province.to_crs(projected_crs)
|
|
|
|
# 查找所有城市 SHP
|
|
city_shp_files = glob.glob(os.path.join(cities_dir, "*.shp"))
|
|
city_shp_files = [f for f in city_shp_files
|
|
if os.path.basename(f) != "province.shp"]
|
|
|
|
if not city_shp_files:
|
|
print(json.dumps({"error": "未找到城市 SHP 文件", "success": False}))
|
|
sys.exit(1)
|
|
|
|
results = []
|
|
|
|
for shp in city_shp_files:
|
|
city_name = os.path.splitext(os.path.basename(shp))[0]
|
|
print(json.dumps({"msg": f"分析城市: {city_name}"}))
|
|
city = gpd.read_file(shp)
|
|
|
|
if city.empty:
|
|
results.append({"city": city_name, "in_province": False,
|
|
"reason": "数据为空", "status": "跳过"})
|
|
continue
|
|
|
|
city = city.to_crs(projected_crs)
|
|
|
|
# 空间判断：城市中心点是否在省界内（使用投影坐标系避免 CRS 警告）
|
|
city_centroid = city.geometry.centroid.iloc[0]
|
|
in_province = province_proj.geometry.contains(city_centroid).any()
|
|
|
|
# 也检查 intersect 作为辅助判断
|
|
intersects = province_proj.geometry.intersects(city.geometry).any()
|
|
|
|
status = "通过 ✅" if in_province else "不通过 ❌"
|
|
|
|
city_name_prop = city.iloc[0].get("name", city_name)
|
|
city_code = city.iloc[0].get("code", "")
|
|
city_pop = city.iloc[0].get("population", 0)
|
|
|
|
results.append({
|
|
"city": city_name,
|
|
"city_code": str(city_code),
|
|
"city_name": str(city_name_prop),
|
|
"population": int(city_pop) if city_pop else 0,
|
|
"in_province": bool(in_province),
|
|
"intersects": bool(intersects),
|
|
"status": status
|
|
})
|
|
print(json.dumps({"msg": f"{city_name}: {status}"}))
|
|
|
|
# 输出分析结果 JSON
|
|
result_path = os.path.join(output_dir, "overlay_results.json")
|
|
with open(result_path, "w", encoding="utf-8") as f:
|
|
json.dump(results, f, ensure_ascii=False, indent=2)
|
|
|
|
summary = {
|
|
"success": True,
|
|
"total_cities": len(results),
|
|
"in_province": sum(1 for r in results if r["in_province"]),
|
|
"not_in_province": sum(1 for r in results if not r["in_province"]),
|
|
"results": results,
|
|
"result_json": result_path
|
|
}
|
|
print(json.dumps(summary))
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|