#!/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) # 确保坐标系一致 target_crs = province.crs or "EPSG:4326" province = province.to_crs(target_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(target_crs) # 空间判断:城市中心点是否在省界内 city_centroid = city.geometry.centroid.iloc[0] in_province = province.geometry.contains(city_centroid).any() # 也检查 intersect 作为辅助判断 intersects = province.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()