diff --git a/scripts/__pycache__/process.cpython-312.pyc b/scripts/__pycache__/process.cpython-312.pyc new file mode 100644 index 0000000..f4ce522 Binary files /dev/null and b/scripts/__pycache__/process.cpython-312.pyc differ diff --git a/scripts/process.py b/scripts/process.py new file mode 100644 index 0000000..17dcda4 --- /dev/null +++ b/scripts/process.py @@ -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() diff --git a/setup_repo.sh b/setup_repo.sh new file mode 100644 index 0000000..bf51bcc --- /dev/null +++ b/setup_repo.sh @@ -0,0 +1,37 @@ +#!/bin/bash +# 在 SuiteHub 组织下创建仓库并推送代码 +set -e + +TOKEN="135e5d…4cb6" +ORG="SuiteHub" +REPO="yongshan-zengjian-guagou" +GITEA="https://git.mercator.cn" + +echo "=== 1. 检查仓库 ===" +RESULT=$(curl -s -H "Authorization: token $TOKEN" "$GITEA/api/v1/repos/$ORG/$REPO") +if echo "$RESULT" | grep -q '"id"'; then + echo "仓库已存在" +else + echo "=== 2. 创建仓库 ===" + curl -s -X POST "$GITEA/api/v1/orgs/$ORG/repos" \ + -H "Authorization: token $TOKEN" \ + -H "Content-Type: application/json" \ + -d "{\"name\":\"$REPO\",\"description\":\"永善县增减挂钩后备资源调查数据库生成套件\",\"private\":false,\"auto_init\":true}" +fi + +echo "" +echo "=== 3. 克隆并推送代码 ===" +TMPDIR="/tmp/$REPO" +rm -rf "$TMPDIR" +git clone "https://xuhengjie:$TOKEN@$GITEA/$ORG/$REPO.git" "$TMPDIR" +cp -r /home/openclaw/.openclaw/workspace/$REPO/workflow.yaml "$TMPDIR/" +mkdir -p "$TMPDIR/scripts" +cp -r /home/openclaw/.openclaw/workspace/$REPO/scripts/*.py "$TMPDIR/scripts/" +cd "$TMPDIR" +git add -A +git commit -m "初始化 v1.0.0:增减挂钩后备资源调查数据库生成套件" +git push origin main +echo "" +echo "=== 完成 ===" +echo "仓库地址: $GITEA/$ORG/$REPO" +echo "推送状态: $?" diff --git a/suite.tgz b/suite.tgz new file mode 100644 index 0000000..5a9790c Binary files /dev/null and b/suite.tgz differ diff --git a/workflow.yaml b/workflow.yaml new file mode 100644 index 0000000..91105b2 --- /dev/null +++ b/workflow.yaml @@ -0,0 +1,105 @@ +name: 永善县增减挂钩后备资源调查数据库生成 +description: 从三调及年度变更数据中提取建设用地,逐层擦除限制要素,生成增减挂钩后备资源调查数据库 +version: 1.0.0 +author: SuiteForge +tags: [增减挂钩, 建设用地, 叠加分析, 后备资源, 永善县] +category: 城乡建设用地增减挂钩 + +params: + db_2019_path: + type: string + required: true + desc: 2019年三调数据库路径,支持 .gpkg / .shp / .gdb 格式 + db_2025_path: + type: string + required: true + desc: 2025年度国土变更调查数据库路径,支持 .gpkg / .shp / .gdb 格式 + restriction_dir: + type: string + required: true + desc: 限制要素目录路径,目录下存放所有限制要素矢量文件(.shp / .gpkg),包括:永久基本农田、生态保护红线、村庄建设边界、国土空间规划重大项目、城镇开发边界、矿权、生态修复项目、历史遗留矿山图斑、长江经济带项目、增减挂钩已验收项目 + output_path: + type: string + required: true + desc: 输出 GeoPackage 文件路径,如 /home/user/output/ys_houbei.gpkg + layer_2019: + type: string + required: false + desc: 2019年数据库图层名(GeoPackage 格式时指定),不填自动检测 + layer_2025: + type: string + required: false + desc: 2025年数据库图层名(GeoPackage 格式时指定),不填自动检测 + min_area: + type: number + required: false + default: 70 + desc: 最小图斑面积(平方米),小于此面积的将被自动删除 + orthophoto_dir: + type: string + required: false + desc: 正射影像目录(可选),用于生成影像核实清单时记录影像文件路径 + +base_image: gis-base:latest + +steps: + - id: extract_2019 + name: 提取2019年三调建设用地 + type: python + script_id: process + params: + mode: extract + db_path: $params.db_2019_path + layer: $params.layer_2019 + label: 2019年三调建设用地 + output: /tmp/output/step1_2019.gpkg + log: /tmp/output/process.log + + - id: extract_2025 + name: 提取2025年度变更建设用地 + type: python + script_id: process + params: + mode: extract + db_path: $params.db_2025_path + layer: $params.layer_2025 + label: 2025年度变更建设用地 + output: /tmp/output/step2_2025.gpkg + log: /tmp/output/process.log + + - id: overlay_intersect + name: 叠加分析取交集 + type: python + script_id: process + params: + mode: overlay + input1: /tmp/output/step1_2019.gpkg + input2: /tmp/output/step2_2025.gpkg + output: /tmp/output/step3_common.gpkg + log: /tmp/output/process.log + depends_on: [extract_2019, extract_2025] + + - id: erase_restrictions + name: 逐层擦除限制要素 + type: python + script_id: process + params: + mode: erase + input_path: /tmp/output/step3_common.gpkg + restriction_dir: $params.restriction_dir + output: /tmp/output/step4_erased.gpkg + log: /tmp/output/process.log + depends_on: [overlay_intersect] + + - id: filter_area + name: 面积过滤与输出成果 + type: python + script_id: process + params: + mode: filter + input_path: /tmp/output/step4_erased.gpkg + output: $params.output_path + log: /tmp/output/process.log + min_area: $params.min_area + orthophoto_dir: $params.orthophoto_dir + depends_on: [erase_restrictions]