初始化 v1.0.0：增减挂钩后备资源调查数据库生成套件

2026-07-20 18:45:38 +08:00
parent a2dcd8b8da
commit bee076504b
5 changed files with 537 additions and 0 deletions
@@ -0,0 +1,395 @@
+#!/usr/bin/env python3
+"""
+永善县增减挂钩后备资源调查数据库生成 — 处理脚本
+
+执行模式:
+  extract_2019      提取2019年三调建设用地
+  extract_2025      提取2025年变更调查建设用地
+  overlay_intersect 叠加分析取交集（两年均为建设用地图斑）
+  erase_restrictions 逐层擦除限制要素
+  filter_area       面积过滤 + 影像核实清单输出
+"""
+import json, sys, os, glob, logging, traceback
+from pathlib import Path
+from datetime import datetime
+
+import geopandas as gpd
+import pandas as pd
+import numpy as np
+
+# ── 三调建设用地编码 ──────────────────────────────────
+CONSTRUCTION_CODES = {'0702': '农村宅基地', '0601': '工业用地', '0602': '采矿用地'}
+CODE_LIST = list(CONSTRUCTION_CODES.keys())
+
+# 常见地类编码字段名（按优先级）
+DLBM_FIELDS = ['DLBM', 'DLBZ', 'DLMC', 'YSDLMC', 'YSDLBZ', 'DLYBM', 'DLMC_1']
+
+
+# ══════════════════════════════════════════════════════
+#  工具函数
+# ══════════════════════════════════════════════════════
+
+def get_logger(log_path=None):
+    """获取/配置日志器"""
+    logger = logging.getLogger('ys_zg')
+    if logger.handlers:
+        return logger
+    logger.setLevel(logging.INFO)
+    fmt = logging.Formatter('%(asctime)s [%(levelname)s] %(message)s', datefmt='%Y-%m-%d %H:%M:%S')
+    sh = logging.StreamHandler(sys.stdout)
+    sh.setFormatter(fmt)
+    logger.addHandler(sh)
+    if log_path:
+        fh = logging.FileHandler(log_path, encoding='utf-8')
+        fh.setFormatter(fmt)
+        logger.addHandler(fh)
+    return logger
+
+
+def read_vector(path, layer=None):
+    """通用矢量读取"""
+    if path.endswith('.gpkg') and layer:
+        return gpd.read_file(path, layer=layer)
+    return gpd.read_file(path)
+
+
+def estimate_utm_crs(gdf):
+    """根据数据范围估算 UTM 投影（永善县一般在 EPSG:32648）"""
+    bounds = gdf.total_bounds
+    center_lon = (bounds[0] + bounds[2]) / 2
+    utm_zone = int((center_lon + 180) / 6) + 1
+    is_north = bounds[1] + bounds[3] >= 0
+    epsg = 32600 + utm_zone if is_north else 32700 + utm_zone
+    logger = get_logger()
+    logger.info(f"  估算投影: lon={center_lon:.2f}, zone={utm_zone}, EPSG:{epsg}")
+    return f"EPSG:{epsg}"
+
+
+def ensure_metric_crs(gdf, logger):
+    """确保数据在投影坐标系下以便计算面积"""
+    if gdf.crs and gdf.crs.is_geographic:
+        utm_crs = estimate_utm_crs(gdf)
+        logger.info(f"  地理坐标系 -> 投影到 {utm_crs}")
+        return gdf.to_crs(utm_crs)
+    return gdf
+
+
+def detect_and_filter_construction(gdf, logger):
+    """
+    自动检测地类编码字段并过滤出建设用地（0702/0601/0602）。
+    返回 (filtered_gdf, used_field, code_counts_dict)
+    """
+    # 1) 按已知字段名搜索
+    for field in DLBM_FIELDS:
+        if field not in gdf.columns:
+            continue
+        logger.info(f"  尝试字段: {field}")
+        try:
+            vals = gdf[field].astype(str).str.strip()
+            mask = vals.str[:4].isin(CODE_LIST)
+            count = mask.sum()
+            if count > 0:
+                result = gdf[mask].copy()
+                result['_code_field'] = field
+                result['_land_code'] = vals[mask].str[:4].values
+                counts = result['_land_code'].value_counts().to_dict()
+                return result, field, counts
+        except Exception:
+            continue
+
+    # 2) 遍历所有字符串字段模糊匹配
+    for field in gdf.select_dtypes(include='object').columns:
+        if field in DLBM_FIELDS:
+            continue  # 已查过
+        try:
+            vals = gdf[field].astype(str).str.strip()
+            mask = vals.str[:4].isin(CODE_LIST)
+            count = mask.sum()
+            if count > 0:
+                result = gdf[mask].copy()
+                result['_code_field'] = field
+                result['_land_code'] = vals[mask].str[:4].values
+                counts = result['_land_code'].value_counts().to_dict()
+                logger.info(f"  模糊匹配字段: {field}, 命中 {count} 条")
+                return result, field, counts
+        except Exception:
+            continue
+
+    return None, None, {}
+
+
+def calc_area(gdf, logger):
+    """计算面积（平方米）"""
+    if 'area_m2' in gdf.columns and gdf['area_m2'].notna().any():
+        return gdf  # 已有面积字段
+    gdf = ensure_metric_crs(gdf, logger)
+    gdf['area_m2'] = gdf.geometry.area.round(2)
+    return gdf
+
+
+# ══════════════════════════════════════════════════════
+#  步骤函数
+# ══════════════════════════════════════════════════════
+
+def step_extract(db_path, layer, output_path, log_path, label):
+    """提取建设用地（step1 / step2 共用）"""
+    logger = get_logger(log_path)
+    logger.info(f"{'='*60}")
+    logger.info(f"提取{label}: {db_path}")
+    logger.info(f"{'='*60}")
+
+    gdf = read_vector(db_path, layer=layer)
+    logger.info(f"  数据读取完成: {len(gdf)} 要素")
+
+    result, field, counts = detect_and_filter_construction(gdf, logger)
+    if result is None or len(result) == 0:
+        logger.error(f"❌ 未能从数据中识别建设用地编码 {CODE_LIST}，请检查字段名")
+        # 输出空 GeoPackage 避免下游步骤报文件缺失
+        empty = gpd.GeoDataFrame(columns=['geometry'], geometry='geometry', crs=gdf.crs)
+        empty.to_file(output_path, driver='GPKG')
+        sys.exit(1)
+
+    # 记录分类统计
+    for code, cnt in counts.items():
+        name = CONSTRUCTION_CODES.get(code, '未知')
+        logger.info(f"  ✅ {code} {name}: {cnt} 条")
+
+    result = calc_area(result, logger)
+    total_area = result['area_m2'].sum()
+    logger.info(f"  合计: {len(result)} 条, 总面积: {total_area:,.2f} m²")
+
+    result.to_file(output_path, driver='GPKG')
+    logger.info(f"  ✅ 结果输出 -> {output_path}")
+    return output_path
+
+
+def step_overlay(input1, input2, output_path, log_path):
+    """叠加分析取交集"""
+    logger = get_logger(log_path)
+    logger.info(f"{'='*60}")
+    logger.info("叠加分析：取 2019 ∩ 2025 交集")
+    logger.info(f"{'='*60}")
+
+    gdf1 = gpd.read_file(input1)
+    gdf2 = gpd.read_file(input2)
+    logger.info(f"  输入1: {len(gdf1)} 要素, CRS={gdf1.crs}")
+    logger.info(f"  输入2: {len(gdf2)} 要素, CRS={gdf2.crs}")
+
+    # 统一 CRS
+    if gdf1.crs != gdf2.crs:
+        gdf2 = gdf2.to_crs(gdf1.crs)
+
+    # 交集
+    result = gpd.overlay(gdf1, gdf2, how='intersection', keep_geom_type=True)
+    result = result.reset_index(drop=True)
+    result = calc_area(result, logger)
+
+    logger.info(f"  交集结果: {len(result)} 条, 总面积: {result['area_m2'].sum():,.2f} m²")
+    result.to_file(output_path, driver='GPKG')
+    logger.info(f"  ✅ 结果输出 -> {output_path}")
+    return output_path
+
+
+def step_erase_restrictions(input_path, restriction_dir, output_path, log_path):
+    """逐层擦除限制要素"""
+    logger = get_logger(log_path)
+    logger.info(f"{'='*60}")
+    logger.info("逐层擦除限制要素")
+    logger.info(f"{'='*60}")
+
+    if not os.path.isdir(restriction_dir):
+        logger.error(f"限制要素目录不存在: {restriction_dir}")
+        sys.exit(1)
+
+    gdf = gpd.read_file(input_path)
+    logger.info(f"  输入要素: {len(gdf)} 条")
+    gdf = calc_area(gdf, logger)
+
+    # 扫描目录下的矢量文件
+    restriction_files = []
+    for ext in ['*.shp', '*.gpkg', '*.geojson', '*.json']:
+        restriction_files.extend(sorted(glob.glob(os.path.join(restriction_dir, ext))))
+
+    if not restriction_files:
+        logger.warning(f"  ⚠️ 未在 {restriction_dir} 中找到矢量文件，跳过擦除")
+        gdf.to_file(output_path, driver='GPKG')
+        return output_path
+
+    logger.info(f"  发现 {len(restriction_files)} 个限制要素文件:")
+    for rf in restriction_files:
+        logger.info(f"    - {os.path.basename(rf)}")
+
+    result = gdf.copy()
+    result = ensure_metric_crs(result, logger)
+
+    for rf in restriction_files:
+        basename = os.path.basename(rf)
+        try:
+            restrict = gpd.read_file(rf)
+            if restrict.crs != result.crs:
+                restrict = restrict.to_crs(result.crs)
+
+            before = len(result)
+            before_area = result['area_m2'].sum()
+            result = gpd.overlay(result, restrict, how='difference', keep_geom_type=True)
+            result = result.reset_index(drop=True)
+            result = calc_area(result, logger)
+            removed = before - len(result)
+            removed_area = before_area - result['area_m2'].sum()
+
+            logger.info(f"    ✅ 擦除 [{basename}]: 移除 {removed} 条({removed_area:,.0f}m²), "
+                        f"剩余 {len(result)} 条({result['area_m2'].sum():,.0f}m²)")
+        except Exception as e:
+            logger.warning(f"    ⚠️ 擦除 [{basename}] 出错: {e}")
+            continue
+
+    result.to_file(output_path, driver='GPKG')
+    logger.info(f"  ✅ 擦除完成: {len(result)} 条, 总面积: {result['area_m2'].sum():,.2f} m²")
+    return output_path
+
+
+def step_filter(input_path, output_path, log_path, min_area=70.0, orthophoto_dir=None):
+    """面积过滤 + 生成影像核实清单"""
+    logger = get_logger(log_path)
+    logger.info(f"{'='*60}")
+    logger.info(f"面积过滤: >= {min_area} 平方米")
+    logger.info(f"{'='*60}")
+
+    gdf = gpd.read_file(input_path)
+    before = len(gdf)
+    gdf = calc_area(gdf, logger)
+
+    result = gdf[gdf['area_m2'] >= min_area].copy().reset_index(drop=True)
+    removed = before - len(result)
+
+    logger.info(f"  过滤前: {before} 条 | 过滤后: {len(result)} 条 | "
+                f"移除 < {min_area}m²: {removed} 条")
+
+    # 添加元数据字段
+    result['fid'] = range(1, len(result) + 1)
+    result['bdlx'] = result.get('_land_code', '')
+    result['mj'] = result['area_m2'].round(2)
+    now_str = datetime.now().strftime('%Y-%m-%d %H:%M:%S')
+    result['scrq'] = now_str
+
+    # 输出最终成果 GeoPackage
+    result.to_file(output_path, driver='GPKG', layer='final_houbei_ziyuan')
+    logger.info(f"  ✅ 最终结果输出 -> {output_path}")
+
+    # ── 输出影像核实 Excel 清单 ──
+    excel_path = output_path.replace('.gpkg', '_影像核实清单.xlsx')
+    if excel_path == output_path:
+        excel_path = os.path.join(os.path.dirname(output_path), '影像核实清单.xlsx')
+
+    export_fields = ['fid', 'bdlx', 'mj']
+    # 保留有意义的属性字段
+    for f in result.columns:
+        if f not in export_fields and f not in ['geometry', '_code_field', '_land_code', 'area_m2', 'fid']:
+            if result[f].dtype in ['object', 'int64', 'float64']:
+                export_fields.append(f)
+
+    export_df = result[export_fields].copy() if len(result) > 0 else pd.DataFrame()
+    export_df['影像核实状态'] = ''
+    export_df['影像核实备注'] = ''
+    export_df['核实人'] = ''
+    export_df['核实日期'] = ''
+
+    # 使用 openpyxl 写 Excel
+    try:
+        export_df.to_excel(excel_path, index=False, engine='openpyxl')
+        logger.info(f"  ✅ 影像核实清单输出 -> {excel_path}")
+    except Exception as e:
+        logger.warning(f"  ⚠️ Excel 输出失败: {e}")
+
+    # 统计输出
+    area_by_type = result.groupby('bdlx')['mj'].agg(['count', 'sum']).reset_index()
+    logger.info(f"  按地类统计:")
+    for _, row in area_by_type.iterrows():
+        name = CONSTRUCTION_CODES.get(str(row['bdlx'])[:4], '其他')
+        logger.info(f"    {name}: {int(row['count'])} 条, {row['sum']:,.0f} m²")
+
+    # 提示影像核实为人工步骤
+    logger.info(f"{'!'*60}")
+    logger.info("⚠️  重要提示：影像核实为人工核查步骤")
+    logger.info(f"    请使用 {excel_path} 清单")
+    logger.info(f"    结合永善县 2023年耕地流出正射影像 和 2026年房体一体正射影像")
+    logger.info(f"    逐一核实图斑：对建筑物完整、年份新的图斑标记删除")
+    logger.info(f"    最终以核实后的调查数据库为准")
+    logger.info(f"{'!'*60}")
+
+    return output_path
+
+
+# ══════════════════════════════════════════════════════
+#  Main
+# ══════════════════════════════════════════════════════
+
+def main():
+    # 参数通过 PARAMS_FILE 环境变量或 sys.argv[1] 传入
+    params = {}
+    if 'PARAMS_FILE' in os.environ:
+        with open(os.environ['PARAMS_FILE'], 'r', encoding='utf-8') as f:
+            params = json.load(f)
+    elif len(sys.argv) > 1:
+        params = json.loads(sys.argv[1])
+    else:
+        # fallback: parse from sys.argv as key=value pairs
+        for arg in sys.argv[1:]:
+            if '=' in arg:
+                k, v = arg.split('=', 1)
+                params[k] = v
+
+    mode = params.get('mode', '')
+    output = os.path.abspath(params.get('output', '/tmp/output/result.gpkg'))
+    log_path = os.path.abspath(params.get('log', '/tmp/output/process.log'))
+
+    # 确保输出目录存在
+    os.makedirs(os.path.dirname(output), exist_ok=True)
+
+    logger = get_logger(log_path)
+
+    try:
+        if mode == 'extract':
+            db_path = os.path.abspath(params.get('db_path', ''))
+            if not db_path:
+                raise ValueError("缺少参数: db_path")
+            layer = params.get('layer', None) or None
+            label = params.get('label', '建设用地')
+            step_extract(db_path, layer, output, log_path, label)
+
+        elif mode == 'overlay':
+            input1 = os.path.abspath(params.get('input1', ''))
+            input2 = os.path.abspath(params.get('input2', ''))
+            step_overlay(input1, input2, output, log_path)
+
+        elif mode == 'erase':
+            input_path = os.path.abspath(params.get('input_path', ''))
+            restrict_dir = os.path.abspath(params.get('restriction_dir', ''))
+            step_erase_restrictions(input_path, restrict_dir, output, log_path)
+
+        elif mode == 'filter':
+            input_path = os.path.abspath(params.get('input_path', ''))
+            min_area = float(params.get('min_area', 70))
+            odir = params.get('orthophoto_dir', None)
+            odir = os.path.abspath(odir) if odir else None
+            step_filter(input_path, output, log_path, min_area, odir)
+
+        else:
+            raise ValueError(f"未知模式: {mode}，支持: extract / overlay / erase / filter")
+
+    except Exception as e:
+        logger.error(f"❌ 执行失败: {e}")
+        logger.error(traceback.format_exc())
+        sys.exit(1)
+
+    # stdout JSON 结果（平台契约）
+    print(json.dumps({
+        "status": "success",
+        "mode": mode,
+        "output": output,
+        "log": log_path
+    }, ensure_ascii=False))
+
+
+if __name__ == '__main__':
+    main()