diff --git a/.trae/rules/geo_rules.md b/.trae/rules/geo_rules.md new file mode 100644 index 0000000..0a24594 --- /dev/null +++ b/.trae/rules/geo_rules.md @@ -0,0 +1,10 @@ +# Geo-Tools 项目特定编码规范 + +## 地理数据处理黄金法则 +1. **投影先行**：任何涉及距离、面积的计算前，必须先检查 CRS，若是地理坐标系(如 EPSG:4326)，必须提示用户或自动调用 `suggest_projected_crs` 转换。 +2. **几何清洗**：任何接收用户输入的 GeoDataFrame 的公开方法，内部必须包含 `gdf = gdf[gdf.is_valid]` 或 `make_valid()` 的处理逻辑，防止脏数据导致底层 C 库崩溃。 +3. **内存控制**：读取大文件函数需要提供 `mask` (裁剪范围) 或 `rows` (读取行数限制) 参数。 + +## 代码风格 +- 空间操作函数返回新的 GeoDataFrame，尽量不修改原始输入数据（Immutable）。 +- 使用 `pathlib.Path` 处理所有文件路径，而不是字符串拼接。 \ No newline at end of file diff --git a/.trae/rules/project_rules.md b/.trae/rules/project_rules.md new file mode 100644 index 0000000..07639a8 --- /dev/null +++ b/.trae/rules/project_rules.md @@ -0,0 +1,8 @@ +# 给 AI 助手看的规则 - 小白专用版 + +## 关于内存（最重要） +1. 凡是读文件的函数，只要没加 `rows` 限制或 `chunk_size` 分块，都要在注释里用【警告】两个字提醒用户注意内存风险。 +2. 不要用 `for` 循环一行行修改 GeoDataFrame，要用 `apply` 或者向量化操作。 + +## 关于报错 +3. 禁止直接把底层库（fiona, shapely）的英文报错抛给用户看。必须包一层 try-except，翻译成中文提示，例如：“文件打不开，检查一下路径是不是有中文？”或者“数据里有烂掉的图形，已自动修复/跳过”。 \ No newline at end of file diff --git a/README.md b/README.md deleted file mode 100644 index 542a57e..0000000 --- a/README.md +++ /dev/null @@ -1,128 +0,0 @@ -# geo_tools - -> 专业地理信息数据处理工具库 —— 基于 geopandas / shapely / fiona - -[![Python](https://img.shields.io/badge/python-3.10%2B-blue)](https://www.python.org) -[![License: MIT](https://img.shields.io/badge/license-MIT-green)](LICENSE) - ---- - -## 功能特性 - -- **统一 IO 接口**：一行代码读写 Shapefile、GeoJSON、GeoPackage、**File Geodatabase (GDB)**、KML、CSV 等格式 -- **核心几何运算**：基于 Shapely 2.x 的缓冲区、集合运算、有效性检查与自动修复 -- **坐标系处理**：重投影、CRS 信息查询、批量坐标转换，内置中国常用 CRS 常量 -- **空间分析**：叠置分析、最近邻、按位置选择、面积加权均值、属性统计汇总 -- **配置驱动**：通过 `.env` 或环境变量控制输出路径、日志级别、默认 CRS 等 -- **栅格预留接口**：为 rasterio 集成预留扩展点 - -## 项目结构 - -``` -geo_tools/ -├── geo_tools/ # 主包 -│ ├── config/ # Pydantic BaseSettings 全局配置 -│ ├── core/ # 核心处理（vector / geometry / projection / raster） -│ ├── io/ # 数据读写（readers / writers，含 GDB） -│ ├── analysis/ # 空间分析（spatial_ops / stats） -│ └── utils/ # 通用工具（logger / validators / config） -├── scripts/ # 独立处理脚本 -├── tests/ # pytest 测试套件 -├── data/sample/ # 示例数据（GeoJSON） -├── output/ # 处理结果输出目录 -├── logs/ # 日志文件目录 -├── docs/ # 文档 -└── pyproject.toml # 项目配置与依赖 -``` - -## 快速开始 - -### 安装依赖 - -```bash -# 推荐使用 conda 安装地理库（避免 GDAL 编译问题） -conda install -c conda-forge geopandas shapely fiona pyproj - -# 然后安装本项目（开发模式） -pip install -e ".[dev]" -``` - -### 基本使用 - -```python -import geo_tools - -# 读取矢量数据（自动识别格式） -gdf = geo_tools.read_vector("data/sample/sample_points.geojson") - -# 读写 File Geodatabase -layers = geo_tools.list_gdb_layers("path/to/data.gdb") -gdf = geo_tools.read_gdb("path/to/data.gdb", layer="my_layer") -geo_tools.write_gdb(gdf, "output/result.gdb", layer="result") - -# 坐标系转换 -gdf_proj = geo_tools.reproject(gdf, "EPSG:3857") - -# 缓冲区分析 -from geo_tools.core.geometry import buffer_geometry -buffered_geom = buffer_geometry(gdf.geometry[0], distance=1000) - -# 空间叠置 -from geo_tools.analysis.spatial_ops import overlay -result = geo_tools.overlay(layer_a, layer_b, how="intersection") - -# 面积加权均值 -from geo_tools.analysis.stats import area_weighted_mean -result = area_weighted_mean(polygon_gdf, value_col="soil_ph", group_col="region") -``` - -### 配置 - -复制 `.env.example` 为 `.env` 并按需修改： - -```bash -GEO_TOOLS_OUTPUT_DIR=D:/output -GEO_TOOLS_DEFAULT_CRS=EPSG:4490 -GEO_TOOLS_LOG_LEVEL=DEBUG -``` - -## 运行测试 - -```bash -# 运行全部测试 -pytest tests/ -v - -# 运行带覆盖率报告 -pytest tests/ -v --cov=geo_tools --cov-report=html -``` - -## 运行示例脚本 - -```bash -python scripts/example_workflow.py -``` - -## GDB 支持说明 - -本项目通过 `fiona>=1.9` 的 `OpenFileGDB` 驱动读写 Esri File Geodatabase（`.gdb`）。 - -| 操作 | 驱动 | 要求 | -|------|------|------| -| 读取 GDB | `OpenFileGDB` | fiona >= 1.9（内置） | -| 写出 GDB | `OpenFileGDB` | fiona >= 1.9（内置） | -| 编辑 GDB（高级） | `FileGDB` | 需要 ESRI FileGDB API | - -```python -# 列出所有图层 -layers = geo_tools.list_gdb_layers("data.gdb") - -# 读取指定图层 -gdf = geo_tools.read_gdb("data.gdb", layer="土地利用", crs="EPSG:4490") - -# 写出到 GDB（新建或追加图层） -geo_tools.write_gdb(result_gdf, "output.gdb", layer="分析结果", mode="w") -``` - -## 许可证 - -MIT License diff --git a/TUTORIAL.md b/TUTORIAL.md new file mode 100644 index 0000000..3dce080 --- /dev/null +++ b/TUTORIAL.md @@ -0,0 +1,468 @@ +# Geo-Tools 使用教程 + +## 1. 安装与环境准备 + +### 1.1 前提条件 + +在安装 Geo-Tools 之前，您需要确保已经安装了以下依赖： + +- Python 3.8 或更高版本 +- pip（Python 包管理工具） + +### 1.2 安装方法 + +1. **克隆项目仓库** + + ```bash + git clone <项目地址> + cd geo_tools + ``` + +2. **安装依赖** + + ```bash + pip install -r requirements.txt + ``` + +3. **安装 Geo-Tools** + + ```bash + pip install -e . + ``` + +### 1.3 验证安装 + +运行以下命令验证 Geo-Tools 是否安装成功： + +```python +from app.io.readers import read_vector +print("Geo-Tools 安装成功！") +``` + +## 2. 5分钟体验 + +### 2.1 读取 Shapefile + +```python +from app.io.readers import read_vector +from pathlib import Path + +# 读取 Shapefile 文件 +# 这里使用示例数据，您可以替换为自己的文件路径 +data_path = Path("data/sample/sample_points.geojson") +gdf = read_vector(data_path) + +# 查看数据基本信息 +print(f"数据形状：{gdf.shape}") +print(f"数据列名：{list(gdf.columns)}") +print(f"坐标系：{gdf.crs}") +``` + +### 2.2 预览数据 + +```python +# 预览前5行数据 +print("\n数据预览：") +print(gdf.head()) + +# 或者使用 rows 参数直接只读取前5行 +print("\n使用 rows 参数预览：") +gdf_preview = read_vector(data_path, rows=5) +print(gdf_preview) +``` + +### 2.3 制作缓冲区 + +```python +from app.core.geometry import buffer_geometry + +# 为每个点创建缓冲区 +gdf["buffer"] = gdf.geometry.apply(lambda geom: buffer_geometry(geom, distance=0.1)) + +# 查看缓冲区结果 +print("\n缓冲区创建完成！") +print(f"原始几何类型：{gdf.geometry.geom_type.unique()}") +print(f"缓冲区几何类型：{gdf['buffer'].geom_type.unique()}") +``` + +### 2.4 保存结果 + +```python +from app.io.writers import write_vector + +# 保存为 GeoJSON 文件 +output_path = Path("output/buffered_points.geojson") +write_vector(gdf, output_path) +print(f"\n结果已保存至：{output_path}") +``` + +## 3. 进阶功能详解 + +### 3.1 大文件处理 + +#### 3.1.1 使用 rows 参数预览数据 + +当处理大文件时，您可以使用 `rows` 参数只读取前几行数据，快速了解数据结构，而不需要加载整个文件： + +```python +from app.io.readers import read_vector +from pathlib import Path + +# 只读取前10行数据进行预览 +large_file_path = Path("path/to/large_file.shp") +gdf_preview = read_vector(large_file_path, rows=10) + +print(f"预览数据包含 {len(gdf_preview)} 条记录") +print(gdf_preview.head()) +print(gdf_preview.columns) +``` + +#### 3.1.2 使用 chunk_size 分块读取 + +对于非常大的文件，您可以使用 `chunk_size` 参数进行分块读取，逐块处理数据，避免内存溢出： + +```python +from app.io.readers import read_vector +from pathlib import Path + +# 分块读取，每块10000条数据 +large_file_path = Path("path/to/large_file.shp") + +# 使用 for 循环处理每个数据块 +for i, chunk in enumerate(read_vector(large_file_path, chunk_size=10000)): + print(f"处理第 {i+1} 块数据，包含 {len(chunk)} 条记录") + + # 在这里进行您的处理操作 + # 例如：计算缓冲区 + chunk["buffer"] = chunk.geometry.apply(lambda geom: buffer_geometry(geom, distance=0.1)) + + # 保存当前块的结果 + output_path = Path(f"output/chunk_{i+1}.geojson") + write_vector(chunk, output_path) + print(f"第 {i+1} 块结果已保存") +``` + +### 3.2 坐标系转换 + +#### 3.2.1 什么时候需要转换坐标系？ + +- 当您需要进行距离、面积计算时，地理坐标系（如 EPSG:4326）的单位是度，不适合直接计算 +- 当您需要将数据与其他不同坐标系的数据集叠加时 +- 当您需要使用特定坐标系的工具或服务时 + +#### 3.2.2 如何转换坐标系 + +```python +from app.io.readers import read_vector +from app.core.projection import reproject +from pathlib import Path + +# 读取数据 +data_path = Path("data/sample/sample_points.geojson") +gdf = read_vector(data_path) +print(f"原始坐标系：{gdf.crs}") + +# 转换到 Web Mercator 坐标系（EPSG:3857） +gdf_3857 = reproject(gdf, "EPSG:3857") +print(f"转换后坐标系：{gdf_3857.crs}") + +# 或者在读取时直接指定目标坐标系 +gdf_direct = read_vector(data_path, crs="EPSG:3857") +print(f"直接指定坐标系：{gdf_direct.crs}") +``` + +### 3.3 空间分析 + +#### 3.3.1 缓冲区分析 + +```python +from app.io.readers import read_vector +from app.core.geometry import buffer_geometry +from pathlib import Path + +# 读取数据 +data_path = Path("data/sample/sample_points.geojson") +gdf = read_vector(data_path) + +# 创建不同距离的缓冲区 +gdf["buffer_05"] = gdf.geometry.apply(lambda geom: buffer_geometry(geom, distance=0.05)) +gdf["buffer_10"] = gdf.geometry.apply(lambda geom: buffer_geometry(geom, distance=0.1)) + +# 保存结果 +output_path = Path("output/buffers.geojson") +write_vector(gdf, output_path) +print(f"缓冲区分析结果已保存至：{output_path}") +``` + +#### 3.3.2 空间叠加分析 + +```python +from app.io.readers import read_vector +from app.analysis.spatial_ops import overlay +from pathlib import Path + +# 读取两个数据集 +points_path = Path("data/sample/sample_points.geojson") +regions_path = Path("data/sample/sample_regions.geojson") + +points = read_vector(points_path) +regions = read_vector(regions_path) + +# 执行空间叠加（交集） +overlay_result = overlay(points, regions, how="intersection") +print(f"叠加结果包含 {len(overlay_result)} 条记录") + +# 保存结果 +output_path = Path("output/overlay_result.geojson") +write_vector(overlay_result, output_path) +print(f"叠加分析结果已保存至：{output_path}") +``` + +#### 3.3.3 最近邻查找 + +```python +from app.core.geometry import distance_between +from app.io.readers import read_vector +from pathlib import Path + +# 读取点数据集 +points_path = Path("data/sample/sample_points.geojson") +gdf = read_vector(points_path) + +# 定义目标点 +target_point = gdf.geometry.iloc[0] +print(f"目标点：{target_point}") + +# 计算每个点到目标点的距离 +gdf["distance"] = gdf.geometry.apply(lambda geom: distance_between(geom, target_point)) + +# 找到最近的点 +nearest_point = gdf.loc[gdf["distance"].idxmin()] +print(f"最近的点：{nearest_point.geometry}") +print(f"距离：{nearest_point['distance']}") +``` + +## 4. 常见问题排查 + +### 4.1 报错“文件打不开”怎么办？ + +**可能原因：** +- 文件路径不存在 +- 文件路径包含中文或特殊字符 +- 文件格式不支持 +- 文件损坏 + +**解决方案：** +1. 检查文件路径是否正确，使用绝对路径 +2. 确保文件路径不包含中文或特殊字符 +3. 确认文件格式是否被支持（GeoJSON、Shapefile、GPKG、GDB等） +4. 尝试使用其他软件打开文件，确认文件是否损坏 + +**示例代码：** + +```python +from app.io.readers import read_vector +from pathlib import Path + +# 使用绝对路径 +try: + file_path = Path("c:/data/my_shapefile.shp") + gdf = read_vector(file_path) + print("文件读取成功！") +except Exception as e: + print(f"文件读取失败：{e}") + # 检查路径是否存在 + if not file_path.exists(): + print("错误：文件路径不存在") + # 检查文件扩展名 + if file_path.suffix not in [".shp", ".geojson", ".gpkg"]: + print("错误：文件格式可能不支持") +``` + +### 4.2 报错“几何无效”怎么办？ + +**可能原因：** +- 几何数据损坏 +- 几何自相交 +- 几何为空 + +**解决方案：** +1. 使用 `fix_geometry` 函数尝试修复无效几何 +2. 使用 `is_valid_geometry` 函数检查几何有效性 +3. 过滤掉无效几何 + +**示例代码：** + +```python +from app.core.geometry import fix_geometry, is_valid_geometry +from app.io.readers import read_vector +from pathlib import Path + +# 读取数据 +data_path = Path("data/sample/sample_points.geojson") +gdf = read_vector(data_path) + +# 检查并修复几何 +print("检查几何有效性...") +valid_count = 0 +fixed_count = 0 + +for i, geom in enumerate(gdf.geometry): + if is_valid_geometry(geom): + valid_count += 1 + else: + # 尝试修复 + fixed_geom = fix_geometry(geom) + if fixed_geom is not None: + gdf.geometry.iloc[i] = fixed_geom + fixed_count += 1 + +print(f"有效几何：{valid_count}") +print(f"修复几何：{fixed_count}") +print(f"无效几何：{len(gdf) - valid_count - fixed_count}") + +# 过滤掉仍然无效的几何 +gdf_valid = gdf[gdf.geometry.apply(is_valid_geometry)] +print(f"过滤后剩余几何：{len(gdf_valid)}") +``` + +### 4.3 内存爆了怎么办？ + +**可能原因：** +- 文件太大，一次性加载到内存 +- 处理过程中创建了过多临时对象 + +**解决方案：** +1. 使用 `rows` 参数预览数据 +2. 使用 `chunk_size` 参数分块读取 +3. 处理完数据后及时释放内存 +4. 考虑使用更高效的数据结构和算法 + +**示例代码：** + +```python +from app.io.readers import read_vector +from app.io.writers import write_vector +from pathlib import Path + +# 分块处理大文件 +large_file = Path("path/to/large_file.shp") +output_file = Path("output/processed_file.geojson") + +# 分块读取并处理 +chunks = [] +for i, chunk in enumerate(read_vector(large_file, chunk_size=10000)): + print(f"处理第 {i+1} 块...") + + # 在这里进行处理 + # 例如：添加面积列 + chunk["area"] = chunk.geometry.area + + chunks.append(chunk) + +# 合并所有块 +import geopandas as gpd +result = gpd.GeoDataFrame(pd.concat(chunks, ignore_index=True)) + +# 保存结果 +write_vector(result, output_file) +print(f"处理完成，结果已保存至：{output_file}") + +# 释放内存 +import gc +gc.collect() +``` + +## 5. 完整案例 + +### 5.1 土地利用数据分析 + +**场景：** 读取土地利用数据，筛选出耕地，计算耕地总面积，然后导出为 GeoJSON 文件。 + +**步骤：** + +1. **读取数据** + + ```python + from app.io.readers import read_vector + from pathlib import Path + + # 读取土地利用数据 + landuse_path = Path("data/landuse.shp") + # 先预览数据结构 + landuse_preview = read_vector(landuse_path, rows=10) + print("数据列名：", list(landuse_preview.columns)) + print("土地利用类型：", landuse_preview["type"].unique()) + + # 读取完整数据 + landuse = read_vector(landuse_path) + print(f"总数据量：{len(landuse)}") + ``` + +2. **筛选耕地** + + ```python + # 假设耕地的类型代码是 "1" 或 "耕地" + # 根据实际数据结构调整条件 + farmland = landuse[landuse["type"] == "耕地"] + print(f"耕地数量：{len(farmland)}") + ``` + +3. **计算耕地总面积** + + ```python + from app.core.projection import reproject + + # 检查坐标系 + print(f"原始坐标系：{farmland.crs}") + + # 如果是地理坐标系，转换到投影坐标系以获得准确的面积 + if farmland.crs and farmland.crs.to_epsg() == 4326: + # 转换到 UTM 坐标系（根据数据所在区域选择合适的 EPSG 代码） + farmland_proj = reproject(farmland, "EPSG:32649") # 示例：UTM 49N + print(f"转换后坐标系：{farmland_proj.crs}") + else: + farmland_proj = farmland + + # 计算面积（单位：平方米） + farmland_proj["area"] = farmland_proj.geometry.area + total_area = farmland_proj["area"].sum() + print(f"耕地总面积：{total_area:.2f} 平方米") + print(f"耕地总面积：{total_area/10000:.2f} 公顷") + ``` + +4. **导出结果** + + ```python + from app.io.writers import write_vector + + # 导出为 GeoJSON + output_path = Path("output/farmland.geojson") + write_vector(farmland_proj, output_path) + print(f"耕地数据已导出至：{output_path}") + + # 导出面积统计 + import pandas as pd + stats = pd.DataFrame({ + "总耕地数量": [len(farmland)], + "总面积（平方米）": [total_area], + "总面积（公顷）": [total_area/10000] + }) + stats_path = Path("output/farmland_stats.csv") + stats.to_csv(stats_path, index=False, encoding="utf-8-sig") + print(f"统计数据已导出至：{stats_path}") + ``` + +## 6. 总结 + +Geo-Tools 是一个功能强大的地理数据处理库，提供了丰富的空间分析工具和便捷的文件 I/O 功能。通过本教程，您应该已经掌握了： + +- 基本的文件读取和写入操作 +- 大文件的分块处理技巧 +- 坐标系转换的方法 +- 常见的空间分析操作 +- 常见问题的排查方法 + +如果您在使用过程中遇到任何问题，可以参考本教程的常见问题排查部分，或者查看项目的详细文档。 + +祝您使用愉快！ \ No newline at end of file diff --git a/app/__init__.py b/app/__init__.py new file mode 100644 index 0000000..161fbe6 --- /dev/null +++ b/app/__init__.py @@ -0,0 +1,63 @@ +""" +geo_tools +~~~~~~~~~ +专业地理信息数据处理工具库。 + +核心依赖：geopandas、shapely、fiona、pyproj。 + +快速开始 +-------- +>>> from geo_tools.io import readers +>>> from geo_tools.core import vector +>>> gdf = readers.read_vector("data/sample/sample_points.geojson") +>>> gdf_proj = vector.reproject(gdf, "EPSG:3857") +>>> print(gdf_proj.crs) + +GDB 读写 +-------- +>>> from geo_tools.io import readers, writers +>>> layers = readers.list_gdb_layers("path/to/data.gdb") +>>> gdf = readers.read_gdb("path/to/data.gdb", layer="my_layer") +>>> writers.write_gdb(gdf, "output/result.gdb", layer="result_layer") + +要素类投影 +---------- +>>> from geo_tools.core import projection +>>> gdf_proj = projection.reproject_gdf(gdf, "EPSG:4490") +>>> gdf_utm = projection.reproject_gdf(gdf, auto_utm=True) +""" + +from importlib.metadata import PackageNotFoundError, version + +# ── 版本 ────────────────────────────────────────────────────────────────────── +try: + __version__ = version("geo-tools") +except PackageNotFoundError: + __version__ = "0.1.0-dev" + +# ── 配置 & 日志 ─────────────────────────────────────────────────────────────── +from .io import readers, writers +from .config.settings import settings +from .utils.logger import get_logger, set_global_level +from .utils.validators import ( + SUPPORTED_VECTOR_EXTENSIONS, + is_supported_vector_format, + is_valid_crs, + validate_crs, + validate_geometry, + validate_vector_path, +) + +__all__ = [ + "__version__", + "settings", + # utils + "get_logger", + "set_global_level", + "is_valid_crs", + "validate_crs", + "validate_geometry", + "is_supported_vector_format", + "validate_vector_path", + "SUPPORTED_VECTOR_EXTENSIONS", +] diff --git a/app/analysis/__init__.py b/app/analysis/__init__.py new file mode 100644 index 0000000..9ff13ca --- /dev/null +++ b/app/analysis/__init__.py @@ -0,0 +1 @@ +"""geo_tools.analysis 包 —— 空间分析层。""" diff --git a/geo_tools/analysis/spatial_ops.py b/app/analysis/spatial_ops.py similarity index 93% rename from geo_tools/analysis/spatial_ops.py rename to app/analysis/spatial_ops.py index cb57d51..b542291 100644 --- a/geo_tools/analysis/spatial_ops.py +++ b/app/analysis/spatial_ops.py @@ -11,7 +11,7 @@ from typing import Any import geopandas as gpd import pandas as pd -from geo_tools.utils.logger import get_logger +from app.utils.logger import get_logger logger = get_logger(__name__) @@ -56,7 +56,7 @@ def buffer_and_overlay( result = gpd.overlay(buffered, target, how=how, keep_geom_type=False) if projected_crs: - result = result.to_crs(original_crs) + result = result.to_crs(original_crs) # type: ignore logger.info("叠置分析完成：%d 条结果", len(result)) return result @@ -77,7 +77,7 @@ def overlay( ``"symmetric_difference"``、``"identity"``。 """ if df1.crs != df2.crs: - df2 = df2.to_crs(df1.crs) + df2 = df2.to_crs(df1.crs) # type: ignore result = gpd.overlay(df1, df2, how=how, keep_geom_type=keep_geom_type) logger.debug("overlay(%s)：%d 条结果", how, len(result)) return result @@ -108,7 +108,7 @@ def nearest_features( 连接了最近 target 属性的 source GDF（可能包含重复行，每行对应一个近邻）。 """ if source.crs != target.crs: - target = target.to_crs(source.crs) + target = target.to_crs(source.crs) # type: ignore result = gpd.sjoin_nearest( source, @@ -141,7 +141,7 @@ def select_by_location( 满足条件的 source 子集。 """ if source.crs != selector.crs: - selector = selector.to_crs(source.crs) + selector = selector.to_crs(source.crs) # type: ignore joined = gpd.sjoin(source, selector, how="inner", predicate=predicate) result = source.loc[source.index.isin(joined.index)].copy() diff --git a/geo_tools/analysis/stats.py b/app/analysis/stats.py similarity index 94% rename from geo_tools/analysis/stats.py rename to app/analysis/stats.py index b9540c4..27a9f55 100644 --- a/geo_tools/analysis/stats.py +++ b/app/analysis/stats.py @@ -10,7 +10,7 @@ import geopandas as gpd import numpy as np import pandas as pd -from geo_tools.utils.logger import get_logger +from app.utils.logger import get_logger logger = get_logger(__name__) @@ -55,7 +55,7 @@ def area_weighted_mean( def _weighted(group: pd.DataFrame) -> float: return float((group[value_col] * group["_area"]).sum() / group["_area"].sum()) - result = gdf.groupby(group_col).apply(_weighted, include_groups=False).rename("area_weighted_mean") + result = gdf.groupby(group_col).apply(_weighted, include_groups=False).rename("area_weighted_mean") # type: ignore[no-untyped-call] area_sum = gdf.groupby(group_col)["_area"].sum().rename("total_area") return pd.concat([result, area_sum], axis=1).reset_index() @@ -97,7 +97,7 @@ def summarize_attributes( subset = df[columns] if group_col is None: - return subset.agg(agg_funcs).T.rename_axis("column").reset_index() + return subset.agg(agg_funcs).T.rename_axis("column").reset_index() # type: ignore[no-untyped-call] df_with_group = df[[group_col] + columns] return df_with_group.groupby(group_col)[columns].agg(agg_funcs).reset_index() @@ -125,7 +125,7 @@ def count_by_polygon( 含 ``count_col`` 列的 polygons 副本。 """ if points.crs != polygons.crs: - points = points.to_crs(polygons.crs) + points = points.to_crs(polygons.crs) # type: ignore joined = gpd.sjoin(points, polygons, how="inner", predicate="within") point_counts = joined.groupby("index_right").size().rename(count_col) diff --git a/geo_tools/config/__init__.py b/app/config/__init__.py similarity index 60% rename from geo_tools/config/__init__.py rename to app/config/__init__.py index cf3de1c..45ee430 100644 --- a/geo_tools/config/__init__.py +++ b/app/config/__init__.py @@ -1,5 +1,5 @@ """geo_tools.config 包 —— 全局配置层。""" -from geo_tools.config.settings import GeoToolsSettings, settings +from app.config.settings import GeoToolsSettings, settings __all__ = ["GeoToolsSettings", "settings"] diff --git a/geo_tools/config/project_enum.py b/app/config/project_enum.py similarity index 100% rename from geo_tools/config/project_enum.py rename to app/config/project_enum.py diff --git a/geo_tools/config/settings.py b/app/config/settings.py similarity index 100% rename from geo_tools/config/settings.py rename to app/config/settings.py diff --git a/app/core/__init__.py b/app/core/__init__.py new file mode 100644 index 0000000..42eaf5f --- /dev/null +++ b/app/core/__init__.py @@ -0,0 +1 @@ +"""geo_tools.core 包 —— 核心地理处理层。""" diff --git a/app/core/geometry.py b/app/core/geometry.py new file mode 100644 index 0000000..54c7cf7 --- /dev/null +++ b/app/core/geometry.py @@ -0,0 +1,469 @@ +""" +geo_tools.core.geometry +~~~~~~~~~~~~~~~~~~~~~~~ +基于 Shapely 2.x 的几何运算工具函数。 +""" + +from __future__ import annotations + +from typing import Literal, Sequence + +import shapely +from shapely.geometry import ( + LinearRing, + LineString, + MultiLineString, + MultiPoint, + MultiPolygon, + Point, + Polygon, +) +from shapely.geometry.base import BaseGeometry + +from app.utils.logger import get_logger +from app.utils.validators import ensure_valid_geometry + +logger = get_logger(__name__) + + +# ── 几何有效性 ──────────────────────────────────────────────────────────────── + +def is_valid_geometry(geom: BaseGeometry | None) -> bool: + """判断几何对象是否有效（非空且通过 Shapely 合法性检查）。 + + Parameters + ---------- + geom: + 输入几何对象，可为 None。 + + Returns + ------- + bool + 如果几何对象有效且非空，返回 True；否则返回 False。 + """ + if geom is None: + return False + return bool(geom.is_valid and not geom.is_empty) + + +def fix_geometry(geom: BaseGeometry | None) -> BaseGeometry | None: + """尝试修复无效几何。 + + 依次尝试： + 1. ``buffer(0)`` — 适合大多数自相交多边形 + 2. ``make_valid``（Shapely 2.x）— 覆盖更多情形 + + Parameters + ---------- + geom: + 输入几何对象，可为 None。 + + Returns + ------- + BaseGeometry | None + 修复后的几何；无法修复时返回 ``None``。 + + Notes + ----- + 对于复杂的无效几何，可能无法完全修复，此时会返回 None。 + """ + if geom is None: + return None + if geom.is_valid: + return geom + + # 方法一：buffer(0) + try: + fixed = geom.buffer(0) + if fixed.is_valid and not fixed.is_empty: + return fixed + except Exception: + pass + + # 方法二：shapely.make_valid（Shapely >= 1.8） + try: + fixed = shapely.make_valid(geom) + if fixed.is_valid and not fixed.is_empty: + return fixed + except Exception: + pass + + logger.warning("无法修复几何：%r", geom.geom_type) + return None + + +def explain_validity(geom: BaseGeometry) -> str: + """返回 Shapely 对该几何的有效性说明（英文）。 + + Parameters + ---------- + geom: + 输入几何对象。 + + Returns + ------- + str + Shapely 生成的有效性说明字符串。 + """ + from shapely.validation import explain_validity as _explain + return _explain(geom) + + +# ── 基础几何运算 ─────────────────────────────────────────────────────────────── + +def buffer_geometry( + geom: BaseGeometry, + distance: float, + cap_style: Literal["round", "square", "flat"] = "round", + join_style: Literal["round", "mitre", "bevel"] = "round", + resolution: int = 16, + verbose: bool = False, +) -> BaseGeometry | None: + """对几何对象执行缓冲区运算。 + + Parameters + ---------- + geom: + 输入几何。 + distance: + 缓冲距离（单位与 CRS 一致；地理坐标系单位为度）。 + cap_style: + 端头样式："round"（圆角）、"square"（方角）、"flat"（平角）（仅线要素有效）。 + join_style: + 转角样式："round"（圆角）、"mitre"（斜角）、"bevel"（尖角）。 + resolution: + 圆弧逼近精度（段数），默认 16。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + 缓冲区运算后的几何对象；失败时返回 None。 + """ + try: + geom = ensure_valid_geometry(geom, verbose) + return geom.buffer(distance, cap_style=cap_style, join_style=join_style, resolution=resolution) + except Exception as e: + logger.warning(f"缓冲区计算失败，已跳过，原因：{str(e)}") + return None + + +def centroid(geom: BaseGeometry, verbose: bool = False) -> Point | None: + """返回几何的质心点。 + + Parameters + ---------- + geom: + 输入几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + Point | None + 几何的质心点；失败时返回 None。 + """ + try: + geom = ensure_valid_geometry(geom, verbose) + return geom.centroid + except Exception as e: + logger.warning(f"质心计算失败，已跳过，原因：{str(e)}") + return None + + +def bounding_box(geom: BaseGeometry, verbose: bool = False) -> Polygon | None: + """返回几何的最小外接矩形（BBOX）为多边形。 + + Parameters + ---------- + geom: + 输入几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + Polygon | None + 最小外接矩形多边形；失败时返回 None。 + """ + try: + geom = ensure_valid_geometry(geom, verbose) + from shapely.geometry import box + return box(*geom.bounds) + except Exception as e: + logger.warning(f"最小外接矩形计算失败，已跳过，原因：{str(e)}") + return None + + +def convex_hull(geom: BaseGeometry, verbose: bool = False) -> BaseGeometry | None: + """返回几何的凸包。 + + Parameters + ---------- + geom: + 输入几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + 几何的凸包；失败时返回 None。 + """ + try: + geom = ensure_valid_geometry(geom, verbose) + return geom.convex_hull + except Exception as e: + logger.warning(f"凸包计算失败，已跳过，原因：{str(e)}") + return None + + +# ── 集合运算 ────────────────────────────────────────────────────────────────── + +def intersect(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> BaseGeometry | None: + """返回两几何的交集。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + 两几何的交集；失败时返回 None。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return geom_a.intersection(geom_b) + except Exception as e: + logger.warning(f"交集计算失败，已跳过，原因：{str(e)}") + return None + + +def union(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> BaseGeometry | None: + """返回两几何的并集。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + 两几何的并集；失败时返回 None。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return geom_a.union(geom_b) + except Exception as e: + logger.warning(f"并集计算失败，已跳过，原因：{str(e)}") + return None + + +def difference(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> BaseGeometry | None: + """返回 ``geom_a`` 减去 ``geom_b`` 的差集。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + geom_a 减去 geom_b 的差集；失败时返回 None。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return geom_a.difference(geom_b) + except Exception as e: + logger.warning(f"差集计算失败，已跳过，原因：{str(e)}") + return None + + +def symmetric_difference(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> BaseGeometry | None: + """返回两几何的对称差集（异或）。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + 两几何的对称差集；失败时返回 None。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return geom_a.symmetric_difference(geom_b) + except Exception as e: + logger.warning(f"对称差集计算失败，已跳过，原因：{str(e)}") + return None + + +def unary_union(geoms: Sequence[BaseGeometry], verbose: bool = False) -> BaseGeometry | None: + """将多个几何合并为一个（等同于逐一 union）。 + + Parameters + ---------- + geoms: + 几何对象序列。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + BaseGeometry | None + 合并后的几何对象；失败时返回 None。 + """ + try: + # 检查几何有效性并尝试修复 + fixed_geoms = [] + for i, geom in enumerate(geoms): + try: + fixed = ensure_valid_geometry(geom, verbose) + if fixed is not None: + fixed_geoms.append(fixed) + else: + # 无法修复的几何跳过 + if verbose: + logger.warning(f"几何对象 {i} 无效且无法修复，已跳过") + except Exception as e: + logger.warning(f"几何对象 {i} 处理失败，已跳过，原因：{str(e)}") + if not fixed_geoms: + logger.warning("没有有效的几何对象可合并") + return None + return shapely.unary_union(fixed_geoms) + except Exception as e: + logger.warning(f"几何合并失败，已跳过，原因：{str(e)}") + return None + + +# ── 空间关系判断 ─────────────────────────────────────────────────────────────── + +def contains(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> bool: + """判断 ``geom_a`` 是否完全包含 ``geom_b``。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + bool + 如果 geom_a 完全包含 geom_b，返回 True；否则返回 False；失败时返回 False。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return bool(geom_a.contains(geom_b)) + except Exception as e: + logger.warning(f"包含关系判断失败，已返回 False，原因：{str(e)}") + return False + + +def within(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> bool: + """判断 ``geom_a`` 是否完全在 ``geom_b`` 内。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + bool + 如果 geom_a 完全在 geom_b 内，返回 True；否则返回 False；失败时返回 False。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return bool(geom_a.within(geom_b)) + except Exception as e: + logger.warning(f"包含于关系判断失败，已返回 False，原因：{str(e)}") + return False + + +def intersects(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> bool: + """判断两几何是否相交（含边界接触）。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + bool + 如果两几何相交，返回 True；否则返回 False；失败时返回 False。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return bool(geom_a.intersects(geom_b)) + except Exception as e: + logger.warning(f"相交关系判断失败，已返回 False，原因：{str(e)}") + return False + + +def distance_between(geom_a: BaseGeometry, geom_b: BaseGeometry, verbose: bool = False) -> float: + """计算两几何间的最小距离（单位与 CRS 一致）。 + + Parameters + ---------- + geom_a: + 第一个几何。 + geom_b: + 第二个几何。 + verbose: + 是否输出修复几何的警告信息，默认 False。 + + Returns + ------- + float + 两几何间的最小距离；失败时返回无穷大。 + """ + try: + geom_a = ensure_valid_geometry(geom_a, verbose) + geom_b = ensure_valid_geometry(geom_b, verbose) + return geom_a.distance(geom_b) + except Exception as e: + logger.warning(f"距离计算失败，已返回无穷大，原因：{str(e)}") + return float('inf') diff --git a/geo_tools/core/projection.py b/app/core/projection.py similarity index 91% rename from geo_tools/core/projection.py rename to app/core/projection.py index cbd424e..58c12c3 100644 --- a/geo_tools/core/projection.py +++ b/app/core/projection.py @@ -9,17 +9,22 @@ from typing import Sequence from pyproj import CRS, Transformer import geopandas as gpd -from geo_tools.utils.logger import get_logger +from app.utils.logger import get_logger logger = get_logger(__name__) def get_crs_info(crs_input: str | int) -> dict[str, str | int | None]: """获取 CRS 的基本信息。 + Parameters + ---------- + crs_input: + EPSG 代码（整数或 ``"EPSG:4326"`` 字符串）或 proj 字符串。 + Returns ------- dict - 包含 ``name``、``epsg``、``unit``、``is_geographic``、``is_projected``。 + 包含 ``name``、``epsg``、``unit``、``is_geographic``、``is_projected``、``datum``。 """ crs = CRS.from_user_input(crs_input) return { @@ -152,7 +157,7 @@ def reproject_gdf( 与 ``auto_crs=True`` 二选一。 auto_crs: 为 ``True`` 时忽略 ``target_crs``，根据数据中心点自动推荐 - CGCS2000 高斯-克吕格 带号（适合面积/距离计算场景）。 + CGCS2000 高斯-克吕格 带号（默认使用 3度分带）。 verbose: 为 ``True`` 时在日志中打印投影前后的 CRS 信息。 @@ -171,11 +176,11 @@ def reproject_gdf( >>> # 指定目标 CRS >>> gdf_proj = reproject_gdf(gdf, "EPSG:4490") - >>> # 自动推荐 CGCS2000 高斯-克吕格 带号（用于面积计算） + >>> # 自动推荐 CGCS2000 高斯-克吕格 带号 默认使用 3度分带 >>> gdf_utm = reproject_gdf(gdf, auto_crs=True) >>> # 配合 GDB 读取 - >>> gdf = read_gdb("data.gdb", layer="图斑") + >>> gdf = geo_tools.readers.read_vector("data.gdb/图斑") >>> gdf_proj = reproject_gdf(gdf, "EPSG:4326") """ @@ -186,9 +191,9 @@ def reproject_gdf( if auto_crs: # 先统一到地理坐标系，再取中心点推荐 CGCS2000 高斯-克吕格 带号 if gdf.crs.is_projected: - center = gdf.to_crs("EPSG:4490").geometry.unary_union.centroid + center = gdf.to_crs("EPSG:4490").geometry.union_all().centroid else: - center = gdf.geometry.unary_union.centroid + center = gdf.geometry.union_all().centroid target_crs = suggest_projected_crs(center.x, center.y) logger.info("auto_crs：自动推荐投影 CRS = %s", target_crs) diff --git a/geo_tools/core/raster.py b/app/core/raster.py similarity index 100% rename from geo_tools/core/raster.py rename to app/core/raster.py diff --git a/geo_tools/core/vector.py b/app/core/vector.py similarity index 92% rename from geo_tools/core/vector.py rename to app/core/vector.py index f215c27..d9c14d8 100644 --- a/geo_tools/core/vector.py +++ b/app/core/vector.py @@ -6,13 +6,13 @@ geo_tools.core.vector from __future__ import annotations -from typing import Any +from typing import Any, Literal import geopandas as gpd import pandas as pd -from geo_tools.utils.logger import get_logger -from geo_tools.utils.validators import validate_geometry +from app.utils.logger import get_logger +from app.utils.validators import validate_geometry logger = get_logger(__name__) @@ -82,7 +82,7 @@ def clip_to_extent( result = gdf.clip(mask) else: if bbox.crs != gdf.crs: - bbox = bbox.to_crs(gdf.crs) + bbox = bbox.to_crs(gdf.crs) # type: ignore result = gdf.clip(bbox) logger.debug("裁切完成：%d → %d 条", len(gdf), len(result)) @@ -119,6 +119,11 @@ def dissolve_by( def explode_multipart(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame: """将多部分几何（MultiPolygon 等）拆分为单部分要素。 + Parameters + ---------- + gdf: + 输入 GeoDataFrame。 + Returns ------- gpd.GeoDataFrame @@ -144,10 +149,10 @@ def drop_invalid_geometries(gdf: gpd.GeoDataFrame, *, fix: bool = False) -> gpd. return gdf if fix: - from geo_tools.core.geometry import fix_geometry + from app.core.geometry import fix_geometry gdf = gdf.copy() mask = ~gdf.geometry.is_valid | gdf.geometry.isna() - gdf.loc[mask, "geometry"] = gdf.loc[mask, "geometry"].apply(fix_geometry) + gdf.loc[mask, "geometry"] = gdf.loc[mask, "geometry"].apply(fix_geometry) # type: ignore logger.info("已修复 %d 个无效几何", stats["invalid"]) else: before = len(gdf) @@ -159,7 +164,7 @@ def drop_invalid_geometries(gdf: gpd.GeoDataFrame, *, fix: bool = False) -> gpd. def spatial_join( left: gpd.GeoDataFrame, right: gpd.GeoDataFrame, - how: str = "left", + how: Literal["left", "right", "inner"] = "left", predicate: str = "intersects", **kwargs: Any, ) -> gpd.GeoDataFrame: @@ -177,7 +182,7 @@ def spatial_join( 空间谓词：``"intersects"``、``"contains"``、``"within"``、``"touches"``。 """ if left.crs != right.crs: - right = right.to_crs(left.crs) + right = right.to_crs(left.crs) # type: ignore result = gpd.sjoin(left, right, how=how, predicate=predicate, **kwargs) logger.debug("空间连接完成：%d 条结果", len(result)) return result diff --git a/app/io/__init__.py b/app/io/__init__.py new file mode 100644 index 0000000..c738e2c --- /dev/null +++ b/app/io/__init__.py @@ -0,0 +1 @@ +"""geo_tools.io 包 —— 数据读写层。""" diff --git a/geo_tools/io/readers.py b/app/io/readers.py similarity index 57% rename from geo_tools/io/readers.py rename to app/io/readers.py index 4f3e0b3..9201804 100644 --- a/geo_tools/io/readers.py +++ b/app/io/readers.py @@ -15,14 +15,15 @@ geo_tools.io.readers from __future__ import annotations +import os from pathlib import Path -from typing import Any +from typing import Any, Generator import fiona import geopandas as gpd -from geo_tools.utils.logger import get_logger -from geo_tools.utils.validators import validate_vector_path +from app.utils.logger import get_logger +from app.utils.validators import validate_vector_path logger = get_logger(__name__) @@ -34,8 +35,10 @@ def read_vector( layer: str | int | None = None, crs: str | int | None = None, encoding: str = "utf-8", + chunk_size: int | None = None, + rows: int | None = None, **kwargs: Any, -) -> gpd.GeoDataFrame: +): """统一的矢量数据读取入口，自动识别文件格式。 Parameters @@ -48,19 +51,46 @@ def read_vector( 读取后强制重投影到目标 CRS（不传则保留原始 CRS）。 encoding: 属性表编码，Shapefile 中文路径常需指定 ``"gbk"``。 + chunk_size: + 分块大小，默认 None（一次性读取全部数据）。 + 【警告】：若不设置 chunk_size，大文件可能会占用大量内存。 + rows: + 限制读取的行数，默认 None（读取全部数据）。 + 用于快速预览数据，避免读取大文件的全部内容。 **kwargs: 透传给 ``geopandas.read_file`` 的额外参数。 Returns ------- - gpd.GeoDataFrame + gpd.GeoDataFrame 或生成器 + 如果 chunk_size 为 None，返回完整的 GeoDataFrame； + 如果设置了 chunk_size，返回一个生成器，每次 yield 一个 GeoDataFrame 块。 + + 示例 + ----- + # 全量读取（老方法） + gdf = read_vector("data.shp") + + # 分块读取（新方法） + for chunk in read_vector("large_data.shp", chunk_size=10000): + # 处理每个数据块 + print(f"处理了 {len(chunk)} 条数据") + # 在这里做你的操作，比如计算、过滤等 + + # 只读取前 5 行数据（预览模式） + gdf_preview = read_vector("large_data.shp", rows=5) + print(gdf_preview.head()) """ + path, layer = _split_gdb_layer(path) path = validate_vector_path(path) suffix = path.suffix.lower() logger.info("读取矢量数据：%s（格式：%s，图层：%s）", path, suffix or "目录", layer) if suffix == ".csv": + # CSV 文件暂时不支持分块读取 + if chunk_size is not None: + logger.warning("CSV 文件暂不支持分块读取，将一次性读取全部数据") return _read_csv_vector(path, crs=crs, **kwargs) # fiona / geopandas 通用读取 @@ -68,14 +98,60 @@ def read_vector( if layer is not None: read_kwargs["layer"] = layer - gdf = gpd.read_file(str(path), **read_kwargs) + # 分块读取模式 + if chunk_size is not None: + def _chunk_generator(): + logger.info("启用分块读取模式，每块 %d 条数据", chunk_size) + try: + # 使用 fiona 打开文件 + with fiona.open(str(path), **read_kwargs) as src: + # 获取坐标系信息 + crs_info = src.crs + # 分块读取 + features = [] + for i, feature in enumerate(src): + # 检查是否达到行数限制 + if rows is not None and i >= rows: + break + + features.append(feature) + if (i + 1) % chunk_size == 0: + # 创建 GeoDataFrame 并设置 CRS + gdf = gpd.GeoDataFrame.from_features(features, crs=crs_info) + # 重投影 + if crs is not None: + gdf = gdf.to_crs(crs) # type: ignore + logger.debug("读取并处理第 %d 块数据，共 %d 条", (i + 1) // chunk_size, len(gdf)) + yield gdf + features = [] + # 处理最后一块 + if features: + gdf = gpd.GeoDataFrame.from_features(features, crs=crs_info) + if crs is not None: + gdf = gdf.to_crs(crs) # type: ignore + logger.debug("读取并处理最后一块数据，共 %d 条", len(gdf)) + yield gdf + except Exception as exc: + raise RuntimeError(f"无法分块读取矢量数据：{exc}") from None + return _chunk_generator() + else: + # 一次性读取模式 + try: + # 添加 rows 参数到读取参数中 + if rows is not None: + read_kwargs["rows"] = rows + logger.info("限制读取行数：%d", rows) + + gdf = gpd.read_file(str(path), **read_kwargs) + except Exception as exc: + raise RuntimeError(f"无法读取矢量数据：{exc}") from None - if crs is not None: - logger.debug("重投影到 %s", crs) - gdf = gdf.to_crs(crs) + if crs is not None: + logger.debug("重投影到 %s", crs) + gdf = gdf.to_crs(crs) # type: ignore - logger.info("读取完成：共 %d 条要素，CRS=%s", len(gdf), gdf.crs) - return gdf + logger.info("读取完成：共 %d 条要素，CRS=%s", len(gdf), gdf.crs) + return gdf # type: ignore # ── GDB 专用 ─────────────────────────────────────────────────────────────────── @@ -125,10 +201,10 @@ def read_gdb( gdf = gpd.read_file(str(gdb_path), layer=layer, encoding=encoding, **kwargs) if crs is not None: - gdf = gdf.to_crs(crs) + gdf = gdf.to_crs(crs) # type: ignore logger.info("GDB 读取完成：%d 条要素，CRS=%s", len(gdf), gdf.crs) - return gdf + return gdf # type: ignore def list_gdb_layers(gdb_path: str | Path) -> list[str]: @@ -154,6 +230,42 @@ def list_gdb_layers(gdb_path: str | Path) -> list[str]: f"原始错误：{exc}" ) from exc +def _split_gdb_layer(path: str | Path) -> tuple[Path, str | None]: + """从完整路径中分离 GDB 数据库路径和图层名。 + + Parameters + ---------- + path: + 完整路径，可以是字符串或 Path 对象。 + + Returns + ------- + tuple[Path, str | None] + (gdb_path, layer_name)，其中 gdb_path 是 GDB 目录路径，layer_name 是图层名，若没有图层名则为 None。 + """ + path_obj = Path(path) + str_path = str(path_obj) + + # 查找 .gdb 的位置 + gdb_pos = str_path.find('.gdb') + + if gdb_pos == -1: + # 如果没有 .gdb，整个路径作为 GDB 路径，没有图层 + return path_obj, None + + # 提取 GDB 路径（包含 .gdb） + gdb_path = str_path[:gdb_pos + 4] + + # 提取图层名（.gdb 之后的部分） + layer_part = str_path[gdb_pos + 4:] + # 去除开头的路径分隔符 + layer_name = layer_part.lstrip(os.sep).lstrip('/').lstrip('\\') + + # 如果没有图层名，返回 None + if not layer_name: + layer_name = None + + return Path(gdb_path), layer_name # ── GPKG 专用 ────────────────────────────────────────────────────────────────── @@ -188,12 +300,23 @@ def read_gpkg( gdf = gpd.read_file(str(gpkg_path), layer=layer, **kwargs) if crs is not None: - gdf = gdf.to_crs(crs) - return gdf + gdf = gdf.to_crs(crs) # type: ignore + return gdf # type: ignore def list_gpkg_layers(gpkg_path: str | Path) -> list[str]: - """列出 GeoPackage 中所有图层名称。""" + """列出 GeoPackage 中所有图层名称。 + + Parameters + ---------- + gpkg_path: + GeoPackage 文件路径，可以是字符串或 Path 对象。 + + Returns + ------- + list[str] + 图层名称列表。 + """ return fiona.listlayers(str(gpkg_path)) diff --git a/geo_tools/io/writers.py b/app/io/writers.py similarity index 99% rename from geo_tools/io/writers.py rename to app/io/writers.py index 63d4619..8e56e20 100644 --- a/geo_tools/io/writers.py +++ b/app/io/writers.py @@ -17,7 +17,7 @@ from typing import Any, Literal import geopandas as gpd -from geo_tools.utils.logger import get_logger +from app.utils.logger import get_logger logger = get_logger(__name__) diff --git a/geo_tools/utils/__init__.py b/app/utils/__init__.py similarity index 73% rename from geo_tools/utils/__init__.py rename to app/utils/__init__.py index b8a88e0..b96baa4 100644 --- a/geo_tools/utils/__init__.py +++ b/app/utils/__init__.py @@ -1,8 +1,8 @@ """geo_tools.utils 包 —— 通用工具函数。""" -from geo_tools.utils.config import load_config, load_json_config, load_toml_config, load_yaml_config -from geo_tools.utils.logger import get_logger, set_global_level -from geo_tools.utils.validators import ( +from app.utils.config import load_config, load_json_config, load_toml_config, load_yaml_config +from app.utils.logger import get_logger, set_global_level +from app.utils.validators import ( SUPPORTED_VECTOR_EXTENSIONS, is_supported_vector_format, is_valid_crs, diff --git a/geo_tools/utils/config.py b/app/utils/config.py similarity index 100% rename from geo_tools/utils/config.py rename to app/utils/config.py diff --git a/geo_tools/utils/logger.py b/app/utils/logger.py similarity index 98% rename from geo_tools/utils/logger.py rename to app/utils/logger.py index c1cd4ea..f2621f8 100644 --- a/geo_tools/utils/logger.py +++ b/app/utils/logger.py @@ -50,7 +50,7 @@ def get_logger( logging.Logger """ # 延迟导入，避免循环依赖 - from geo_tools.config.settings import settings as _settings + from app.config.settings import settings as _settings if level is None: level = _settings.log_level diff --git a/geo_tools/utils/validators.py b/app/utils/validators.py similarity index 79% rename from geo_tools/utils/validators.py rename to app/utils/validators.py index 32c0bad..f75a5d6 100644 --- a/geo_tools/utils/validators.py +++ b/app/utils/validators.py @@ -73,11 +73,9 @@ def validate_geometry(gdf: "gpd.GeoDataFrame", *, raise_on_invalid: bool = False dict 包含 ``total``、``valid``、``invalid``、``null`` 计数。 """ - import geopandas as gpd # noqa: F811 - null_count = gdf.geometry.isna().sum() non_null = gdf.geometry.dropna() - invalid_mask = ~non_null.is_valid + invalid_mask = ~non_null.is_valid # type: ignore invalid_count = int(invalid_mask.sum()) valid_count = len(non_null) - invalid_count @@ -114,7 +112,18 @@ SUPPORTED_VECTOR_EXTENSIONS: dict[str, str] = { def is_supported_vector_format(path: str | Path) -> bool: - """判断路径是否为已知的矢量格式。""" + """判断路径是否为已知的矢量格式。 + + Parameters + ---------- + path: + 输入路径，可以是字符串或 Path 对象。 + + Returns + ------- + bool + 如果路径是支持的矢量格式，返回 True；否则返回 False。 + """ path = Path(path) suffix = path.suffix.lower() # .gdb 可能是目录（FileGDB） @@ -143,3 +152,32 @@ def validate_vector_path(path: str | Path) -> Path: f"支持的格式：{list(SUPPORTED_VECTOR_EXTENSIONS.keys())}" ) return path + + +def ensure_valid_geometry(geom: 'BaseGeometry', verbose: bool = False) -> 'BaseGeometry': + """确保几何对象有效，无效时尝试修复。 + + Parameters + ---------- + geom: + 输入几何。 + verbose: + 是否输出修复几何的警告信息。 + + Returns + ------- + BaseGeometry + 有效的几何对象（可能是修复后的）。 + """ + from app.core.geometry import fix_geometry + from app.utils.logger import get_logger + + logger = get_logger(__name__) + + if not geom.is_valid: + fixed_geom = fix_geometry(geom) + if fixed_geom is not None: + if verbose: + logger.warning("几何对象无效，已自动修复") + return fixed_geom + return geom diff --git a/docs/projection_guide.md b/docs/projection_guide.md deleted file mode 100644 index bd11818..0000000 --- a/docs/projection_guide.md +++ /dev/null @@ -1,248 +0,0 @@ -# geo_tools.core.projection 使用说明 - -> 坐标系查询、坐标转换、投影推荐工具，基于 [pyproj](https://pyproj4.github.io/pyproj/)。 - ---- - -## 导入方式 - -```python -# 推荐：从顶层包导入 -from geo_tools.core.projection import ( - get_crs_info, - crs_to_epsg, - transform_coordinates, - transform_point, - suggest_projected_crs, - WGS84, CGCS2000, WEB_MERCATOR, CGCS2000_UTM_50N, -) - -# 或直接通过 geo_tools 导入 -import geo_tools -``` - ---- - -## CRS 常量 - -模块内置了中国地理信息处理中最常用的 CRS 快捷常量，可直接作为参数传入所有函数： - -| 常量名 | EPSG | 说明 | -|--------|------|------| -| `WGS84` | `EPSG:4326` | WGS84 地理坐标系（经纬度，最通用） | -| `CGCS2000` | `EPSG:4490` | 中国国家大地坐标系 2000（经纬度） | -| `WEB_MERCATOR` | `EPSG:3857` | Web Mercator 投影（网络地图常用，单位：米） | -| `CGCS2000_UTM_50N` | `EPSG:4508` | CGCS2000 / 3° 高斯-克吕格 50 带（单位：米） | - -```python -from geo_tools.core.projection import WGS84, CGCS2000 - -# 直接用常量替代字符串 -gdf = gdf.to_crs(CGCS2000) -``` - ---- - -## 函数说明 - -### `get_crs_info(crs_input)` — 查询 CRS 信息 - -返回坐标系的详细描述字典，方便快速了解一个未知 EPSG 的含义。 - -**参数** -- `crs_input`：EPSG 代码字符串（如 `"EPSG:4523"`）、整数（如 `4523`）或 proj 字符串。 - -**返回值**（`dict`） - -| 键 | 含义 | -|----|------| -| `name` | 坐标系名称 | -| `epsg` | EPSG 整数编号（无法识别时为 `None`） | -| `unit` | 坐标单位（`degree` / `metre`） | -| `is_geographic` | 是否为地理坐标系（经纬度） | -| `is_projected` | 是否为投影坐标系（平面直角） | -| `datum` | 基准面名称 | - -```python -from geo_tools.core.projection import get_crs_info - -# 查询读取到的 GDB 数据的 CRS 含义 -info = get_crs_info("EPSG:4523") -print(info) -# { -# 'name': 'CGCS2000 / 3-degree Gauss-Kruger zone 45', -# 'epsg': 4523, -# 'unit': 'metre', -# 'is_geographic': False, -# 'is_projected': True, -# 'datum': 'China Geodetic Coordinate System 2000' -# } - -# 直接传整数 -info = get_crs_info(32650) -print(info["name"]) # WGS 84 / UTM zone 50N -``` - ---- - -### `crs_to_epsg(crs_input)` — 获取 EPSG 编号 - -将任意 CRS 描述转为整数 EPSG 编号，无法识别时返回 `None`（不抛异常）。 - -```python -from geo_tools.core.projection import crs_to_epsg - -epsg = crs_to_epsg("EPSG:4490") -print(epsg) # 4490 - -epsg = crs_to_epsg("WGS 84") -print(epsg) # 4326 - -epsg = crs_to_epsg("invalid_crs") -print(epsg) # None -``` - ---- - -### `transform_coordinates(xs, ys, source_crs, target_crs)` — 批量坐标转换 - -将一组坐标点从源坐标系批量转换到目标坐标系，返回转换后的 `(xs, ys)` 列表。 - -**参数** -- `xs`：X 坐标序列（地理 CRS 时为**经度**） -- `ys`：Y 坐标序列（地理 CRS 时为**纬度**） -- `source_crs`：源坐标系 -- `target_crs`：目标坐标系 -- `always_xy`（关键字参数）：强制按 (经度/X, 纬度/Y) 顺序处理，默认 `True`，**建议不修改** - -```python -from geo_tools.core.projection import transform_coordinates, WGS84, WEB_MERCATOR - -# 将北京、上海、广州的 WGS84 经纬度转为 Web Mercator 米制坐标 -lons = [116.4074, 121.4737, 113.2644] -lats = [39.9042, 31.2304, 23.1291] - -xs, ys = transform_coordinates(lons, lats, WGS84, WEB_MERCATOR) -print(xs) # [12959618.8, 13521606.3, 12608870.0]（单位：米） -print(ys) # [4859767.2, 3649094.2, 2641877.0] - -# 国家坐标系转换：CGCS2000 经纬度 → CGCS2000 3° 高斯带（50带） -from geo_tools.core.projection import CGCS2000, CGCS2000_UTM_50N -xs_proj, ys_proj = transform_coordinates(lons, lats, CGCS2000, CGCS2000_UTM_50N) -``` - ---- - -### `transform_point(x, y, source_crs, target_crs)` — 单点坐标转换 - -`transform_coordinates` 的单点版本，直接返回 `(x, y)` 元组。 - -```python -from geo_tools.core.projection import transform_point, WGS84, CGCS2000 - -# 单点：WGS84 → CGCS2000（两者数值非常接近，差异在毫米级） -x, y = transform_point(116.4074, 39.9042, WGS84, CGCS2000) -print(f"CGCS2000 坐标：经度={x:.6f}, 纬度={y:.6f}") - -# 单点：经纬度 → 投影坐标（米） -from geo_tools.core.projection import WEB_MERCATOR -mx, my = transform_point(116.4074, 39.9042, WGS84, WEB_MERCATOR) -print(f"墨卡托坐标：X={mx:.2f}m, Y={my:.2f}m") -``` - ---- - -### `suggest_projected_crs(lon, lat)` — 自动推荐投影 CRS - -根据数据中心坐标（WGS84 经纬度）自动推荐适合**面积/距离计算**的 UTM 投影带，避免在地理坐标系下计算面积出错。 - -**参数** -- `lon`：中心经度（WGS84） -- `lat`：中心纬度（WGS84，北半球为正） - -**返回值**：EPSG 代码字符串，如 `"EPSG:32650"` - -```python -from geo_tools.core.projection import suggest_projected_crs - -# 云南马关县（约 104.4°E, 23.0°N） -proj_crs = suggest_projected_crs(lon=104.4, lat=23.0) -print(proj_crs) # EPSG:32648 (WGS84 UTM zone 48N) - -# 北京（116.4°E, 39.9°N） -proj_crs = suggest_projected_crs(lon=116.4, lat=39.9) -print(proj_crs) # EPSG:32650 (WGS84 UTM zone 50N) - -# 实际场景：读取 GDB 后用推荐的投影计算面积 -import geo_tools - -gdf = geo_tools.read_gdb("data.gdb", layer="图斑") -cx, cy = gdf.geometry.unary_union.centroid.x, gdf.geometry.unary_union.centroid.y - -# 如果数据是投影坐标系（单位：米），先转到地理坐标系再推荐 -if gdf.crs.is_projected: - cx, cy = geo_tools.transform_point(cx, cy, gdf.crs, "EPSG:4326") - -proj_crs = suggest_projected_crs(cx, cy) -gdf_proj = geo_tools.reproject(gdf, proj_crs) # 重投影 -gdf_proj = geo_tools.add_area_column(gdf_proj) # 计算面积（单位：m²） -``` - ---- - -## 常见场景示例 - -### 场景一：不认识数据的 CRS，先查一下 - -```python -import geo_tools - -gdf = geo_tools.read_gdb("临时数据库.gdb", layer="马关综合后图斑") -# 读取完成：CRS=EPSG:4523 - -info = geo_tools.get_crs_info(gdf.crs) -print(info["name"]) # CGCS2000 / 3-degree Gauss-Kruger zone 45 -print(info["unit"]) # metre（投影坐标系，单位是米） -print(info["is_projected"]) # True -``` - -### 场景二：统一坐标系后叠置分析 - -```python -import geo_tools -from geo_tools.core.projection import CGCS2000 - -layer_a = geo_tools.read_gdb("a.gdb", layer="林地") # EPSG:4523 -layer_b = geo_tools.read_vector("b.geojson") # EPSG:4326 - -# 统一到 CGCS2000 地理坐标系后再做叠置 -layer_a = geo_tools.reproject(layer_a, CGCS2000) -layer_b = geo_tools.reproject(layer_b, CGCS2000) - -result = geo_tools.overlay(layer_a, layer_b, how="intersection") -``` - -### 场景三：在地理坐标系数据上正确计算面积 - -```python -import geo_tools -from geo_tools.core.projection import suggest_projected_crs - -gdf = geo_tools.read_vector("data.geojson") # EPSG:4326，单位是度 - -# 自动推荐合适的投影 -proj = suggest_projected_crs(lon=105.0, lat=25.0) # 云贵地区 - -gdf = geo_tools.add_area_column(gdf, projected_crs=proj) -print(gdf["area_m2"].describe()) -``` - ---- - -## 注意事项 - -> [!WARNING] -> 在**地理坐标系**（EPSG:4326 / 4490）下直接调用 `geometry.area` 得到的是"平方度"，**不是平方米**，面积计算会严重失真。始终用 `add_area_column()` 或先 `reproject()` 到投影坐标系后再计算。 - -> [!NOTE] -> `WGS84`（EPSG:4326）与 `CGCS2000`（EPSG:4490）的坐标数值差异极小（通常 < 1 米），在普通精度的分析中可视为等价，但正式国家项目中必须使用 CGCS2000。 diff --git a/geo_tools/__init__.py b/geo_tools/__init__.py deleted file mode 100644 index d0739a4..0000000 --- a/geo_tools/__init__.py +++ /dev/null @@ -1,40 +0,0 @@ -""" -geo_tools -~~~~~~~~~ -专业地理信息数据处理工具库。 - -核心依赖：geopandas、shapely、fiona、pyproj。 - -快速开始 --------- ->>> import geo_tools ->>> gdf = geo_tools.read_vector("data/sample/sample_points.geojson") ->>> gdf_proj = geo_tools.reproject(gdf, "EPSG:3857") ->>> print(gdf_proj.crs) - -GDB 读写 --------- ->>> layers = geo_tools.list_gdb_layers("path/to/data.gdb") ->>> gdf = geo_tools.read_gdb("path/to/data.gdb", layer="my_layer") ->>> geo_tools.write_gdb(gdf, "output/result.gdb", layer="result_layer") -""" - -from importlib.metadata import PackageNotFoundError, version - -# ── 版本 ────────────────────────────────────────────────────────────────────── -try: - __version__ = version("geo-tools") -except PackageNotFoundError: - __version__ = "0.1.0-dev" - -# ── 配置 & 日志 ─────────────────────────────────────────────────────────────── -from geo_tools.config.settings import settings -from geo_tools.utils.logger import get_logger, set_global_level - -__all__ = [ - "__version__", - "settings", - # utils - "get_logger", - "set_global_level" -] diff --git a/geo_tools/analysis/__init__.py b/geo_tools/analysis/__init__.py deleted file mode 100644 index ae7beb7..0000000 --- a/geo_tools/analysis/__init__.py +++ /dev/null @@ -1,25 +0,0 @@ -"""geo_tools.analysis 包 —— 空间分析层。""" - -from geo_tools.analysis.spatial_ops import ( - buffer_and_overlay, - nearest_features, - overlay, - select_by_location, -) -from geo_tools.analysis.stats import ( - area_weighted_mean, - count_by_polygon, - summarize_attributes, -) - -__all__ = [ - # spatial_ops - "buffer_and_overlay", - "overlay", - "nearest_features", - "select_by_location", - # stats - "area_weighted_mean", - "summarize_attributes", - "count_by_polygon", -] diff --git a/geo_tools/core/__init__.py b/geo_tools/core/__init__.py deleted file mode 100644 index 6562bde..0000000 --- a/geo_tools/core/__init__.py +++ /dev/null @@ -1,80 +0,0 @@ -"""geo_tools.core 包 —— 核心地理处理层。""" - -from geo_tools.core.geometry import ( - buffer_geometry, - bounding_box, - centroid, - contains, - convex_hull, - difference, - distance_between, - explain_validity, - fix_geometry, - intersect, - intersects, - is_valid_geometry, - symmetric_difference, - unary_union, - union, - within, -) -from geo_tools.core.projection import ( - CGCS2000, - CGCS2000_UTM_50N, - WEB_MERCATOR, - WGS84, - crs_to_epsg, - get_crs_info, - suggest_projected_crs, - transform_coordinates, - transform_point, -) -from geo_tools.core.vector import ( - add_area_column, - clip_to_extent, - dissolve_by, - drop_invalid_geometries, - explode_multipart, - reproject, - set_crs, - spatial_join, -) - -__all__ = [ - # geometry - "is_valid_geometry", - "fix_geometry", - "explain_validity", - "buffer_geometry", - "centroid", - "bounding_box", - "convex_hull", - "intersect", - "union", - "difference", - "symmetric_difference", - "unary_union", - "contains", - "within", - "intersects", - "distance_between", - # projection - "WGS84", - "CGCS2000", - "WEB_MERCATOR", - "CGCS2000_UTM_50N", - "get_crs_info", - "crs_to_epsg", - "transform_coordinates", - "transform_point", - "suggest_projected_crs", - # vector - "reproject", - "set_crs", - "clip_to_extent", - "dissolve_by", - "explode_multipart", - "drop_invalid_geometries", - "spatial_join", - "add_area_column", -] diff --git a/geo_tools/core/geometry.py b/geo_tools/core/geometry.py deleted file mode 100644 index 5fcf449..0000000 --- a/geo_tools/core/geometry.py +++ /dev/null @@ -1,169 +0,0 @@ -""" -geo_tools.core.geometry -~~~~~~~~~~~~~~~~~~~~~~~ -基于 Shapely 2.x 的几何运算工具函数。 -""" - -from __future__ import annotations - -from typing import Sequence - -import shapely -from shapely.geometry import ( - LinearRing, - LineString, - MultiLineString, - MultiPoint, - MultiPolygon, - Point, - Polygon, -) -from shapely.geometry.base import BaseGeometry - -from geo_tools.utils.logger import get_logger - -logger = get_logger(__name__) - - -# ── 几何有效性 ──────────────────────────────────────────────────────────────── - -def is_valid_geometry(geom: BaseGeometry | None) -> bool: - """判断几何对象是否有效（非空且通过 Shapely 合法性检查）。""" - if geom is None: - return False - return bool(geom.is_valid and not geom.is_empty) - - -def fix_geometry(geom: BaseGeometry | None) -> BaseGeometry | None: - """尝试修复无效几何。 - - 依次尝试： - 1. ``buffer(0)`` — 适合大多数自相交多边形 - 2. ``make_valid``（Shapely 2.x）— 覆盖更多情形 - - Returns - ------- - BaseGeometry | None - 修复后的几何；无法修复时返回 ``None``。 - """ - if geom is None: - return None - if geom.is_valid: - return geom - - # 方法一：buffer(0) - try: - fixed = geom.buffer(0) - if fixed.is_valid and not fixed.is_empty: - return fixed - except Exception: - pass - - # 方法二：shapely.make_valid（Shapely >= 1.8） - try: - fixed = shapely.make_valid(geom) - if fixed.is_valid and not fixed.is_empty: - return fixed - except Exception: - pass - - logger.warning("无法修复几何：%r", geom.geom_type) - return None - - -def explain_validity(geom: BaseGeometry) -> str: - """返回 Shapely 对该几何的有效性说明（英文）。""" - from shapely.validation import explain_validity as _explain - return _explain(geom) - - -# ── 基础几何运算 ─────────────────────────────────────────────────────────────── - -def buffer_geometry( - geom: BaseGeometry, - distance: float, - cap_style: int = 1, - join_style: int = 1, - resolution: int = 16, -) -> BaseGeometry: - """对几何对象执行缓冲区运算。 - - Parameters - ---------- - geom: - 输入几何。 - distance: - 缓冲距离（单位与 CRS 一致；地理坐标系单位为度）。 - cap_style: - 端头样式：1=圆形，2=平头，3=方头（仅线要素有效）。 - join_style: - 转角样式：1=圆角，2=斜角，3=尖角。 - resolution: - 圆弧逼近精度（段数），默认 16。 - """ - return geom.buffer(distance, cap_style=cap_style, join_style=join_style, resolution=resolution) - - -def centroid(geom: BaseGeometry) -> Point: - """返回几何的质心点。""" - return geom.centroid - - -def bounding_box(geom: BaseGeometry) -> Polygon: - """返回几何的最小外接矩形（BBOX）为多边形。""" - from shapely.geometry import box - return box(*geom.bounds) - - -def convex_hull(geom: BaseGeometry) -> BaseGeometry: - """返回几何的凸包。""" - return geom.convex_hull - - -# ── 集合运算 ────────────────────────────────────────────────────────────────── - -def intersect(geom_a: BaseGeometry, geom_b: BaseGeometry) -> BaseGeometry: - """返回两几何的交集。""" - return geom_a.intersection(geom_b) - - -def union(geom_a: BaseGeometry, geom_b: BaseGeometry) -> BaseGeometry: - """返回两几何的并集。""" - return geom_a.union(geom_b) - - -def difference(geom_a: BaseGeometry, geom_b: BaseGeometry) -> BaseGeometry: - """返回 ``geom_a`` 减去 ``geom_b`` 的差集。""" - return geom_a.difference(geom_b) - - -def symmetric_difference(geom_a: BaseGeometry, geom_b: BaseGeometry) -> BaseGeometry: - """返回两几何的对称差集（异或）。""" - return geom_a.symmetric_difference(geom_b) - - -def unary_union(geoms: Sequence[BaseGeometry]) -> BaseGeometry: - """将多个几何合并为一个（等同于逐一 union）。""" - return shapely.unary_union(list(geoms)) - - -# ── 空间关系判断 ─────────────────────────────────────────────────────────────── - -def contains(geom_a: BaseGeometry, geom_b: BaseGeometry) -> bool: - """判断 ``geom_a`` 是否完全包含 ``geom_b``。""" - return bool(geom_a.contains(geom_b)) - - -def within(geom_a: BaseGeometry, geom_b: BaseGeometry) -> bool: - """判断 ``geom_a`` 是否完全在 ``geom_b`` 内。""" - return bool(geom_a.within(geom_b)) - - -def intersects(geom_a: BaseGeometry, geom_b: BaseGeometry) -> bool: - """判断两几何是否相交（含边界接触）。""" - return bool(geom_a.intersects(geom_b)) - - -def distance_between(geom_a: BaseGeometry, geom_b: BaseGeometry) -> float: - """计算两几何间的最小距离（单位与 CRS 一致）。""" - return float(geom_a.distance(geom_b)) diff --git a/geo_tools/io/__init__.py b/geo_tools/io/__init__.py deleted file mode 100644 index ff9c62c..0000000 --- a/geo_tools/io/__init__.py +++ /dev/null @@ -1,31 +0,0 @@ -"""geo_tools.io 包 —— 数据读写层。""" - -from geo_tools.io.readers import ( - list_gdb_layers, - list_gpkg_layers, - read_csv_points, - read_gdb, - read_gpkg, - read_vector, -) -from geo_tools.io.writers import ( - write_csv, - write_gdb, - write_gpkg, - write_vector, -) - -__all__ = [ - # readers - "read_vector", - "read_gdb", - "list_gdb_layers", - "read_gpkg", - "list_gpkg_layers", - "read_csv_points", - # writers - "write_vector", - "write_gdb", - "write_gpkg", - "write_csv", -] diff --git a/scripts/example_workflow.py b/scripts/example_workflow.py index 8d7533b..f30a4eb 100644 --- a/scripts/example_workflow.py +++ b/scripts/example_workflow.py @@ -17,8 +17,8 @@ from pathlib import Path import sys sys.path.insert(0, str(Path(__file__).parent.parent)) -import geo_tools -from geo_tools.utils.logger import get_logger +import app +from app.utils.logger import get_logger logger = get_logger("example_workflow") @@ -29,37 +29,37 @@ OUTPUT_DIR.mkdir(exist_ok=True) def main() -> None: logger.info("=" * 60) - logger.info("geo_tools 端到端工作流示例 v%s", geo_tools.__version__) + logger.info("geo_tools 端到端工作流示例 v%s", app.__version__) logger.info("=" * 60) # ── 1. 读取示例点数据 ────────────────────────────────────────── logger.info("\n[步骤 1] 读取示例点数据（GeoJSON）") - points = geo_tools.read_vector(DATA_DIR / "sample_points.geojson") + points = app.read_vector(DATA_DIR / "sample_points.geojson") logger.info(" 读取完成：%d 条要素，CRS=%s", len(points), points.crs) logger.info(" 字段：%s", list(points.columns)) # ── 2. 读取示例面数据 ────────────────────────────────────────── logger.info("\n[步骤 2] 读取示例区域多边形（GeoJSON）") - regions = geo_tools.read_vector(DATA_DIR / "sample_regions.geojson") + regions = app.read_vector(DATA_DIR / "sample_regions.geojson") logger.info(" 区域列表：%s", regions["name"].tolist()) # ── 3. 数据校验 ─────────────────────────────────────────────── logger.info("\n[步骤 3] 几何有效性校验") - stats = geo_tools.validate_geometry(points) + stats = app.validate_geometry(points) logger.info(" 点数据校验结果：%s", stats) - stats = geo_tools.validate_geometry(regions) + stats = app.validate_geometry(regions) logger.info(" 面数据校验结果：%s", stats) # ── 4. 坐标系信息 ───────────────────────────────────────────── logger.info("\n[步骤 4] 查询 CRS 信息") - crs_info = geo_tools.get_crs_info("EPSG:4326") + crs_info = app.get_crs_info("EPSG:4326") logger.info(" WGS84 信息：%s", crs_info) - proj_crs = geo_tools.suggest_projected_crs(116.4, 39.9) + proj_crs = app.suggest_projected_crs(116.4, 39.9) logger.info(" 北京适合的投影 CRS：%s", proj_crs) # ── 5. 重投影 ───────────────────────────────────────────────── logger.info("\n[步骤 5] 重投影到 Web Mercator（用于可视化）") - points_3857 = geo_tools.reproject(points, "EPSG:3857") + points_3857 = app.reproject(points, "EPSG:3857") logger.info(" 重投影完成：CRS=%s", points_3857.crs) # ── 6. 面积加权均值 ─────────────────────────────────────────── @@ -68,31 +68,31 @@ def main() -> None: points_buffered = points.to_crs("EPSG:3857").copy() points_buffered["geometry"] = points_buffered.geometry.buffer(100_000) # 100km缓冲 points_buffered = points_buffered.to_crs("EPSG:4326") - from geo_tools.analysis.stats import area_weighted_mean + from app.analysis.stats import area_weighted_mean aw_result = area_weighted_mean(points_buffered, value_col="value") logger.info(" 全局面积加权均值：%.4f", aw_result["area_weighted_mean"]) # ── 7. 按位置选择 ───────────────────────────────────────────── logger.info("\n[步骤 7] 按位置选择：筛选华南区域内的城市") hua_nan = regions[regions["name"] == "华南"] - points_in_huanan = geo_tools.select_by_location(points, hua_nan, predicate="intersects") + points_in_huanan = app.select_by_location(points, hua_nan, predicate="intersects") logger.info(" 华南区域内的城市：%s", points_in_huanan["name"].tolist()) # ── 8. 统计汇总 ─────────────────────────────────────────────── logger.info("\n[步骤 8] 属性统计汇总") - from geo_tools.analysis.stats import summarize_attributes + from app.analysis.stats import summarize_attributes summary = summarize_attributes(points, columns=["value"], group_col="category") logger.info(" 按分类汇总：\n%s", summary.to_string(index=False)) # ── 9. 写出结果 ─────────────────────────────────────────────── logger.info("\n[步骤 9] 写出处理结果") out_geojson = OUTPUT_DIR / "result_points_3857.geojson" - geo_tools.write_vector(points_3857, out_geojson) + app.write_vector(points_3857, out_geojson) logger.info(" GeoJSON 写出：%s", out_geojson) out_gpkg = OUTPUT_DIR / "results.gpkg" - geo_tools.write_gpkg(points, out_gpkg, layer="original_points") - geo_tools.write_gpkg(regions, out_gpkg, layer="regions", mode="a") + app.write_gpkg(points, out_gpkg, layer="original_points") + app.write_gpkg(regions, out_gpkg, layer="regions", mode="a") logger.info(" GPKG 写出（2 图层）：%s", out_gpkg) logger.info("\n" + "=" * 60) diff --git a/scripts/其他工具/A耕作层厚度栅格制作_新1.py b/scripts/其他工具/A耕作层厚度栅格制作_新1.py new file mode 100644 index 0000000..5506c35 --- /dev/null +++ b/scripts/其他工具/A耕作层厚度栅格制作_新1.py @@ -0,0 +1,188 @@ +import os +from pathlib import Path +import time + +import geopandas as gpd +from geopandas.io import file +import pandas as pd +import numpy as np + +def assign_gzchd_flexible_v2(soil_prop, point_path, polygon_path, output_path): + print("正在读取数据...") + points = gpd.read_file(point_path) + polygons = gpd.read_file(polygon_path) + + # 1. 坐标系转换 + if points.crs != polygons.crs: + print(f"坐标系不一致，正在转换点数据...") + points = points.to_crs(polygons.crs) + + # 2. 预处理, 判断样点是否存在TZ字段，如果不存在，则用TDLYLX字段代替,并将其转为字符串类型,如果两个字段都不存在，则报错 + if 'TZ' not in points.columns: + if 'TDLYLX' in points.columns: + points['TZ'] = points['TDLYLX'].astype(str).str.strip() + else: + raise ValueError("点要素类中不存在TZ或TDLYLX字段，无法进行匹配！") + else: + points['TZ'] = points['TZ'].astype(str).str.strip() + + polygons['TZ'] = polygons['TZ'].astype(str).str.strip() + + # 确保 GZCHD 是数值类型，避免合并时类型冲突 + points[soil_prop] = pd.to_numeric(points[soil_prop], errors='coerce').fillna(0) + + if soil_prop in polygons.columns: + polygons = polygons.drop(columns=[soil_prop]) + + # 辅助函数：按指定字段分组进行最近点匹配 + def match_by_attribute(poly_gdf, pt_gdf, attr_name, suffix): + if attr_name not in poly_gdf.columns or attr_name not in pt_gdf.columns: + return None, [] + + poly_sub = poly_gdf[poly_gdf[attr_name].notna()].copy() + point_sub = pt_gdf[pt_gdf[attr_name].notna()].copy() + + if poly_sub.empty or point_sub.empty: + return None, [] + + poly_sub[attr_name] = poly_sub[attr_name].astype(str).str.strip() + point_sub[attr_name] = point_sub[attr_name].astype(str).str.strip() + + common_values = set(poly_sub[attr_name].unique()) & set(point_sub[attr_name].unique()) + if not common_values: + return None, [] + + matched_parts = [] + matched_ids = [] + for value in common_values: + poly_part = poly_sub[poly_sub[attr_name] == value].copy() + point_part = point_sub[point_sub[attr_name] == value][[soil_prop, 'geometry']].copy() + if poly_part.empty or point_part.empty: + continue + + matched_part = gpd.sjoin_nearest(poly_part, point_part, how='left', rsuffix=suffix) + matched_part = matched_part[~matched_part.index.duplicated(keep='first')] + if not matched_part.empty: + matched_parts.append(matched_part) + matched_ids.extend(matched_part.index.tolist()) + + if matched_parts: + return pd.concat(matched_parts), matched_ids + return None, [] + + matched_results = [] + matched_indices = [] + + # --- 第一步：按相同 TZ 匹配 --- + print("步骤 1: 正在匹配相同 TZ 的最近点...") + first_matched, first_ids = match_by_attribute(polygons, points, 'TZ', '_p1') + if first_matched is not None and not first_matched.empty: + matched_results.append(first_matched) + matched_indices.extend(first_ids) + + # --- 第二步：按 TS 匹配未匹配面 --- + unmatched_mask = ~polygons.index.isin(matched_indices) + remaining_polygons = polygons[unmatched_mask].copy() + print(f"步骤 2: 正在为 {len(remaining_polygons)} 个要素匹配 TS 最近点...") + + if not remaining_polygons.empty: + ts_matched, ts_ids = match_by_attribute(remaining_polygons, points, 'TS', '_p_ts') + if ts_matched is not None and not ts_matched.empty: + matched_results.append(ts_matched) + matched_indices.extend(ts_ids) + remaining_polygons = polygons[~polygons.index.isin(matched_indices)].copy() + print(f"已匹配 TS: {len(ts_ids)} 个要素，剩余 {len(remaining_polygons)} 个。") + else: + print("未匹配到 TS 类型要素，继续下一步。") + + # --- 第三步：按 YL 匹配未匹配面 --- + if not remaining_polygons.empty: + print(f"步骤 3: 正在为 {len(remaining_polygons)} 个要素匹配 YL 最近点...") + yl_matched, yl_ids = match_by_attribute(remaining_polygons, points, 'YL', '_p_yl') + if yl_matched is not None and not yl_matched.empty: + matched_results.append(yl_matched) + matched_indices.extend(yl_ids) + remaining_polygons = polygons[~polygons.index.isin(matched_indices)].copy() + print(f"已匹配 YL: {len(yl_ids)} 个要素，剩余 {len(remaining_polygons)} 个。") + else: + print("未匹配到 YL 类型要素，继续下一步。") + + # --- 第四步：按 TL 匹配未匹配面 --- + if not remaining_polygons.empty: + print(f"步骤 4: 正在为 {len(remaining_polygons)} 个要素匹配 TL 最近点...") + tl_matched, tl_ids = match_by_attribute(remaining_polygons, points, 'TL', '_p_tl') + if tl_matched is not None and not tl_matched.empty: + matched_results.append(tl_matched) + matched_indices.extend(tl_ids) + remaining_polygons = polygons[~polygons.index.isin(matched_indices)].copy() + print(f"已匹配 TL: {len(tl_ids)} 个要素，剩余 {len(remaining_polygons)} 个。") + else: + print("未匹配到 TL 类型要素，继续全局最近点。") + + else: + print("没有未匹配的面要素，跳过 TS/YL/TL 匹配。") + + # --- 最后：全局最近点匹配剩余面 --- + unmatched_mask = ~polygons.index.isin(matched_indices) + remaining_polygons = polygons[unmatched_mask].copy() + print(f"最后一步: 正在为 {len(remaining_polygons)} 个要素匹配全局最近点...") + + if not remaining_polygons.empty: + point_pool = points[[soil_prop, 'geometry']] + step_final = gpd.sjoin_nearest(remaining_polygons, point_pool, how='left', rsuffix='_p2') + step_final = step_final[~step_final.index.duplicated(keep='first')] + matched_results.append(step_final) + + # --- 第三步：稳健合并 --- + print("正在合并数据...") + # 过滤掉列表中可能的 None 或空 DataFrame，防止 FutureWarning + to_concat = [res for res in matched_results if res is not None and not res.empty] + + if to_concat: + final_gdf = pd.concat(to_concat) + else: + # 如果没有任何匹配结果，返回带空 GZCHD 的原面要素 + final_gdf = polygons.copy() + final_gdf[soil_prop] = 0 + + # --- 4. 清理与保存 --- + # 删除临时列 + cols_to_drop = [ + c for c in final_gdf.columns + if 'index_' in c + or '_p1' in c + or '_p2' in c + or '_p_ts' in c + or '_p_yl' in c + or '_p_tl' in c + ] + final_gdf = final_gdf.drop(columns=cols_to_drop) + + # 强制去重复列名 + final_gdf = final_gdf.loc[:, ~final_gdf.columns.duplicated()] + + # 填充空值并确保类型一致 + if soil_prop in final_gdf.columns: + final_gdf[soil_prop] = final_gdf[soil_prop].fillna(0) + else: + final_gdf[soil_prop] = 0 + + print(f"正在保存结果至: {output_path}") + final_gdf.to_file(output_path, encoding='utf-8') + print("处理完成！") + + return final_gdf + +if __name__ == "__main__": + # 遍历文件夹中所有样点shp文件，并进行处理 + shp_file = r"E:\@三普属性图出图\测试\YXTCHD.shp" # 样点数据文件夹 + dltb_file = r"E:\@三普属性图出图\广西天峨县\@基础数据\土壤类型图\土壤类型图.shp" # 耕地图斑 + output_folder = r"E:\@三普属性图出图\广西天峨县" # 输出文件夹 + + assign_gzchd_flexible_v2( + soil_prop= "YXTCHD", # 耕地层厚度字段名 + point_path=shp_file, # 样点数据 + polygon_path= dltb_file, + output_path=fr'{output_folder}\YXTCHD.shp' # 输出文件 + ) + time.sleep(2) # 防止文件读写冲突 \ No newline at end of file diff --git a/scripts/其他工具/样点剔除统计表格.py b/scripts/其他工具/样点剔除统计表格.py new file mode 100644 index 0000000..c6960a6 --- /dev/null +++ b/scripts/其他工具/样点剔除统计表格.py @@ -0,0 +1,404 @@ +import pandas as pd +import numpy as np +import os +import geopandas as gpd +from openpyxl import Workbook +from openpyxl.styles import Alignment, Font, Border, Side +from openpyxl.utils import get_column_letter + +# 定义指标代码与单位的对应关系 +INDICATOR_UNITS = { + # 基本指标 + 'PH': ('pH', '-'), + 'ECA': ('交换性钙', 'cmol(½Ca²⁺)/kg'), + 'EMG': ('交换性镁', 'cmol(½Mg²⁺)/kg'), + 'TN': ('全氮', 'g/kg'), + 'TP': ('全磷', 'g/kg'), + 'TK': ('全钾', 'g/kg'), + 'AS1': ('有效硫', 'mg/kg'), + 'AB': ('有效硼', 'mg/kg'), + 'AP': ('有效磷', 'mg/kg'), + 'AFE': ('有效铁', 'mg/kg'), + 'ACU': ('有效铜', 'mg/kg'), + 'AZN': ('有效锌', 'mg/kg'), + 'AMN': ('有效锰', 'mg/kg'), + 'OM': ('有机质', 'g/kg'), + 'GZCHD': ('耕层厚度', 'cm'), + 'AK': ('速效钾', 'mg/kg'), + 'CEC': ('阳离子交换量', 'cmol/kg'), + # 特殊指标 - 根据文件名对应字段 + 'FL': ('粉粒', '%'), + 'NL': ('黏粒', '%'), + 'SL': ('砂粒', '%'), + 'TRRZPJZ': ('土壤容重', 'g/cm³'), + 'TRZD': ('土壤质地', '分类'), + # 其他可能指标 + 'AMO': ('有效钼', 'mg/kg'), + 'TSE': ('全硒', 'mg/kg'), + 'YXTCHD': ('有效土层厚度', 'cm') +} + +# 文件名到字段的映射 +FILENAME_TO_FIELD = { + '粉粒': 'FL', + '黏粒': 'NL', + '砂粒': 'SL', + '表层容重': 'TRRZPJZ', + '土壤质地十二级分类': 'TRZD', + '双江县YXTCHD': 'YXTCHD' +} + +# 扩展字段别名映射，支持更多pH字段名 +FIELD_ALIASES = { + 'PH': ['pH', 'PH', 'ph'], # 支持pH的各种大小写形式 + 'ECA': ['交换性钙', 'ECA'], + 'EMG': ['交换性镁', 'EMG'], + 'TN': ['全氮', 'TN'], + 'TP': ['全磷', 'TP'], + 'TK': ['全钾', 'TK'], + 'AS1': ['有效硫', 'AS1'], + 'AB': ['有效硼', 'AB'], + 'AP': ['有效磷', 'AP'], + 'AFE': ['有效铁', 'AFE'], + 'ACU': ['有效铜', 'ACU'], + 'AZN': ['有效锌', 'AZN'], + 'AMN': ['有效锰', 'AMN'], + 'OM': ['有机质', 'OM'], + 'GZCHD': ['耕层厚度', 'GZCHD'], + 'AK': ['速效钾', 'AK'], + 'CEC': ['阳离子交换量', 'CEC'], + 'FL': ['粉粒', 'FL'], + 'NL': ['黏粒', 'NL'], + 'SL': ['砂粒', 'SL'], + 'TRRZPJZ': ['土壤容重', 'TRRZPJZ'], + 'TRZD': ['土壤质地', 'TRZD'], + 'AMO': ['有效钼', 'AMO'], + 'TSE': ['全硒', 'TSE'], + 'YXTCHD': ['有效土层厚度', 'YXTCHD'] +} + + +def find_shapefiles(folder_path): + """在文件夹中递归查找所有的Shapefile文件""" + shapefiles = [] + + for root, dirs, files in os.walk(folder_path): + for file in files: + if file.lower().endswith('.shp'): + shapefiles.append(os.path.join(root, file)) + + return shapefiles + + +def read_shapefile_data(shapefile_path): + """读取Shapefile数据并返回属性表""" + try: + print(f" 读取Shapefile: {os.path.basename(shapefile_path)}") + gdf = gpd.read_file(shapefile_path, encoding='utf-8') + + print(f" 要素数量: {len(gdf)}") + print(f" 属性字段: {list(gdf.columns)}") + + return gdf + except Exception as e: + print(f" 读取Shapefile失败: {e}") + try: + gdf = gpd.read_file(shapefile_path, encoding='gbk') + print(f" 使用GBK编码成功读取") + return gdf + except: + return None + + +def get_indicator_data(gdf, filename): + """从GeoDataFrame中获取指标数据，使用统一字段匹配逻辑""" + indicator_data = {} + + basename = os.path.basename(filename).replace('.shp', '') + + # 1. 首先尝试文件名映射 + target_field = None + if basename in FILENAME_TO_FIELD: + target_field = FILENAME_TO_FIELD[basename] + if target_field in gdf.columns: + indicator_data[target_field] = gdf[target_field] + print(f" 通过文件名映射找到字段: {target_field}") + else: + # 尝试通过别名查找 + for indicator_code in INDICATOR_UNITS.keys(): + if target_field == indicator_code: + for alias in FIELD_ALIASES.get(indicator_code, []): + if alias in gdf.columns: + indicator_data[indicator_code] = gdf[alias] + print(f" 通过文件名映射+别名找到字段: {alias} -> {indicator_code}") + break + + # 2. 如果没有通过文件名找到，尝试直接匹配所有指标和别名 + if not indicator_data: + for indicator_code in INDICATOR_UNITS.keys(): + # 先尝试直接匹配指标代码 + if indicator_code in gdf.columns: + indicator_data[indicator_code] = gdf[indicator_code] + print(f" 直接匹配字段: {indicator_code}") + continue + + # 再尝试匹配别名 + aliases = FIELD_ALIASES.get(indicator_code, []) + for alias in aliases: + if alias in gdf.columns: + indicator_data[indicator_code] = gdf[alias] + print(f" 通过别名匹配: {alias} -> {indicator_code}") + break + + # 3. 额外检查：如果文件名包含特定关键词，尝试匹配 + if not indicator_data: + filename_lower = basename.lower() + for indicator_code, (chinese_name, unit) in INDICATOR_UNITS.items(): + if indicator_code.lower() in filename_lower or chinese_name in filename_lower: + # 尝试匹配指标代码或中文名 + if indicator_code in gdf.columns: + indicator_data[indicator_code] = gdf[indicator_code] + print(f" 通过文件名关键词匹配: {indicator_code}") + break + elif chinese_name in gdf.columns: + indicator_data[indicator_code] = gdf[chinese_name] + print(f" 通过文件名关键词匹配中文名: {chinese_name} -> {indicator_code}") + break + + return indicator_data + + +def get_combined_stats_from_folder(folder_path, folder_name="数据"): + """从文件夹中所有shapefile合并统计指定指标""" + shapefiles = find_shapefiles(folder_path) + + if not shapefiles: + print(f" 未找到Shapefile文件") + return pd.DataFrame() + + print(f" 找到 {len(shapefiles)} 个Shapefile文件") + + all_data = {code: [] for code in INDICATOR_UNITS.keys()} + + for i, shp_file in enumerate(shapefiles, 1): + print(f"\n [{i}] 处理文件: {os.path.basename(shp_file)}") + gdf = read_shapefile_data(shp_file) + + if gdf is not None: + indicator_data = get_indicator_data(gdf, shp_file) + + for indicator_code, data_series in indicator_data.items(): + if indicator_code in all_data: + # 转换为数值类型，处理可能的非数值数据 + try: + data_series = pd.to_numeric(data_series, errors='coerce') + valid_data = data_series.dropna() + if len(valid_data) > 0: + all_data[indicator_code].extend(valid_data.tolist()) + print(f" 提取 {indicator_code}: {len(valid_data)} 个值") + except Exception as e: + print(f" 处理 {indicator_code} 数据时出错: {e}") + + # 计算每个指标的合并统计 + stats_list = [] + + for indicator_code, (chinese_name, unit) in INDICATOR_UNITS.items(): + data_list = all_data.get(indicator_code, []) + if not data_list: + continue + + data_series = pd.Series(data_list) + # 过滤极端值（可选，根据实际需求调整） + data_series = data_series[(data_series >= 0) | pd.isna(data_series)] + + if len(data_series) == 0: + continue + + # 关键修复1：计算总体标准差（ddof=0），而不是默认的样本标准差（ddof=1） + std_dev = data_series.std(ddof=0) + mean_val = data_series.mean() + + # 关键修复2：优化变异系数计算 + if abs(mean_val) < 1e-8: # 均值接近0时 + cv_value = 0.0 + else: + # CV = (标准差 / 均值) * 100，保留2位小数 + cv_value = round((std_dev / mean_val) * 100, 2) + + stats = { + '指标代码': indicator_code, + '指标': chinese_name, + '单位': unit, + '样点数': int(len(data_series)), + 'Min': round(float(data_series.min()), 2), + 'Max': round(float(data_series.max()), 2), + 'Mean': round(float(mean_val), 2), + 'Std': round(float(std_dev), 2), # 使用总体标准差 + 'CV': cv_value + } + stats_list.append(stats) + print(f" 统计 {chinese_name}({indicator_code}): {len(data_series)} 个样点") + + if stats_list: + stats_df = pd.DataFrame(stats_list) + stats_df = stats_df.sort_values('指标') + print(f"\n 总共统计到 {len(stats_df)} 个指标") + return stats_df + + print(" 未找到任何指标数据") + return pd.DataFrame() + + +def create_statistics_excel(before_folder, after_folder, output_path): + """创建融合的统计表格，在剔除后表格前加一列剔除前样点数和剔除样点数""" + workbook = Workbook() + + # 移除默认sheet + if 'Sheet' in workbook.sheetnames: + default_sheet = workbook['Sheet'] + workbook.remove(default_sheet) + + # 定义样式 + thin_border = Border( + left=Side(style='thin'), + right=Side(style='thin'), + top=Side(style='thin'), + bottom=Side(style='thin') + ) + + print("=" * 60) + print("开始分析样点数据") + print("=" * 60) + + # 分析剔除前数据 + before_stats = None + if os.path.exists(before_folder): + print(f"\n[1] 分析剔除前数据:") + print(f"文件夹路径: {before_folder}") + before_stats = get_combined_stats_from_folder(before_folder, "剔除前") + if not before_stats.empty: + print(f"✓ 剔除前统计完成: {len(before_stats)} 个指标") + else: + print("✗ 剔除前未找到指定指标数据") + before_stats = None + else: + print(f"✗ 剔除前文件夹不存在: {before_folder}") + before_stats = None + + # 分析剔除后数据 + if os.path.exists(after_folder): + print(f"\n[2] 分析剔除后数据:") + print(f"文件夹路径: {after_folder}") + after_stats = get_combined_stats_from_folder(after_folder, "剔除后") + + if not after_stats.empty: + # 创建融合的统计工作表 + sheet_combined = workbook.create_sheet(title="样点统计") + + # 新的表头：指标, 单位, 剔除前样点数, 剔除样点数, 剔除后样点数, Min, Max, Mean, Std, CV + combined_headers = ['指标', '单位', '剔除前样点数', '剔除样点数', '剔除后样点数', 'Min', 'Max', 'Mean', + 'Std', 'CV'] + for col, header in enumerate(combined_headers, 1): + cell = sheet_combined.cell(row=1, column=col, value=header) + cell.alignment = Alignment(horizontal='center', vertical='center') + cell.font = Font(bold=True) + cell.border = thin_border + + # 写入数据 + for row_idx, (index, after_row) in enumerate(after_stats.iterrows(), start=2): + # 查找对应的剔除前数据 + before_sample_count = 0 + before_row = None + if before_stats is not None: + # 首先尝试通过指标代码匹配 + matching_rows = before_stats[before_stats['指标代码'] == after_row['指标代码']] + if not matching_rows.empty: + before_row = matching_rows.iloc[0] + else: + # 如果指标代码匹配失败，尝试通过指标名称匹配 + matching_rows = before_stats[before_stats['指标'] == after_row['指标']] + if not matching_rows.empty: + before_row = matching_rows.iloc[0] + + if before_row is not None: + before_sample_count = int(before_row['样点数']) + + # 计算剔除样点数 + after_sample_count = int(after_row['样点数']) + abnormal_count = max(0, before_sample_count - after_sample_count) + + # 写入数据 + sheet_combined.cell(row=row_idx, column=1, value=after_row['指标']) # 指标 + sheet_combined.cell(row=row_idx, column=2, value=after_row['单位']) # 单位 + sheet_combined.cell(row=row_idx, column=3, value=before_sample_count) # 剔除前样点数 + sheet_combined.cell(row=row_idx, column=4, value=abnormal_count) # 剔除样点数 + sheet_combined.cell(row=row_idx, column=5, value=after_sample_count) # 剔除后样点数 + sheet_combined.cell(row=row_idx, column=6, value=after_row['Min']) # Min + sheet_combined.cell(row=row_idx, column=7, value=after_row['Max']) # Max + sheet_combined.cell(row=row_idx, column=8, value=after_row['Mean']) # Mean + sheet_combined.cell(row=row_idx, column=9, value=after_row['Std']) # Std + sheet_combined.cell(row=row_idx, column=10, value=after_row['CV']) # CV + + # 设置所有单元格的样式 + for col_idx in range(1, 11): + cell = sheet_combined.cell(row=row_idx, column=col_idx) + cell.alignment = Alignment(horizontal='center', vertical='center') + cell.border = thin_border + + # 如果剔除了样点，高亮显示剔除样点数列 + if abnormal_count > 0: + cell = sheet_combined.cell(row=row_idx, column=4) # 剔除样点数列 + cell.font = Font(bold=True, color="FF0000") # 红色加粗 + + # 调整列宽 + combined_column_widths = { + '指标': 15, + '单位': 12, + '剔除前样点数': 12, + '剔除样点数': 12, + '剔除后样点数': 12, + 'Min': 10, + 'Max': 10, + 'Mean': 10, + 'Std': 10, + 'CV': 10 + } + + for col_idx, col_name in enumerate(combined_headers, 1): + column_letter = get_column_letter(col_idx) + if col_name in combined_column_widths: + sheet_combined.column_dimensions[column_letter].width = combined_column_widths[col_name] + + print(f"\n✓ 融合统计完成: {len(after_stats)} 个指标") + + # 输出匹配信息 + if before_stats is not None: + print(f" 剔除前找到 {len(before_stats)} 个指标") + print(f" 剔除后找到 {len(after_stats)} 个指标") + print( + f" 成功匹配 {len([i for i in range(2, len(after_stats) + 2) if sheet_combined.cell(row=i, column=3).value > 0])} 个指标的剔除前数据") + else: + print("✗ 剔除后未找到指定指标数据") + sheet_combined = workbook.create_sheet(title="样点统计") + sheet_combined.cell(row=1, column=1, value="未找到指定指标数据") + else: + print(f"✗ 剔除后文件夹不存在: {after_folder}") + sheet_combined = workbook.create_sheet(title="样点统计") + sheet_combined.cell(row=1, column=1, value="剔除后文件夹不存在") + + # 保存文件 + workbook.save(output_path) + print(f"\n" + "=" * 60) + print(f"文件保存成功: {output_path}") + print("=" * 60) + + +# ================ 使用示例 ================ +if __name__ == "__main__": + # 方式1: 处理单个样点数据对 + before_folder = r"D:\a陆平\1.5实验数据\12月新版实验室数据\云南省实验室数据成果1218\永仁县20260127" + after_folder = r"D:\a陆平\1.5实验数据\12月新版实验室数据\云南省实验室数据成果1218\永仁县20260127剔除后" + output_path = r"D:\a陆平\1.5实验数据\12月新版实验室数据\云南省实验室数据成果1218\永仁县样点统计结果.xlsx" + + # 执行 + create_statistics_excel(before_folder, after_folder, output_path) \ No newline at end of file diff --git a/scripts/其他工具/面积加权均值_华南.py b/scripts/其他工具/面积加权均值_华南.py new file mode 100644 index 0000000..23177db --- /dev/null +++ b/scripts/其他工具/面积加权均值_华南.py @@ -0,0 +1,1050 @@ +# 生成完整的exactextract面积加权计算Python脚本 +# -*- coding: utf-8 -*- +""" +土壤属性栅格数据面积加权统计脚本 +基于exactextract库实现多属性栅格的面积加权平均值计算 +最终输出格式与土壤属性图斑数据表一致 +""" + +import re +import sys + +import exactextract as ee +import geopandas as gpd +import rasterio +import pandas as pd +import numpy as np +from pathlib import Path +import warnings +from pyproj import CRS +warnings.filterwarnings('ignore') + +project_root = Path(__file__).parent.parent.parent +sys.path.insert(0, str(project_root)) + +import app + + +def init_logger(): + """初始化日志输出，便于跟踪处理过程""" + import logging + logging.basicConfig( + level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + datefmt='%Y-%m-%d %H:%M:%S' + ) + return logging.getLogger(__name__) + +def validate_data(vector_path, raster_files): + """ + 验证输入数据的有效性 + :param vector_path: 矢量图斑文件路径 + :param raster_files: 栅格文件字典 + :return: 验证结果（布尔值） + """ + logger = init_logger() + + # 验证矢量文件是否存在 + # if not Path(vector_path).exists(): + # logger.error(f"矢量文件不存在：{vector_path}") + # return False + + # 验证栅格文件是否存在 + for attr_name, raster_path in raster_files.items(): + if not Path(raster_path).exists(): + logger.error(f"栅格文件不存在：{attr_name} -> {raster_path}") + return False + + # 验证矢量文件格式 + try: + # gdf = gpd.read_file(vector_path) + gdf = app.readers.read_vector(vector_path) + if gdf.geometry.type.unique()[0] != 'Polygon' and gdf.geometry.type.unique()[0] != 'MultiPolygon': + logger.error("矢量文件必须是Polygon类型（面要素）") + return False + except Exception as e: + logger.error(f"矢量文件读取失败：{str(e)}") + return False + + # 验证栅格文件格式 + try: + test_raster = next(iter(raster_files.values())) + with rasterio.open(test_raster) as src: + if src.count != 1: + logger.error("每个栅格文件必须是单波段（每个属性单独一个栅格）") + return False + except Exception as e: + logger.error(f"栅格文件读取失败：{str(e)}") + return False + + logger.info("所有输入数据验证通过") + return True + +def extract_epsg(crs_obj): + """ + 终极EPSG提取器：从正常方法到暴力正则 + """ + if crs_obj is None: + return None + + # 1. 尝试常规方法 + epsg = crs_obj.to_epsg() + if epsg: + return str(epsg) + + # 2. 降低匹配严格度尝试 (pyproj 内置的容错方法) + epsg_loose = crs_obj.to_epsg(min_confidence=25) + if epsg_loose: + return str(epsg_loose) + + # 3. 终极必杀：正则匹配 WKT 文本 + try: + wkt_str = crs_obj.to_wkt() + # 匹配文本中的 AUTHORITY["EPSG","4523"] 或 ID["EPSG", 4523] + matches = re.findall(r'(?:AUTHORITY|ID)\["EPSG",\s*"?(\d+)"?\]', wkt_str) + if matches: + # WKT中会有多个EPSG(比如单位、椭球体)，最后一个通常才是整个坐标系的EPSG + return matches[-1] + except Exception as e: + pass + + return None + +def standardize_crs(gdf, raster_path): + """ + 标准化矢量与栅格的坐标参考系统（CRS） + """ + logger = init_logger() # 请确保你有 init_logger + + with rasterio.open(raster_path) as src: + # 兼容旧版本 rasterio 的读取方式 + raster_crs = CRS.from_user_input(src.crs) + + vector_crs = CRS.from_user_input(gdf.crs) + + # 使用终极提取器获取 EPSG + v_epsg = extract_epsg(vector_crs) + r_epsg = extract_epsg(raster_crs) + + print(f"当前矢量CRS EPSG：{v_epsg}") + print(f"目标栅格CRS EPSG：{r_epsg}") + + # 判断逻辑：先看数学是否相等，再看强制提取的EPSG编号是否相等 + is_same_crs = False + if vector_crs.equals(raster_crs): + is_same_crs = True + elif v_epsg and r_epsg and v_epsg == r_epsg: + is_same_crs = True + + # 执行转换（若真的不同才转） + if not is_same_crs: + logger.warning("矢量与栅格CRS不一致，正在进行转换...") + gdf = gdf.to_crs(raster_crs) + new_epsg = extract_epsg(CRS.from_user_input(gdf.crs)) + logger.info(f"CRS转换完成，新矢量CRS EPSG：{new_epsg}") + else: + logger.info(f"矢量与栅格CRS一致 (均判定为 EPSG:{v_epsg})，无需转换。") + + return gdf + +def calculate_area_weighted_stats(vector_path, raster_files, output_path): + """ + 核心函数：计算面积加权统计值 + :param vector_path: 矢量图斑文件路径 + :param raster_files: 栅格文件字典（键：属性名，值：栅格路径） + :param output_path: 结果输出路径（Excel文件） + :return: 统计结果DataFrame + """ + logger = init_logger() + logger.info("开始面积加权统计计算") + + # 1. 加载矢量数据 + logger.info(f"加载矢量数据：{vector_path}") + gdf = app.readers.read_vector(vector_path) + + # 2. 标准化CRS + test_raster = next(iter(raster_files.values())) + gdf = standardize_crs(gdf, test_raster) + + # 中断程序 + # sys.exit(0) + + # 3. 初始化结果DataFrame（保留矢量中的关键属性） + logger.info("初始化结果数据结构") + # 基础字段列表（与土壤属性图斑表格式对齐） + # TODO + base_fields = ["BSM","YSDM","TBBH","DLBM", "DLMC","TL", "YL", "TS", "TZ","ZLDWDM","ZLDWMC","TBMJ","geometry"] + + # 检查矢量中是否包含必要字段，若不包含则创建空字段 + result_gdf = gpd.GeoDataFrame() + for field in base_fields: + if field in gdf.columns: + result_gdf[field] = gdf[field] + else: + result_gdf[field] = np.nan + logger.warning(f"矢量中缺少'{field}'字段，将生成空值") + + # 设置几何列以保持GeoDataFrame + result_gdf = gpd.GeoDataFrame(result_gdf, geometry=gdf.geometry) + + # 4. 计算图斑面积（转换为亩） + # logger.info("计算图斑面积（单位：亩）") + # # 计算平方米面积（根据CRS单位自动适应） + # gdf["area_sqm"] = gdf.geometry.area + # # 转换为亩（1亩 ≈ 666.6667平方米） + # result_gdf["面积亩"] = gdf["area_sqm"] * 0.0015 + # # 保留6位小数，与示例数据格式一致 + # result_gdf["面积亩"] = result_gdf["面积亩"].round(6) + + # 5. 对每个土壤属性进行面积加权平均计算 + logger.info("开始处理土壤属性栅格（面积加权平均）") + for attr_idx, (attr_name, raster_path) in enumerate(raster_files.items(), 1): + total_attrs = len(raster_files) + logger.info(f"处理进度：{attr_idx}/{total_attrs} - 属性：{attr_name}") + + try: + # 使用exactextract计算面积加权平均 + # weights="area"：按矢量与栅格的交集面积进行加权,如果是TRZD等分类属性，则不使用面积加权，直接计算众数 + if attr_name in ["TRZD"]: + stats = ee.exact_extract( + raster_path, + gdf, + ["majority"], # 计算众数 + output="pandas" # 输出为DataFrame格式 + ) + else: + stats = ee.exact_extract( + raster_path, + gdf, + ["mean"], # 计算平均值 + output="pandas" # 输出为DataFrame格式 + ) + + # 将统计结果添加到结果DataFrame + if stats is None: + logger.warning(f"{attr_name}计算结果为空，可能无交集区域") + result_gdf[attr_name] = np.nan + continue + else: + # 确保 stats 为 pandas.DataFrame，以便使用字符串索引 + if not isinstance(stats, pd.DataFrame): + stats = pd.DataFrame(stats) + + # 保留4位小数，确保数据精度 + if attr_name in ["TRZD"]: + result_gdf[attr_name] = stats["majority"] + else: + result_gdf[attr_name] = stats["mean"].round(2) + + # 处理可能的空值（无交集区域） + if result_gdf[attr_name].isnull().sum() > 0: + null_count = result_gdf[attr_name].isnull().sum() + logger.warning(f"{attr_name}存在{null_count}个空值（图斑与栅格无交集）") + # 用0填充空值（可根据业务需求调整） + result_gdf[attr_name] = result_gdf[attr_name].fillna(0) + + except Exception as e: + logger.error(f"{attr_name}处理失败：{str(e)}") + # 失败时填充空值，避免整个程序崩溃 + result_gdf[attr_name] = np.nan + + # 6. 添加属性分级字段（根据业务规则实现） + logger.info("添加土壤属性分级字段") + result_gdf = add_attribute_classification(result_gdf) + + # 7. 整理最终字段顺序（与示例表格完全对齐） + logger.info("整理输出字段顺序") + # TODO + final_columns = [ + "BSM","YSDM","TBBH","DLBM", "DLMC","TL", "YL", "TS", "TZ","ZLDWDM","ZLDWMC","TBMJ", + "SL", "FL", "NL", "TRZD", "GZCHD", "TRRZ", "SWXDT", "DBLSFD", + "PH", "CEC", "OM", "TN", "TP", "TK", "AP", "AK", "ECA", "EMG", + "AS1", "AFE", "AMN", "ACU","AZN","AB","AMO", + "SL_DJ", "FL_DJ", "NL_DJ", "TRZD_DJ", "GZCHD_DJ", "TRRZ_DJ", "SWXDT_DJ", "DBLSFD_DJ", + "PH_DJ", "CEC_DJ", "OM_DJ", "TN_DJ", "TP_DJ", "TK_DJ", "AP_DJ", "AK_DJ", "ECA_DJ", "EMG_DJ", + "AS1_DJ", "AFE_DJ", "AMN_DJ", "ACU_DJ","AZN_DJ","AB_DJ","AMO_DJ", + "SL_DJZY", "FL_DJZY", "NL_DJZY", "TRZD_DJZY", "GZCHD_DJZY", "TRRZ_DJZY", "SWXDT_DJZY", "DBLSFD_DJZY", + "PH_DJZY", "CEC_DJZY", "OM_DJZY", "TN_DJZY", "TP_DJZY", "TK_DJZY", "AP_DJZY", "AK_DJZY", "ECA_DJZY", "EMG_DJZY", + "AS1_DJZY", "AFE_DJZY", "AMN_DJZY", "ACU_DJZY","AZN_DJZY","AB_DJZY","AMO_DJZY", + "geometry" # 保留几何列以保持GeoDataFrame + ] + + # 补充缺失的字段（如乡代码等） + for col in final_columns: + if col not in result_gdf.columns: + result_gdf[col] = np.nan + + # 按最终顺序排列字段 + result_gdf = result_gdf[final_columns] + + # 8. 导出结果到Excel + logger.info(f"导出结果到：{output_path}") + # 使用openpyxl引擎支持.xlsx格式 + # result_gdf.to_excel(output_path, index=False, engine="openpyxl") + result_gdf.to_file(output_path, index=False, driver="ESRI Shapefile") # 导出为GeoPackage格式，保留空间信息 + logger.info(f"结果导出完成，共生成{len(result_gdf)}条记录") + + return result_gdf + +def add_attribute_classification(df): + """ + 添加土壤属性分级字段（根据常见土壤分类标准实现） + 可根据实际业务需求调整分级阈值 + :param df: 包含原始属性的DataFrame + :return: 包含分级字段和分级值域的DataFrame + """ + # todo + # 1. 酸碱度分级（pH值） - 按PH标准 + def classify_ph(ph): + """ + 土壤pH值分级（第三次全国土壤普查标准） + 标准等级: + 等级一: 6.0～7.0 + 等级二: 7.0～7.5, 5.5～6.0 + 等级三: 7.5～8.0, 5.0～5.5 + 等级四: 8.0～8.5, 4.5～5.0 + 等级五: ＞8.5, ≤4.5 + """ + if ph > 8.5 or ph <= 4.5: + return 5 # 等级五: ＞8.5, ≤4.5 + elif (8.0 < ph <= 8.5) or (4.5 < ph <= 5.0): + return 4 # 等级四: 8.0～8.5, 4.5～5.0 + elif (7.5 < ph <= 8.0) or (5.0 < ph <= 5.5): + return 3 # 等级三: 7.5～8.0, 5.0～5.5 + elif (7.0 < ph <= 7.5) or (5.5 < ph <= 6.0): + return 2 # 等级二: 7.0～7.5, 5.5～6.0 + elif 6.0 < ph <= 7.0: + return 1 # 等级一: 6.0～7.0 + else: + return None # 异常值 + + # 2. 有机质分级（单位：g/kg） - 按OM标准 + def classify_organic(organic): + """ + 土壤有机质分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞35.0 + 等级二: 25.0～35.0 + 等级三: 15.0～25.0 + 等级四: 10.0～15.0 + 等级五: ≤10.0 + """ + if organic > 35.0: + return 1 # 等级一: ＞35.0 + elif 30.0 < organic <= 35.0: + return 2 # 等级二: 25.0～35.0 + elif 20.0 < organic <= 30.0: + return 3 # 等级三: 15.0～25.0 + elif 10.0 < organic <= 20.0: + return 4 # 等级四: 10.0～15.0 + elif organic <= 10.0: + return 5 # 等级五: ≤10.0 + else: + return None # 异常值 + + # 3. 阳离子交换量分级（单位：cmol/kg） - 按CEC标准 + def classify_cation(cation): + """ + 土壤阳离子交换量分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞20.0 + 等级二: 15.0～20.0 + 等级三: 10.0～15.0 + 等级四: 5.0～10.0 + 等级五: ≤5.0 + """ + if cation > 20.0: + return 1 # 等级一: ＞20.0 + elif 15.0 < cation <= 20.0: + return 2 # 等级二: 15.0～20.0 + elif 10.0 < cation <= 15.0: + return 3 # 等级三: 10.0～15.0 + elif 5.0 < cation <= 10.0: + return 4 # 等级四: 5.0～10.0 + elif cation <= 5.0: + return 5 # 等级五: ≤5.0 + else: + return None # 异常值 + + # 4. 有效磷分级（单位：mg/kg） - 按AP标准 + def classify_available_p(p): + """ + 土壤有效磷分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞40.0 + 等级二: 20.0～40.0 + 等级三: 10.0～20.0 + 等级四: 5.0～10.0 + 等级五: ≤5.0 + """ + if p > 40.0: + return 1 # 等级一: ＞40.0 + elif 20.0 < p <= 40.0: + return 2 # 等级二: 25.0～40.0 + elif 10.0 < p <= 20.0: + return 3 # 等级三: 15.0～25.0 + elif 5.0 < p <= 10.0: + return 4 # 等级四: 5.0～15.0 + elif p <= 5.0: + return 5 # 等级五: ≤5.0 + else: + return None # 异常值 + + # 5. 速效钾分级（单位：mg/kg） - 按AK标准 + def classify_available_k(k): + """ + 土壤速效钾分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞150 + 等级二: 100～150 + 等级三: 75～100 + 等级四: 50～75 + 等级五: ≤50 + """ + if k > 150: + return 1 # 等级一: ＞150 + elif 100 < k <= 150: + return 2 # 等级二: 100～150 + elif 75 < k <= 100: + return 3 # 等级三: 75～100 + elif 50 < k <= 75: + return 4 # 等级四: 50～75 + elif k <= 50: + return 5 # 等级五: ≤50 + else: + return None # 异常值 + + # 6. 耕层厚度分级（单位：cm） - 按GZCHD标准 + def classify_soil_depth(depth): + """ + 土壤耕作层厚度分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞25.0 + 等级二: 20.0～25.0 + 等级三: 15.0～20.0 + 等级四: 10.0～15.0 + 等级五: ≤10.0 + """ + if depth > 25.0: + return 1 # 等级一: ＞25.0 + elif 20.0 < depth <= 25.0: + return 2 # 等级二: 20.0～25.0 + elif 15.0 < depth <= 20.0: + return 3 # 等级三: 15.0～20.0 + elif 10.0 < depth <= 15.0: + return 4 # 等级四: 10.0～15.0 + elif depth <= 10.0: + return 5 # 等级五: ≤10.0 + else: + return None # 异常值 + + # 7. 土壤容重分级（单位：g/cm³） - 按TRRZ标准 + def classify_bulk_density(density): + """ + 土壤容重分级（第三次全国土壤普查标准） + 标准等级: + 等级一: 1.00～1.20 + 等级二: 1.20～1.30 + 等级三: 1.30～1.40, 0.90～1.00 + 等级四: 1.40～1.50 + 等级五: ＞1.50, ≤0.90 + """ + if 1.00 < density <= 1.20: + return 1 # 等级一: 1.00～1.20 + elif 1.20 < density <= 1.30: + return 2 # 等级二: 1.20～1.30 + elif (1.30 < density <= 1.40) or (0.90 < density <= 1.00): + return 3 # 等级三: 1.30～1.40, 0.90～1.00 + elif 1.40 < density <= 1.50: + return 4 # 等级四: 1.40～1.50 + elif density > 1.50 or density <= 0.90: + return 5 # 等级五: ＞1.50, ≤0.90 + else: + return None # 异常值 + + # 8. 全氮分级（单位：g/kg） - 按TN标准 + def classify_total_n(n): + """ + 土壤全氮分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞2.00 + 等级二: 1.50～2.00 + 等级三: 1.00～1.50 + 等级四: 0.50～1.00 + 等级五: ≤0.50 + """ + if n > 2.00: + return 1 # 等级一: ＞2.00 + elif 1.50 < n <= 2.00: + return 2 # 等级二: 1.50～2.00 + elif 1.00 < n <= 1.50: + return 3 # 等级三: 1.00～1.50 + elif 0.50 < n <= 1.00: + return 4 # 等级四: 0.50～1.00 + elif n <= 0.50: + return 5 # 等级五: ≤0.50 + else: + return None # 异常值 + + # 9. 全磷分级（单位：g/kg） - 按TP标准 + def classify_total_p(p): + """ + 土壤全磷分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞1.50 + 等级二: 1.00～1.50 + 等级三: 0.60～1.00 + 等级四: 0.40～0.60 + 等级五: ≤0.40 + """ + if p > 1.50: + return 1 # 等级一: ＞1.50 + elif 1.00 < p <= 1.50: + return 2 # 等级二: 1.00～1.50 + elif 0.60 < p <= 1.00: + return 3 # 等级三: 0.60～1.00 + elif 0.40 <= p <= 0.60: + return 4 # 等级四: 0.40～0.60 + elif p <= 0.40: + return 5 # 等级五: ≤0.40 + else: + return None # 异常值 + + # 10. 全钾分级（单位：g/kg） - 按TK标准 + def classify_total_k(k): + """ + 土壤全钾分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞20.0 + 等级二: 15.0～20.0 + 等级三: 10.0～15.0 + 等级四: 5.0～10.0 + 等级五: ≤5.0 + """ + if k > 20.0: + return 1 # 等级一: ＞20.0 + elif 15.0 < k <= 20.0: + return 2 # 等级二: 15.0～20.0 + elif 10.0 < k <= 15.0: + return 3 # 等级三: 10.0～15.0 + elif 5.0 < k <= 10.0: + return 4 # 等级四: 5.0～10.0 + elif k <= 5.0: + return 5 # 等级五: ≤5.0 + else: + return None # 异常值 + + # 11. 有效铁分级（单位：mg/kg） + def classify_available_fe(fe): + """ + 土壤有效铁分级（第三次全国土壤普查标准） + + 参数: + fe: 有效铁含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if fe <= 3.0: + return 5 # 等级五: ≤3.0 + elif 3.0 < fe <= 8.0: + return 4 # 等级四: 3.0～8.0 + elif 8.0 < fe <= 10.0: + return 3 # 等级三: 8.0～10.0 + elif 10.0 < fe <= 20.0: + return 2 # 等级二: 10.0～20.0 + else: # fe > 20.0 + return 1 # 等级一: ＞20.0 + + # 12. 有效锌分级（单位：mg/kg） + def classify_available_zn(zn): + """ + 土壤有效锌分级（第三次全国土壤普查标准） + + 参数: + zn: 有效锌含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if zn <= 0.50: + return 5 # 等级五: ≤0.50 + elif 0.50 < zn <= 1.00: + return 4 # 等级四: 0.50～1.00 + elif 1.00 < zn <= 2.00: + return 3 # 等级三: 1.00～2.00 + elif 2.00 < zn <= 3.00: + return 2 # 等级二: 2.00～3.00 + else: # zn > 3.00 + return 1 # 等级一: ＞3.00 + + # 13. 有效锰分级（单位：mg/kg） + def classify_available_mn(mn): + """ + 土壤有效锰分级（第三次全国土壤普查标准） + + 参数: + mn: 有效锰含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if mn <= 5.0: + return 5 # 等级五: ≤5.0 + elif 5.0 < mn <= 10.0: + return 4 # 等级四: 5.0～10.0 + elif 10.0 < mn <= 20.0: + return 3 # 等级三: 10.0～20.0 + elif 20.0 < mn <= 30.0: + return 2 # 等级二: 20.0～30.0 + else: # mn > 30.0 + return 1 # 等级一: ＞30.0 + + # 14. 有效铜分级（单位：mg/kg） + def classify_available_cu(cu): + """ + 土壤有效铜分级（第三次全国土壤普查标准） + + 参数: + cu: 有效铜含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if cu <= 0.20: + return 5 # 等级五: ≤0.20 + elif 0.20 < cu <= 0.50: + return 4 # 等级四: 0.20～0.50 + elif 0.50 < cu <= 1.00: + return 3 # 等级三: 0.50～1.00 + elif 1.00 < cu <= 1.80: + return 2 # 等级二: 1.00～1.80 + else: # cu > 1.80 + return 1 # 等级一: ＞1.80 + + # 15. 有效硼分级（单位：mg/kg） + def classify_available_b(b): + """ + 土壤有效硼分级（第三次全国土壤普查标准） + + 参数: + b: 有效硼含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if b <= 0.20: + return 5 # 等级五: ≤0.20 + elif 0.20 < b <= 0.50: + return 4 # 等级四: 0.20～0.50 + elif 0.50 < b <= 1.00: + return 3 # 等级三: 0.50～1.00 + elif 1.00 < b <= 2.00: + return 2 # 等级二: 1.00～2.00 + else: # b > 2.00 + return 1 # 等级一: ＞2.00 + + # 16. 有效钼分级（单位：mg/kg） + def classify_available_mo(mo): + """ + 土壤有效钼分级（第三次全国土壤普查标准） + + 参数: + mo: 有效钼含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if mo <= 0.05: + return 5 # 等级五: ≤0.05 + elif 0.05 < mo <= 0.10: + return 4 # 等级四: 0.05～0.10 + elif 0.10 < mo <= 0.15: + return 3 # 等级三: 0.10～0.15 + elif 0.15 < mo <= 0.20: + return 2 # 等级二: 0.15～0.20 + else: # mo > 0.20 + return 1 # 等级一: ＞0.20 + + # 17. 有效硫分级（单位：mg/kg） + def classify_available_s(s): + """ + 土壤有效硫分级（第三次全国土壤普查标准） + + 参数: + s: 有效硫含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if s <= 10.0: + return 5 # 等级五: ≤10.0 + elif 10.0 < s <= 20.0: + return 4 # 等级四: 10.0～20.0 + elif 20.0 < s <= 30.0: + return 3 # 等级三: 20.0～30.0 + elif 30.0 < s <= 40.0: + return 2 # 等级二: 30.0～40.0 + else: # s > 40.0 + return 1 # 等级一: ＞40.0 + + # 18. 交换性钙分级（单位：cmol(½Ca²⁺)/kg） + def classify_exchangeable_ca(ca): + """ + 土壤交换性钙分级（第三次全国土壤普查标准） + + 参数: + ca: 交换性钙含量 (cmol(½Ca²⁺)/kg) + + 返回: + 分级等级 (1-5) + """ + if ca <= 0.25: + return 5 # 等级五: ≤0.25 + elif 0.25 < ca <= 1.00: + return 4 # 等级四: 0.25～1.00 + elif 1.00 < ca <= 2.50: + return 3 # 等级三: 1.00～2.50 + elif 2.50 < ca <= 4.99: + return 2 # 等级二: 2.50～4.99 + else: # ca > 4.99 + return 1 # 等级一: ＞4.99 + + # 19. 交换性镁分级（单位：cmol(½Mg²⁺)/kg） + def classify_exchangeable_mg(mg): + """ + 土壤交换性镁分级（第三次全国土壤普查标准） + + 参数: + mg: 交换性镁含量 (cmol(½Mg²⁺)/kg) + + 返回: + 分级等级 (1-5) + """ + if mg <= 0.21: + return 5 # 等级五: ≤0.21 + elif 0.21 < mg <= 0.41: + return 4 # 等级四: 0.21～0.41 + elif 0.41 < mg <= 0.82: + return 3 # 等级三: 0.41～0.82 + elif 0.82 < mg <= 1.64: + return 2 # 等级二: 0.82～1.64 + else: # mg > 1.64 + return 1 # 等级一: ＞1.64 + + # 20. 全硒分级（单位：mg/kg） + def classify_total_se(se): + """ + 土壤全硒分级（第三次全国土壤普查标准） + + 参数: + se: 全硒含量 (mg/kg) + + 返回: + 分级等级 (1-4) + """ + if se <= 0.17: + return 4 # 等级四: ≤0.17 + elif 0.17 < se <= 0.40: + return 3 # 等级三: 0.17～0.40 + elif 0.40 < se <= 3.00: + return 2 # 等级二: 0.40～3.00 + else: # se > 3.00 + return 1 # 等级一: ＞3.00 + + # 21. 粉粒含量分级（单位：%） + def classify_silt(silt): + """ + 土壤粉粒含量分级（第三次全国土壤普查标准） + + 参数: + silt: 粉粒含量 (%) + + 返回: + 分级等级 (1-5) + """ + if silt > 75: + return 5 # 等级五: ＞75 + elif 45 < silt <= 75: + return 4 # 等级四: 45～75 + elif 30 < silt <= 45: + return 3 # 等级三: 30～45 + elif 15 < silt <= 30: + return 2 # 等级二: 15～30 + else: # silt <= 15 + return 1 # 等级一: ≤15 + + # 22. 黏粒含量分级（单位：%） + def classify_clay(clay): + """ + 土壤黏粒含量分级（第三次全国土壤普查标准） + + 参数: + clay: 黏粒含量 (%) + + 返回: + 分级等级 (1-5) + """ + if clay > 65: + return 5 # 等级五: ＞65 + elif 45 < clay <= 65: + return 4 # 等级四: 45～65 + elif 25 < clay <= 45: + return 3 # 等级三: 25～45 + elif 15 < clay <= 25: + return 2 # 等级二: 15～25 + else: # clay <= 15 + return 1 # 等级一: ≤15 + + # 23. 砂粒含量分级（单位：%） + def classify_sand(sand): + """ + 土壤砂粒含量分级（第三次全国土壤普查标准） + + 参数: + sand: 砂粒含量 (%) + + 返回: + 分级等级 (1-5) + """ + if sand > 85: + return 5 # 等级五: ＞85 + elif 55 < sand <= 85: + return 4 # 等级四: 55～85 + elif 40 < sand <= 55: + return 3 # 等级三: 40～55 + elif 30 < sand <= 40: + return 2 # 等级二: 30～40 + else: # sand <= 30 + return 1 # 等级一: ≤30 + + # 24. 有效土层厚度分级（单位：cm） + def classify_yxtchd(depth): + """ + 土壤有效土层厚度分级（第三次全国土壤普查标准） + + 参数: + depth: 有效土层厚度 (cm) + + 返回: + 分级等级 (1-5) + """ + if depth <= 40: + return 5 # 等级五: ≤40 + elif 40 < depth <= 60: + return 4 # 等级四: 40～60 + elif 60 < depth <= 80: + return 3 # 等级三: 60～80 + elif 80 < depth <= 100: + return 2 # 等级二: 80～100 + else: # depth > 100 + return 1 # 等级一: ＞100 + + # 25. 土壤质地 + def classify_trzd(trzd): + """ + 土壤质地 + + 参数: + trzd: 土壤质地分类（1-5） + + 返回: + 分级等级 (1-5) + """ + trzd = round(trzd, 0) + return trzd + # if trzd == 5: + # return 5 # 等级五: ≤40 + # elif trzd == 4: + # return 4 # 等级四: 40～60 + # elif trzd == 3: + # return 3 # 等级三: 60～80 + # elif trzd == 2: + # return 2 # 等级二: 80～100 + # else: # depth > 100 + # return 1 # 等级一: ＞100 + + # 26. 地表砾石丰度分级（单位：%） + def classify_lsfd(gravel): + """ + 土壤地表砾石丰度分级（第三次全国土壤普查标准） + + 参数: + gravel: 地表砾石丰度 (%) + + 返回: + 分级等级 (1-5) + """ + if gravel <= 0: + return 1 # 等级一: ≤0 + elif 0 < gravel <= 5: + return 2 # 等级二: 0～5 + elif 5 < gravel <= 15: + return 3 # 等级三: 5～15 + elif 15 < gravel <= 50: + return 4 # 等级三: 15～50 + else: # gravel > 50 + return 5 # 等级五: ＞50 + + # 应用分级函数（只处理非空值） + if "PH" in df.columns: + df["PH_DJ"] = df["PH"].apply(lambda x: classify_ph(x) if pd.notna(x) else np.nan) + df["PH_DJZY"] = "等级一: 6.0～7.0,等级二: 7.0～7.5, 5.5～6.0,等级三: 7.5～8.0, 5.0～5.5,等级四: 8.0～8.5, 4.5～5.0,等级五: ＞8.5, ≤4.5" + if "OM" in df.columns: + df["OM_DJ"] = df["OM"].apply(lambda x: classify_organic(x) if pd.notna(x) else np.nan) + df["OM_DJZY"] = "等级一: ＞35,等级二: 30～35,等级三: 20～30,等级四: 10～20,等级五: ≤10" + if "CEC" in df.columns: + df["CEC_DJ"] = df["CEC"].apply(lambda x: classify_cation(x) if pd.notna(x) else np.nan) + df["CEC_DJZY"] = "等级一: ＞20,等级二: 15～20,等级三: 10～15,等级四: 5～10,等级五: ≤5" + if "AP" in df.columns: + df["AP_DJ"] = df["AP"].apply(lambda x: classify_available_p(x) if pd.notna(x) else np.nan) + df["AP_DJZY"] = "等级一: ＞40,等级二: 20～40,等级三: 10～20,等级四: 5～10,等级五: ≤5" + if "AK" in df.columns: + df["AK_DJ"] = df["AK"].apply(lambda x: classify_available_k(x) if pd.notna(x) else np.nan) + df["AK_DJZY"] = "等级一: ＞150,等级二: 100～150,等级三: 75～100,等级四: 50～75,等级五: ≤50" + if "GZCHD" in df.columns: + df["GZCHD_DJ"] = df["GZCHD"].apply(lambda x: classify_soil_depth(x) if pd.notna(x) else np.nan) + df["GZCHD_DJZY"] = "等级一: ＞25,等级二: 20～25,等级三: 15～20,等级四: 10～15,等级五: ≤10" + if "TRRZ" in df.columns: + df["TRRZ_DJ"] = df["TRRZ"].apply(lambda x: classify_bulk_density(x) if pd.notna(x) else np.nan) + df["TRRZ_DJZY"] = "等级一: ＞1.00～1.20,等级二: 1.20～1.30,等级三: 1.30～1.40,0.90～1.0,等级四: 1.40～1.50,等级五: ＞1.50,≤0.90" + if "TN" in df.columns: + df["TN_DJ"] = df["TN"].apply(lambda x: classify_total_n(x) if pd.notna(x) else np.nan) + df["TN_DJZY"] = "等级一: ＞2.0,等级二: 1.5～2.0,等级三: 1.0～1.5,等级四: 0.5～1.0,等级五: ≤0.5" + if "TP" in df.columns: + df["TP_DJ"] = df["TP"].apply(lambda x: classify_total_p(x) if pd.notna(x) else np.nan) + df["TP_DJZY"] = "等级一: ＞1.5,等级二: 1.0～1.5,等级三: 0.6～1.0,等级四: 0.4～0.6,等级五: ≤0.4" + if "TK" in df.columns: + df["TK_DJ"] = df["TK"].apply(lambda x: classify_total_k(x) if pd.notna(x) else np.nan) + df["TK_DJZY"] = "等级一: ＞20.0,等级二: 15.0～20.0,等级三: 10.0～15.0,等级四: 5.0～10.0,等级五: ≤5.0" + if "AFE" in df.columns: + df["AFE_DJ"] = df["AFE"].apply(lambda x: classify_available_fe(x) if pd.notna(x) else np.nan) + df["AFE_DJZY"] = "等级一: ＞20,等级二: 10～20,等级三: 8～10,等级四: 3～8,等级五: ≤3" + if "AZN" in df.columns: + df["AZN_DJ"] = df["AZN"].apply(lambda x: classify_available_zn(x) if pd.notna(x) else np.nan) + df["AZN_DJZY"] = "等级一: ＞3.0,等级二: 2.0～3.0,等级三: 1.0～2.0,等级四: 0.5～1.0,等级五: ≤0.5" + if "AMN" in df.columns: + df["AMN_DJ"] = df["AMN"].apply(lambda x: classify_available_mn(x) if pd.notna(x) else np.nan) + df["AMN_DJZY"] = "等级一: ＞30.0,等级二: 20.0～30.0,等级三: 10.0～20.0,等级四: 5.0～10.0,等级五: ≤5.0" + if "ACU" in df.columns: + df["ACU_DJ"] = df["ACU"].apply(lambda x: classify_available_cu(x) if pd.notna(x) else np.nan) + df["ACU_DJZY"] = "等级一: ＞1.80,等级二: 1.00～1.80,等级三: 0.50～1.00,等级四: 0.2～0.5,等级五: ≤0.2" + if "AB" in df.columns: + df["AB_DJ"] = df["AB"].apply(lambda x: classify_available_b(x) if pd.notna(x) else np.nan) + df["AB_DJZY"] = "等级一: ＞2.0,等级二: 1.0～2.0,等级三: 0.5～1.0,等级四: 0.2～0.5,等级五: ≤0.2" + if "AMO" in df.columns: + df["AMO_DJ"] = df["AMO"].apply(lambda x: classify_available_mo(x) if pd.notna(x) else np.nan) + df["AMO_DJZY"] = "等级一: ＞0.20,等级二: 0.15～0.20,等级三: 0.10～0.15,等级四: 0.05～0.10,等级五: ≤0.05" + if "AS1" in df.columns: + df["AS1_DJ"] = df["AS1"].apply(lambda x: classify_available_s(x) if pd.notna(x) else np.nan) + df["AS1_DJZY"] = "等级一: ＞40.0,等级二: 30.0～40.0,等级三: 20.0～30.0,等级四: 10.0～20.0,等级五: ≤10.0" + if "ECA" in df.columns: + df["ECA_DJ"] = df["ECA"].apply(lambda x: classify_exchangeable_ca(x) if pd.notna(x) else np.nan) + df["ECA_DJZY"] = "等级一: ＞20,等级二: 15～30.0,等级三: 10.0～15.0,等级四: 5～10.0,等级五: ≤5.0" + if "EMG" in df.columns: + df["EMG_DJ"] = df["EMG"].apply(lambda x: classify_exchangeable_mg(x) if pd.notna(x) else np.nan) + df["EMG_DJZY"] = "等级一: ＞1.64,等级二: 0.82～1.64,等级三: 0.41～0.82,等级四: 0.21～0.41,等级五: ≤0.21" + if "TSE" in df.columns: + df["TSE_DJ"] = df["TSE"].apply(lambda x: classify_total_se(x) if pd.notna(x) else np.nan) + df["TSE_DJZY"] = "等级一: ＞3.0,等级二: 0.4～3.0,等级三: 0.17～0.4,等级四:≤0.17" + if "FL" in df.columns: + df["FL_DJ"] = df["FL"].apply(lambda x: classify_silt(x) if pd.notna(x) else np.nan) + df["FL_DJZY"] = "等级一: ≤15,等级二: 15～30,等级三: 30～45,等级四: 45～75,等级五: ＞75" + if "NL" in df.columns: + df["NL_DJ"] = df["NL"].apply(lambda x: classify_clay(x) if pd.notna(x) else np.nan) + df["NL_DJZY"] = "等级一: ≤15,等级二: 15～25,等级三: 25～45,等级四: 45～65,等级五: ＞65" + if "SL" in df.columns: + df["SL_DJ"] = df["SL"].apply(lambda x: classify_sand(x) if pd.notna(x) else np.nan) + df["SL_DJZY"] = "等级一: ≤30,等级二: 30～40,等级三: 40～55,等级四: 55～85,等级五: ＞85" + if "YXTCHD" in df.columns: + df["YXTCHD_DJ"] = df["YXTCHD"].apply(lambda x: classify_yxtchd(x) if pd.notna(x) else np.nan) + df["YXTCHD_DJZY"] = "等级一: ＞100,等级二: 80～100,等级三: 60～80,等级四: 40～60,等级五: ≤40" + if "TRZD" in df.columns: + df["TRZD_DJ"] = df["TRZD"].apply(lambda x: classify_trzd(x) if pd.notna(x) else np.nan) + df["TRZD_DJZY"] = "粉(砂)质黏壤土: 1,粉(砂)质黏土: 2,粉(砂)质壤土: 3,黏壤土: 4,黏土: 5,壤土: 6,壤质黏土: 7,砂土及壤质砂土: 8,砂质黏壤土: 9,砂质黏土: 10,砂质壤土: 11,重黏土: 12" + if "DBLSFD" in df.columns: + df["DBLSFD_DJ"] = df["DBLSFD"].apply(lambda x: classify_lsfd(x) if pd.notna(x) else np.nan) + df["DBLSFD_DJZY"] = "等级一: =0,等级二: 0～5,等级三: 5.0～15.0,等级四: 15～50,等级五: ≥50" + + return df + +def main(): + """ + 主函数：程序入口 + 用户需根据实际情况修改以下参数 + """ + logger = init_logger() + logger.info("="*50) + logger.info("土壤属性栅格面积加权统计程序启动") + logger.info("="*50) + + # -------------------------- + # 用户配置区域（必须修改！） + # -------------------------- + # 1. 矢量图斑文件路径（支持Shapefile、GeoPackage等格式） + # TODO + VECTOR_PATH = r"F:\@user\missum\Documents\ArcGIS\Projects\MyProject\MyProject.gdb\西畴县耕园林草其他_SpatialJoin" # 示例："D:/data/土壤图斑.shp" + + # 2. 土壤属性栅格文件配置（键：属性名称，值：栅格文件路径） + # 注意：属性名称必须与最终表格列名一致 + # TODO + RASTER_FILES = { + # "GZCHD": r"E:\@三普\云南丘北县\基础数据\三普栅格\GZCHD.tif", + # "TRRZ": r"E:\@三普\云南丘北县\基础数据\三普栅格\TRRZ.tif", + # "SL": r"E:\@三普\云南丘北县\基础数据\三普栅格\SL.tif", + # "FL": r"E:\@三普\云南丘北县\基础数据\三普栅格\FL.tif", + # "NL": r"E:\@三普\云南丘北县\基础数据\三普栅格\NL.tif", + # "PH": r"E:\@三普\云南丘北县\基础数据\三普栅格\PH.tif", + # "CEC": r"E:\@三普\云南丘北县\基础数据\三普栅格\CEC.tif", + # "OM": r"E:\@三普\云南丘北县\基础数据\三普栅格\OM.tif", + # "TN": r"E:\@三普\云南丘北县\基础数据\三普栅格\TN.tif", + # "TP": r"E:\@三普\云南丘北县\基础数据\三普栅格\TP.tif", + # "TK": r"E:\@三普\云南丘北县\基础数据\三普栅格\TK.tif", + # "AP": r"E:\@三普\云南丘北县\基础数据\三普栅格\AP.tif", + # "AK": r"E:\@三普\云南丘北县\基础数据\三普栅格\AK.tif", + # "AFE": r"E:\@三普\云南丘北县\基础数据\三普栅格\AFE.tif", + # "AZN": r"E:\@三普\云南丘北县\基础数据\三普栅格\AZN.tif", + # "AMN": r"E:\@三普\云南丘北县\基础数据\三普栅格\AMN.tif", + # "ACU": r"E:\@三普\云南丘北县\基础数据\三普栅格\ACU.tif", + # "AB": r"E:\@三普\云南丘北县\基础数据\三普栅格\AB.tif", + # "AMO": r"E:\@三普\云南丘北县\基础数据\三普栅格\AMO.tif", + # "AS1": r"E:\@三普\云南丘北县\基础数据\三普栅格\AS1.tif", + # "ECA": r"E:\@三普\云南丘北县\基础数据\三普栅格\ECA.tif", + # "EMG": r"E:\@三普\云南丘北县\基础数据\三普栅格\EMG.tif", + # "YXTCHD": r"E:\@三普\云南丘北县\基础数据\三普栅格\YXTCHD.tif", + "TRZD": r"E:\@三普\云南西畴县\基础数据\三普栅格\投影后\TRZD.tif", + # "DBLSFD": r"E:\@三普\云南丘北县\基础数据\三普栅格\LSFD.tif", + # "TSE": r"E:\@三普\云南丘北县\基础数据\三普栅格\TSE.tif" + } + + # 3. 结果输出路径（Excel文件） + OUTPUT_PATH = r"E:\@三普\@临时文件夹\西畴县TRSX11.shp" # 示例："D:/result/结果.shp" + + # -------------------------- + # 程序执行流程（无需修改） + # -------------------------- + try: + # 1. 数据验证 + if not validate_data(VECTOR_PATH, RASTER_FILES): + logger.error("数据验证失败，程序终止") + return + + # 2. 执行面积加权统计 + result_gdf = calculate_area_weighted_stats(VECTOR_PATH, RASTER_FILES, OUTPUT_PATH) + + # 3. 显示结果预览 + # logger.info("\\n结果预览（前3行）：") + # print(result_gdf.head(3).to_string(index=False)) + + logger.info("\\n" + "="*50) + logger.info("程序执行完成！") + logger.info(f"结果文件：{OUTPUT_PATH}") + logger.info("="*50) + + except Exception as e: + logger.error(f"程序执行出错：{str(e)}", exc_info=True) + logger.error("程序异常终止") + +if __name__ == "__main__": + # 启动主程序 + main() diff --git a/scripts/其他工具/面积加权均值_西南.py b/scripts/其他工具/面积加权均值_西南.py new file mode 100644 index 0000000..76c2207 --- /dev/null +++ b/scripts/其他工具/面积加权均值_西南.py @@ -0,0 +1,932 @@ +# 生成完整的exactextract面积加权计算Python脚本 +# -*- coding: utf-8 -*- +""" +土壤属性栅格数据面积加权统计脚本 +基于exactextract库实现多属性栅格的面积加权平均值计算 +最终输出格式与土壤属性图斑数据表一致 +""" + +import exactextract as ee +import geopandas as gpd +import rasterio +import pandas as pd +import numpy as np +from pathlib import Path +import warnings +warnings.filterwarnings('ignore') + +def init_logger(): + """初始化日志输出，便于跟踪处理过程""" + import logging + logging.basicConfig( + level=logging.INFO, + format='%(asctime)s - %(levelname)s - %(message)s', + datefmt='%Y-%m-%d %H:%M:%S' + ) + return logging.getLogger(__name__) + +def validate_data(vector_path, raster_files): + """ + 验证输入数据的有效性 + :param vector_path: 矢量图斑文件路径 + :param raster_files: 栅格文件字典 + :return: 验证结果（布尔值） + """ + logger = init_logger() + + # 验证矢量文件是否存在 + if not Path(vector_path).exists(): + logger.error(f"矢量文件不存在：{vector_path}") + return False + + # 验证栅格文件是否存在 + for attr_name, raster_path in raster_files.items(): + if not Path(raster_path).exists(): + logger.error(f"栅格文件不存在：{attr_name} -> {raster_path}") + return False + + # 验证矢量文件格式 + try: + gdf = gpd.read_file(vector_path) + if gdf.geometry.type.unique()[0] != 'Polygon': + logger.error("矢量文件必须是Polygon类型（面要素）") + return False + except Exception as e: + logger.error(f"矢量文件读取失败：{str(e)}") + return False + + # 验证栅格文件格式 + try: + test_raster = next(iter(raster_files.values())) + with rasterio.open(test_raster) as src: + if src.count != 1: + logger.error("每个栅格文件必须是单波段（每个属性单独一个栅格）") + return False + except Exception as e: + logger.error(f"栅格文件读取失败：{str(e)}") + return False + + logger.info("所有输入数据验证通过") + return True + +def standardize_crs(gdf, raster_path): + """ + 标准化矢量与栅格的坐标参考系统（CRS） + :param gdf: 矢量GeoDataFrame + :param raster_path: 任意一个栅格文件路径（用于获取目标CRS） + :return: 标准化后的GeoDataFrame + """ + logger = init_logger() + + # 获取栅格CRS + with rasterio.open(raster_path) as src: + raster_crs = src.crs + raster_crs_str = src.crs.to_string() + + # 获取矢量CRS + vector_crs = gdf.crs + vector_crs_str = gdf.crs.to_string() + + logger.info(f"当前矢量CRS：{vector_crs_str}") + logger.info(f"目标栅格CRS：{raster_crs_str}") + + # 若CRS不一致，进行转换 + if vector_crs != raster_crs: + logger.warning("矢量与栅格CRS不一致，正在进行转换...") + gdf = gdf.to_crs(raster_crs) + logger.info(f"CRS转换完成，新矢量CRS：{gdf.crs.to_string()}") + + return gdf + +def calculate_area_weighted_stats(vector_path, raster_files, output_path): + """ + 核心函数：计算面积加权统计值 + :param vector_path: 矢量图斑文件路径 + :param raster_files: 栅格文件字典（键：属性名，值：栅格路径） + :param output_path: 结果输出路径（Excel文件） + :return: 统计结果DataFrame + """ + logger = init_logger() + logger.info("开始面积加权统计计算") + + # 1. 加载矢量数据 + logger.info(f"加载矢量数据：{vector_path}") + gdf = gpd.read_file(vector_path) + + # 2. 标准化CRS + test_raster = next(iter(raster_files.values())) + gdf = standardize_crs(gdf, test_raster) + + # 3. 初始化结果DataFrame（保留矢量中的关键属性） + logger.info("初始化结果数据结构") + # 基础字段列表（与土壤属性图斑表格式对齐） + # TODO + base_fields = ["FID","DM","XZM","QSDWDM","QSDWMC","TL", "YL", "TS", "TZ", "DLBM", "DLMC"] + + # 检查矢量中是否包含必要字段，若不包含则创建空字段 + result_df = pd.DataFrame() + for field in base_fields: + if field in gdf.columns: + result_df[field] = gdf[field] + else: + result_df[field] = np.nan + logger.warning(f"矢量中缺少'{field}'字段，将生成空值") + + # 4. 计算图斑面积（转换为亩） + logger.info("计算图斑面积（单位：亩）") + # 计算平方米面积（根据CRS单位自动适应） + gdf["area_sqm"] = gdf.geometry.area + # 转换为亩（1亩 ≈ 666.6667平方米） + result_df["面积亩"] = gdf["area_sqm"] * 0.0015 + # 保留6位小数，与示例数据格式一致 + result_df["面积亩"] = result_df["面积亩"].round(6) + + # 5. 对每个土壤属性进行面积加权平均计算 + logger.info("开始处理土壤属性栅格（面积加权平均）") + for attr_idx, (attr_name, raster_path) in enumerate(raster_files.items(), 1): + total_attrs = len(raster_files) + logger.info(f"处理进度：{attr_idx}/{total_attrs} - 属性：{attr_name}") + + try: + # 使用exactextract计算面积加权平均 + # weights="area"：按矢量与栅格的交集面积进行加权 + stats = ee.exact_extract( + raster_path, + gdf, + ["mean"], # 计算平均值 + output="pandas" # 输出为DataFrame格式 + ) + + # 将统计结果添加到结果DataFrame + if stats is None: + logger.warning(f"{attr_name}计算结果为空，可能无交集区域") + result_df[attr_name] = np.nan + continue + else: + # 确保 stats 为 pandas.DataFrame，以便使用字符串索引 + if not isinstance(stats, pd.DataFrame): + stats = pd.DataFrame(stats) + + # 保留4位小数，确保数据精度 + result_df[attr_name] = stats["mean"].round(4) + + # 处理可能的空值（无交集区域） + if result_df[attr_name].isnull().sum() > 0: + null_count = result_df[attr_name].isnull().sum() + logger.warning(f"{attr_name}存在{null_count}个空值（图斑与栅格无交集）") + # 用0填充空值（可根据业务需求调整） + result_df[attr_name] = result_df[attr_name].fillna(0) + + except Exception as e: + logger.error(f"{attr_name}处理失败：{str(e)}") + # 失败时填充空值，避免整个程序崩溃 + result_df[attr_name] = np.nan + + # 6. 添加属性分级字段（根据业务规则实现） + logger.info("添加土壤属性分级字段") + result_df = add_attribute_classification(result_df) + + # 7. 整理最终字段顺序（与示例表格完全对齐） + logger.info("整理输出字段顺序") + # TODO + final_columns = [ + "FID","DM", "XZM", "QSDWDM", "QSDWMC", + "TL", "YL", "TS", "TZ", "DLBM", "DLMC", + "耕层厚度", "土壤容重", "砂粒", "粉粒", "黏粒", "酸碱度", "阳离子", + "有机质", "全氮", "全磷", "全钾", "有效磷", "速效钾", "有效铁", "有效锰", + "有效铜", "有效锌", "有效硼", "有效钼", "有效硫", "交换性钙", "交换性镁", "全硒", + "有效土层厚度", "土壤质地", + "耕层厚度分级", "土壤容重分级", "砂粒分级", "粉粒分级", "黏粒分级", "酸碱度分级", "阳离子分级", + "有机质分级", "全氮分级", "全磷分级", "全钾分级", "有效铁分级", "速效钾分级", "有效铁分级", "有效锰分级", + "有效铜分级", "有效锌分级", "有效硼分级", "有效钼分级", "有效硫分级", "交换性钙分级", "交换性镁分级", "全硒分级", + "有效土层厚度分级", "土壤质地分级", + "面积亩" + ] + + # 补充缺失的字段（如乡代码等） + for col in final_columns: + if col not in result_df.columns: + result_df[col] = np.nan + + # 按最终顺序排列字段 + result_df = result_df[final_columns] + + # 8. 导出结果到Excel + logger.info(f"导出结果到：{output_path}") + # 使用openpyxl引擎支持.xlsx格式 + result_df.to_excel(output_path, index=False, engine="openpyxl") + logger.info(f"结果导出完成，共生成{len(result_df)}条记录") + + return result_df + +def add_attribute_classification(df): + """ + 添加土壤属性分级字段（根据常见土壤分类标准实现） + 可根据实际业务需求调整分级阈值 + :param df: 包含原始属性的DataFrame + :return: 包含分级字段的DataFrame + """ + # todo + # 1. 酸碱度分级（pH值） - 按PH标准 + def classify_ph(ph): + """ + 土壤pH值分级（第三次全国土壤普查标准） + 标准等级: + 等级一: 6.0～7.0 + 等级二: 7.0～7.5, 5.5～6.0 + 等级三: 7.5～8.0, 5.0～5.5 + 等级四: 8.0～8.5, 4.5～5.0 + 等级五: ＞8.5, ≤4.5 + """ + if ph > 8.5 or ph <= 4.5: + return 5 # 等级五: ＞8.5, ≤4.5 + elif (8.0 < ph <= 8.5) or (4.5 < ph <= 5.0): + return 4 # 等级四: 8.0～8.5, 4.5～5.0 + elif (7.5 < ph <= 8.0) or (5.0 < ph <= 5.5): + return 3 # 等级三: 7.5～8.0, 5.0～5.5 + elif (7.0 < ph <= 7.5) or (5.5 < ph <= 6.0): + return 2 # 等级二: 7.0～7.5, 5.5～6.0 + elif 6.0 < ph <= 7.0: + return 1 # 等级一: 6.0～7.0 + else: + return None # 异常值 + + # 2. 有机质分级（单位：g/kg） - 按OM标准 + def classify_organic(organic): + """ + 土壤有机质分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞35.0 + 等级二: 25.0～35.0 + 等级三: 15.0～25.0 + 等级四: 10.0～15.0 + 等级五: ≤10.0 + """ + if organic > 35.0: + return 1 # 等级一: ＞35.0 + elif 25.0 < organic <= 35.0: + return 2 # 等级二: 25.0～35.0 + elif 15.0 < organic <= 25.0: + return 3 # 等级三: 15.0～25.0 + elif 10.0 < organic <= 15.0: + return 4 # 等级四: 10.0～15.0 + elif organic <= 10.0: + return 5 # 等级五: ≤10.0 + else: + return None # 异常值 + + # 3. 阳离子交换量分级（单位：cmol/kg） - 按CEC标准 + def classify_cation(cation): + """ + 土壤阳离子交换量分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞30.0 + 等级二: 20.0～30.0 + 等级三: 15.0～20.0 + 等级四: 10.0～15.0 + 等级五: ≤10.0 + """ + if cation > 30.0: + return 1 # 等级一: ＞30.0 + elif 20.0 < cation <= 30.0: + return 2 # 等级二: 20.0～30.0 + elif 15.0 < cation <= 20.0: + return 3 # 等级三: 15.0～20.0 + elif 10.0 < cation <= 15.0: + return 4 # 等级四: 10.0～15.0 + elif cation <= 10.0: + return 5 # 等级五: ≤10.0 + else: + return None # 异常值 + + # 4. 有效磷分级（单位：mg/kg） - 按AP标准 + def classify_available_p(p): + """ + 土壤有效磷分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞40.0 + 等级二: 25.0～40.0 + 等级三: 15.0～25.0 + 等级四: 5.0～15.0 + 等级五: ≤5.0 + """ + if p > 40.0: + return 1 # 等级一: ＞40.0 + elif 25.0 < p <= 40.0: + return 2 # 等级二: 25.0～40.0 + elif 15.0 < p <= 25.0: + return 3 # 等级三: 15.0～25.0 + elif 5.0 < p <= 15.0: + return 4 # 等级四: 5.0～15.0 + elif p <= 5.0: + return 5 # 等级五: ≤5.0 + else: + return None # 异常值 + + # 5. 速效钾分级（单位：mg/kg） - 按AK标准 + def classify_available_k(k): + """ + 土壤速效钾分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞150 + 等级二: 100～150 + 等级三: 75～100 + 等级四: 50～75 + 等级五: ≤50 + """ + if k > 150: + return 1 # 等级一: ＞150 + elif 100 < k <= 150: + return 2 # 等级二: 100～150 + elif 75 < k <= 100: + return 3 # 等级三: 75～100 + elif 50 < k <= 75: + return 4 # 等级四: 50～75 + elif k <= 50: + return 5 # 等级五: ≤50 + else: + return None # 异常值 + + # 6. 耕层厚度分级（单位：cm） - 按GZCHD标准 + def classify_soil_depth(depth): + """ + 土壤耕作层厚度分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞25.0 + 等级二: 20.0～25.0 + 等级三: 15.0～20.0 + 等级四: 10.0～15.0 + 等级五: ≤10.0 + """ + if depth > 25.0: + return 1 # 等级一: ＞25.0 + elif 20.0 < depth <= 25.0: + return 2 # 等级二: 20.0～25.0 + elif 15.0 < depth <= 20.0: + return 3 # 等级三: 15.0～20.0 + elif 10.0 < depth <= 15.0: + return 4 # 等级四: 10.0～15.0 + elif depth <= 10.0: + return 5 # 等级五: ≤10.0 + else: + return None # 异常值 + + # 7. 土壤容重分级（单位：g/cm³） - 按TRRZ标准 + def classify_bulk_density(density): + """ + 土壤容重分级（第三次全国土壤普查标准） + 标准等级: + 等级一: 1.10～1.25 + 等级二: 1.25～1.35, 1.00～1.10 + 等级三: 1.35～1.45 + 等级四: 1.45～1.55, 0.90～1.00 + 等级五: ＞1.55, ≤0.90 + """ + if 1.10 < density <= 1.25: + return 1 # 等级一: 1.10～1.25 + elif (1.25 < density <= 1.35) or (1.00 < density <= 1.10): + return 2 # 等级二: 1.25～1.35, 1.00～1.10 + elif 1.35 < density <= 1.45: + return 3 # 等级三: 1.35～1.45 + elif (1.45 < density <= 1.55) or (0.90 < density <= 1.00): + return 4 # 等级四: 1.45～1.55, 0.90～1.00 + elif density > 1.55 or density <= 0.90: + return 5 # 等级五: ＞1.55, ≤0.90 + else: + return None # 异常值 + + # 8. 全氮分级（单位：g/kg） - 按TN标准 + def classify_total_n(n): + """ + 土壤全氮分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞2.00 + 等级二: 1.50～2.00 + 等级三: 1.00～1.50 + 等级四: 0.50～1.00 + 等级五: ≤0.50 + """ + if n > 2.00: + return 1 # 等级一: ＞2.00 + elif 1.50 < n <= 2.00: + return 2 # 等级二: 1.50～2.00 + elif 1.00 < n <= 1.50: + return 3 # 等级三: 1.00～1.50 + elif 0.50 < n <= 1.00: + return 4 # 等级四: 0.50～1.00 + elif n <= 0.50: + return 5 # 等级五: ≤0.50 + else: + return None # 异常值 + + # 9. 全磷分级（单位：g/kg） - 按TP标准 + def classify_total_p(p): + """ + 土壤全磷分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞1.00 + 等级二: 0.80～1.00 + 等级三: 0.60～0.80 + 等级四: 0.40～0.60 + 等级五: ≤0.40 + """ + if p > 1.00: + return 1 # 等级一: ＞1.00 + elif 0.80 < p <= 1.00: + return 2 # 等级二: 0.80～1.00 + elif 0.60 < p <= 0.80: + return 3 # 等级三: 0.60～0.80 + elif 0.40 < p <= 0.60: + return 4 # 等级四: 0.40～0.60 + elif p <= 0.40: + return 5 # 等级五: ≤0.40 + else: + return None # 异常值 + + # 10. 全钾分级（单位：g/kg） - 按TK标准 + def classify_total_k(k): + """ + 土壤全钾分级（第三次全国土壤普查标准） + 标准等级: + 等级一: ＞20.0 + 等级二: 15.0～20.0 + 等级三: 10.0～15.0 + 等级四: 5.0～10.0 + 等级五: ≤5.0 + """ + if k > 20.0: + return 1 # 等级一: ＞20.0 + elif 15.0 < k <= 20.0: + return 2 # 等级二: 15.0～20.0 + elif 10.0 < k <= 15.0: + return 3 # 等级三: 10.0～15.0 + elif 5.0 < k <= 10.0: + return 4 # 等级四: 5.0～10.0 + elif k <= 5.0: + return 5 # 等级五: ≤5.0 + else: + return None # 异常值 + + # 11. 有效铁分级（单位：mg/kg） + def classify_available_fe(fe): + """ + 土壤有效铁分级（第三次全国土壤普查标准） + + 参数: + fe: 有效铁含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if fe <= 3.0: + return 5 # 等级五: ≤3.0 + elif 3.0 < fe <= 5.0: + return 4 # 等级四: 3.0～5.0 + elif 5.0 < fe <= 10.0: + return 3 # 等级三: 5.0～10.0 + elif 10.0 < fe <= 20.0: + return 2 # 等级二: 10.0～20.0 + else: # fe > 20.0 + return 1 # 等级一: ＞20.0 + + # 12. 有效锌分级（单位：mg/kg） + def classify_available_zn(zn): + """ + 土壤有效锌分级（第三次全国土壤普查标准） + + 参数: + zn: 有效锌含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if zn <= 0.20: + return 5 # 等级五: ≤0.20 + elif 0.20 < zn <= 0.50: + return 4 # 等级四: 0.20～0.50 + elif 0.50 < zn <= 1.00: + return 3 # 等级三: 0.50～1.00 + elif 1.00 < zn <= 3.00: + return 2 # 等级二: 1.00～3.00 + else: # zn > 3.00 + return 1 # 等级一: ＞3.00 + + # 13. 有效锰分级（单位：mg/kg） + def classify_available_mn(mn): + """ + 土壤有效锰分级（第三次全国土壤普查标准） + + 参数: + mn: 有效锰含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if mn <= 1.0: + return 5 # 等级五: ≤1.0 + elif 1.0 < mn <= 5.0: + return 4 # 等级四: 1.0～5.0 + elif 5.0 < mn <= 15.0: + return 3 # 等级三: 5.0～15.0 + elif 15.0 < mn <= 30.0: + return 2 # 等级二: 15.0～30.0 + else: # mn > 30.0 + return 1 # 等级一: ＞30.0 + + # 14. 有效铜分级（单位：mg/kg） + def classify_available_cu(cu): + """ + 土壤有效铜分级（第三次全国土壤普查标准） + + 参数: + cu: 有效铜含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if cu <= 0.20: + return 5 # 等级五: ≤0.20 + elif 0.20 < cu <= 0.50: + return 4 # 等级四: 0.20～0.50 + elif 0.50 < cu <= 1.00: + return 3 # 等级三: 0.50～1.00 + elif 1.00 < cu <= 2.00: + return 2 # 等级二: 1.00～2.00 + else: # cu > 2.00 + return 1 # 等级一: ＞2.00 + + # 15. 有效硼分级（单位：mg/kg） + def classify_available_b(b): + """ + 土壤有效硼分级（第三次全国土壤普查标准） + + 参数: + b: 有效硼含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if b <= 0.20: + return 5 # 等级五: ≤0.20 + elif 0.20 < b <= 0.50: + return 4 # 等级四: 0.20～0.50 + elif 0.50 < b <= 0.80: + return 3 # 等级三: 0.50～0.80 + elif 0.80 < b <= 1.00: + return 2 # 等级二: 0.80～1.00 + else: # b > 1.00 + return 1 # 等级一: ＞1.00 + + # 16. 有效钼分级（单位：mg/kg） + def classify_available_mo(mo): + """ + 土壤有效钼分级（第三次全国土壤普查标准） + + 参数: + mo: 有效钼含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if mo <= 0.05: + return 5 # 等级五: ≤0.05 + elif 0.05 < mo <= 0.10: + return 4 # 等级四: 0.05～0.10 + elif 0.10 < mo <= 0.15: + return 3 # 等级三: 0.10～0.15 + elif 0.15 < mo <= 0.20: + return 2 # 等级二: 0.15～0.20 + else: # mo > 0.20 + return 1 # 等级一: ＞0.20 + + # 17. 有效硫分级（单位：mg/kg） + def classify_available_s(s): + """ + 土壤有效硫分级（第三次全国土壤普查标准） + + 参数: + s: 有效硫含量 (mg/kg) + + 返回: + 分级等级 (1-5) + """ + if s <= 10.0: + return 5 # 等级五: ≤10.0 + elif 10.0 < s <= 20.0: + return 4 # 等级四: 10.0～20.0 + elif 20.0 < s <= 30.0: + return 3 # 等级三: 20.0～30.0 + elif 30.0 < s <= 40.0: + return 2 # 等级二: 30.0～40.0 + else: # s > 40.0 + return 1 # 等级一: ＞40.0 + + # 18. 交换性钙分级（单位：cmol(½Ca²⁺)/kg） + def classify_exchangeable_ca(ca): + """ + 土壤交换性钙分级（第三次全国土壤普查标准） + + 参数: + ca: 交换性钙含量 (cmol(½Ca²⁺)/kg) + + 返回: + 分级等级 (1-5) + """ + if ca <= 1.00: + return 5 # 等级五: ≤1.00 + elif 1.00 < ca <= 2.50: + return 4 # 等级四: 1.00～2.50 + elif 2.50 < ca <= 4.99: + return 3 # 等级三: 2.50～4.99 + elif 4.99 < ca <= 7.49: + return 2 # 等级二: 4.99～7.49 + else: # ca > 7.49 + return 1 # 等级一: ＞7.49 + + # 19. 交换性镁分级（单位：cmol(½Mg²⁺)/kg） + def classify_exchangeable_mg(mg): + """ + 土壤交换性镁分级（第三次全国土壤普查标准） + + 参数: + mg: 交换性镁含量 (cmol(½Mg²⁺)/kg) + + 返回: + 分级等级 (1-5) + """ + if mg <= 0.41: + return 5 # 等级五: ≤0.41 + elif 0.41 < mg <= 0.82: + return 4 # 等级四: 0.41～0.82 + elif 0.82 < mg <= 1.23: + return 3 # 等级三: 0.82～1.23 + elif 1.23 < mg <= 1.64: + return 2 # 等级二: 1.23～1.64 + else: # mg > 1.64 + return 1 # 等级一: ＞1.64 + + # 20. 全硒分级（单位：mg/kg） + def classify_total_se(se): + """ + 土壤全硒分级（第三次全国土壤普查标准） + + 参数: + se: 全硒含量 (mg/kg) + + 返回: + 分级等级 (1-4) + """ + if se <= 0.17: + return 4 # 等级四: ≤0.17 + elif 0.17 < se <= 0.40: + return 3 # 等级三: 0.17～0.40 + elif 0.40 < se <= 3.00: + return 2 # 等级二: 0.40～3.00 + else: # se > 3.00 + return 1 # 等级一: ＞3.00 + + # 21. 粉粒含量分级（单位：%） + def classify_silt(silt): + """ + 土壤粉粒含量分级（第三次全国土壤普查标准） + + 参数: + silt: 粉粒含量 (%) + + 返回: + 分级等级 (1-5) + """ + if silt > 75: + return 5 # 等级五: ＞75 + elif 45 < silt <= 75: + return 4 # 等级四: 45～75 + elif 30 < silt <= 45: + return 3 # 等级三: 30～45 + elif 15 < silt <= 30: + return 2 # 等级二: 15～30 + else: # silt <= 15 + return 1 # 等级一: ≤15 + + # 22. 黏粒含量分级（单位：%） + def classify_clay(clay): + """ + 土壤黏粒含量分级（第三次全国土壤普查标准） + + 参数: + clay: 黏粒含量 (%) + + 返回: + 分级等级 (1-5) + """ + if clay > 65: + return 5 # 等级五: ＞65 + elif 45 < clay <= 65: + return 4 # 等级四: 45～65 + elif 25 < clay <= 45: + return 3 # 等级三: 25～45 + elif 15 < clay <= 25: + return 2 # 等级二: 15～25 + else: # clay <= 15 + return 1 # 等级一: ≤15 + + # 23. 砂粒含量分级（单位：%） + def classify_sand(sand): + """ + 土壤砂粒含量分级（第三次全国土壤普查标准） + + 参数: + sand: 砂粒含量 (%) + + 返回: + 分级等级 (1-5) + """ + if sand > 85: + return 5 # 等级五: ＞85 + elif 55 < sand <= 85: + return 4 # 等级四: 55～85 + elif 40 < sand <= 55: + return 3 # 等级三: 40～55 + elif 30 < sand <= 40: + return 2 # 等级二: 30～40 + else: # sand <= 30 + return 1 # 等级一: ≤30 + + # 24. 有效土层厚度分级（单位：cm） + def classify_yxtchd(depth): + """ + 土壤有效土层厚度分级（第三次全国土壤普查标准） + + 参数: + depth: 有效土层厚度 (cm) + + 返回: + 分级等级 (1-5) + """ + if depth <= 40: + return 5 # 等级五: ≤40 + elif 40 < depth <= 60: + return 4 # 等级四: 40～60 + elif 60 < depth <= 80: + return 3 # 等级三: 60～80 + elif 80 < depth <= 100: + return 2 # 等级二: 80～100 + else: # depth > 100 + return 1 # 等级一: ＞100 + + # 25. 土壤质地 + def classify_trzd(trzd): + """ + 土壤质地 + + 参数: + trzd: 土壤质地分类（1-5） + + 返回: + 分级等级 (1-5) + """ + trzd = round(trzd, 0) + if trzd == 5: + return 5 # 等级五: ≤40 + elif trzd == 4: + return 4 # 等级四: 40～60 + elif trzd == 3: + return 3 # 等级三: 60～80 + elif trzd == 2: + return 2 # 等级二: 80～100 + else: # depth > 100 + return 1 # 等级一: ＞100 + + + # 应用分级函数（只处理非空值） + if "酸碱度" in df.columns: + df["酸碱度分级"] = df["酸碱度"].apply(lambda x: classify_ph(x) if pd.notna(x) else np.nan) + if "有机质" in df.columns: + df["有机质分级"] = df["有机质"].apply(lambda x: classify_organic(x) if pd.notna(x) else np.nan) + if "阳离子" in df.columns: + df["阳离子分级"] = df["阳离子"].apply(lambda x: classify_cation(x) if pd.notna(x) else np.nan) + if "有效磷" in df.columns: + df["有效磷分级"] = df["有效磷"].apply(lambda x: classify_available_p(x) if pd.notna(x) else np.nan) + if "速效钾" in df.columns: + df["速效钾分级"] = df["速效钾"].apply(lambda x: classify_available_k(x) if pd.notna(x) else np.nan) + if "耕层厚度" in df.columns: + df["耕层厚度分级"] = df["耕层厚度"].apply(lambda x: classify_soil_depth(x) if pd.notna(x) else np.nan) + if "土壤容重" in df.columns: + df["土壤容重分级"] = df["土壤容重"].apply(lambda x: classify_bulk_density(x) if pd.notna(x) else np.nan) + if "全氮" in df.columns: + df["全氮分级"] = df["全氮"].apply(lambda x: classify_total_n(x) if pd.notna(x) else np.nan) + if "全磷" in df.columns: + df["全磷分级"] = df["全磷"].apply(lambda x: classify_total_p(x) if pd.notna(x) else np.nan) + if "全钾" in df.columns: + df["全钾分级"] = df["全钾"].apply(lambda x: classify_total_k(x) if pd.notna(x) else np.nan) + if "有效铁" in df.columns: + df["有效铁分级"] = df["有效铁"].apply(lambda x: classify_available_fe(x) if pd.notna(x) else np.nan) + if "有效锌" in df.columns: + df["有效锌分级"] = df["有效锌"].apply(lambda x: classify_available_zn(x) if pd.notna(x) else np.nan) + if "有效锰" in df.columns: + df["有效锰分级"] = df["有效锰"].apply(lambda x: classify_available_mn(x) if pd.notna(x) else np.nan) + if "有效铜" in df.columns: + df["有效铜分级"] = df["有效铜"].apply(lambda x: classify_available_cu(x) if pd.notna(x) else np.nan) + if "有效硼" in df.columns: + df["有效硼分级"] = df["有效硼"].apply(lambda x: classify_available_b(x) if pd.notna(x) else np.nan) + if "有效钼" in df.columns: + df["有效钼分级"] = df["有效钼"].apply(lambda x: classify_available_mo(x) if pd.notna(x) else np.nan) + if "有效硫" in df.columns: + df["有效硫分级"] = df["有效硫"].apply(lambda x: classify_available_s(x) if pd.notna(x) else np.nan) + if "交换性钙" in df.columns: + df["交换性钙分级"] = df["交换性钙"].apply(lambda x: classify_exchangeable_ca(x) if pd.notna(x) else np.nan) + if "交换性镁" in df.columns: + df["交换性镁分级"] = df["交换性镁"].apply(lambda x: classify_exchangeable_mg(x) if pd.notna(x) else np.nan) + if "全硒" in df.columns: + df["全硒分级"] = df["全硒"].apply(lambda x: classify_total_se(x) if pd.notna(x) else np.nan) + if "粉粒" in df.columns: + df["粉粒分级"] = df["粉粒"].apply(lambda x: classify_silt(x) if pd.notna(x) else np.nan) + if "黏粒" in df.columns: + df["黏粒分级"] = df["黏粒"].apply(lambda x: classify_clay(x) if pd.notna(x) else np.nan) + if "砂粒" in df.columns: + df["砂粒分级"] = df["砂粒"].apply(lambda x: classify_sand(x) if pd.notna(x) else np.nan) + if "有效土层厚度" in df.columns: + df["有效土层厚度分级"] = df["有效土层厚度"].apply(lambda x: classify_yxtchd(x) if pd.notna(x) else np.nan) + if "土壤质地" in df.columns: + df["土壤质地分级"] = df["土壤质地"].apply(lambda x: classify_trzd(x) if pd.notna(x) else np.nan) + + return df + +def main(): + """ + 主函数：程序入口 + 用户需根据实际情况修改以下参数 + """ + logger = init_logger() + logger.info("="*50) + logger.info("土壤属性栅格面积加权统计程序启动") + logger.info("="*50) + + # -------------------------- + # 用户配置区域（必须修改！） + # -------------------------- + # 1. 矢量图斑文件路径（支持Shapefile、GeoPackage等格式） + # TODO + VECTOR_PATH = r"D:\工作\三普成果编制\出图数据\北海\三普栅格\DL_ALL.shp" # 示例："D:/data/土壤图斑.shp" + + # 2. 土壤属性栅格文件配置（键：属性名称，值：栅格文件路径） + # 注意：属性名称必须与最终表格列名一致 + # TODO + RASTER_FILES = { + "耕层厚度": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\GZCHD.tif", # 示例："D:/data/耕层厚度.tif" + "土壤容重": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\TRRZ.tif", # 示例："D:/data/土壤容重.tif" + "砂粒": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\SL.tif", # 示例："D:/data/砂粒含量.tif" + "粉粒": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\FL.tif", # 示例："D:/data/粉粒含量.tif" + "黏粒": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\NL.tif", # 示例："D:/data/黏粒含量.tif" + "酸碱度": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\PH.tif", # 示例："D:/data/pH值.tif" + "阳离子": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\CEC.tif", # 示例："D:/data/阳离子交换量.tif" + "有机质": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\OM.tif", # 示例："D:/data/有机质含量.tif" + "全氮": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\TN.tif", # 示例："D:/data/全氮含量.tif" + "全磷": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\TP.tif", # 示例："D:/data/全磷含量.tif" + "全钾": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\TK.tif", # 示例："D:/data/全钾含量.tif" + "有效磷": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AP.tif", # 示例："D:/data/有效磷含量.tif" + "速效钾": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AK.tif", # 示例："D:/data/速效钾含量.tif" + # "有效铁": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AFE.tif", # 示例："D:/data/有效铁含量.tif" + "有效锌": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AZN.tif", # 示例："D:/data/有效锌含量.tif" + "有效锰": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AMN.tif", # 示例："D:/data/有效锰含量.tif" + "有效铜": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\ACU.tif", # 示例："D:/data/有效铜含量.tif" + "有效硼": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AB.tif", # 示例："D:/data/有效硼含量.tif" + "有效钼": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AMO.tif", # 示例："D:/data/有效钼含量.tif" + "有效硫": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\AS1.tif", # 示例："D:/data/有效硫含量.tif" + "交换性钙": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\ECA.tif", # 示例："D:/data/交换性钙含量.tif" + "交换性镁": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\EMG.tif", # 示例："D:/data/交换性镁含量.tif" + "全硒": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\TSE.tif", # 示例："D:/data/全硒含量.tif" + "有效土层厚度": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\YXTCHD.tif", # 示例："D:/data/有效土层厚度.tif" + "土壤质地": r"D:\工作\三普成果编制\出图数据\北海\三普栅格\TRZD.tif", # 示例："D:/data/土壤质地.tif" + } + + # 3. 结果输出路径（Excel文件） + OUTPUT_PATH = "土壤属性图斑数据_面积加权结果.xlsx" # 示例："D:/result/结果.xlsx" + + # -------------------------- + # 程序执行流程（无需修改） + # -------------------------- + try: + # 1. 数据验证 + if not validate_data(VECTOR_PATH, RASTER_FILES): + logger.error("数据验证失败，程序终止") + return + + # 2. 执行面积加权统计 + result_df = calculate_area_weighted_stats(VECTOR_PATH, RASTER_FILES, OUTPUT_PATH) + + # 3. 显示结果预览 + logger.info("\\n结果预览（前3行）：") + print(result_df.head(3).to_string(index=False)) + + logger.info("\\n" + "="*50) + logger.info("程序执行完成！") + logger.info(f"结果文件：{OUTPUT_PATH}") + logger.info("="*50) + + except Exception as e: + logger.error(f"程序执行出错：{str(e)}", exc_info=True) + logger.error("程序异常终止") + +if __name__ == "__main__": + # 启动主程序 + main() diff --git a/tests/test1.py b/tests/test1.py index 2b2ece4..397b4e1 100644 --- a/tests/test1.py +++ b/tests/test1.py @@ -7,14 +7,21 @@ from pathlib import Path # 添加项目根目录到路径 sys.path.insert(0, str(Path(__file__).parent.parent)) -import geo_tools +import app -gdb_path = r"E:\@三普\@临时文件夹\临时数据库.gdb" +gdb_path = r"E:\@三普\@临时文件夹\临时数据库.gdb\马关综合后图斑" +shp_path = r"E:\@三普\@临时文件夹\新建文件夹\靖西二普样点\AK.shp" +shp = r"E:\@三普\@临时文件夹\新建文件夹\靖西二普样点" + +gdb = r"E:\@三普\@临时文件夹\新建文件地理数据库.gdb" # 列出图层 -# layers = geo_tools.list_gdb_layers(gdb_path) +# layers = geo_tools.readers.list_gdb_layers(gdb_path) # print(layers) # 读取图层 -gdf = geo_tools.read_gdb(gdb_path, layer="马关综合后图斑") -print(gdf.crs) \ No newline at end of file +gdf = app.readers.read_vector(gdb_path) +print(gdf.crs) +# 获取几何类型 +print(gdf.head()) + diff --git a/tests/test_analysis.py b/tests/test_analysis.py index a86a0b5..269d8cf 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -4,8 +4,8 @@ import pytest import geopandas as gpd from shapely.geometry import Point, Polygon -from geo_tools.analysis.spatial_ops import overlay, select_by_location -from geo_tools.analysis.stats import area_weighted_mean, count_by_polygon, summarize_attributes +from app.analysis.spatial_ops import overlay, select_by_location +from app.analysis.stats import area_weighted_mean, count_by_polygon, summarize_attributes class TestOverlay: diff --git a/tests/test_geometry.py b/tests/test_geometry.py index a1cd50d..a242d8a 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -3,8 +3,8 @@ import pytest from shapely.geometry import LineString, Point, Polygon -import geo_tools -from geo_tools.core.geometry import ( +import app +from app.core.geometry import ( buffer_geometry, bounding_box, centroid, diff --git a/tests/test_io.py b/tests/test_io.py index a628874..859b98b 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -4,8 +4,8 @@ import pytest import geopandas as gpd from pathlib import Path -from geo_tools.io.readers import read_vector, read_gpkg, list_gpkg_layers, read_csv_points -from geo_tools.io.writers import write_vector, write_gpkg, write_csv +from app.io.readers import read_vector, read_gpkg, list_gpkg_layers, read_csv_points +from app.io.writers import write_vector, write_gpkg, write_csv class TestReadVector: @@ -61,7 +61,7 @@ class TestWriteReadRoundtrip: assert "geometry" in df.columns # 存在 WKT 几何列 assert len(df) == len(sample_points_gdf) # 行数一致 # 再用 read_csv_points 以 WKT 模式读回 - from geo_tools.io.readers import _read_csv_vector + from app.io.readers import _read_csv_vector from pathlib import Path gdf_back = _read_csv_vector(Path(out), wkt_col="geometry") assert len(gdf_back) == len(sample_points_gdf) diff --git a/tests/test_proj.py b/tests/test_proj.py index 5d95a8d..0108e62 100644 --- a/tests/test_proj.py +++ b/tests/test_proj.py @@ -3,19 +3,31 @@ import os os.environ["OGR_ORGANIZE_POLYGONS"] = "SKIP" from pathlib import Path +import geopandas as gpd # 添加项目根目录到路径 project_root = Path(__file__).parent.parent sys.path.insert(0, str(project_root)) -import geo_tools -from geo_tools.core import projection -from geo_tools.config.project_enum import CRS +import app +from app.core import projection +from app.config.project_enum import CRS -info = projection.get_crs_info(CRS.CGCS2000_6_DEGREE_ZONE_18.value) -print(info) -print(type(CRS.CGCS2000_3_DEGREE_ZONE_27)) +input_folder = r"E:\@三普\@二普和测土配方样点处理\云南省" +output_folder = r"E:\@三普\@二普和测土配方样点处理\云南省_投影转换" -# aa = geo_tools.read_vector(r"E:\@三普\@临时文件夹\样点异常值剔除\容县\异常样点数据\AB_outliers.shp") -# projection.reproject_gdf(aa,CRS.CGCS2000_3_DEGREE_ZONE_37).to_file(r"E:\@三普\@临时文件夹\样点异常值剔除\容县\AB_ou.shp") \ No newline at end of file +for root, dirs, files in os.walk(input_folder): + for file in files: + if file.endswith(".shp"): + input_path = os.path.join(root, file) + output_path = os.path.join(output_folder, os.path.relpath(input_path, input_folder)) + os.makedirs(os.path.dirname(output_path), exist_ok=True) + + gdf = gpd.read_file(input_path) + if gdf.crs is not None: + reproject_gdf = projection.reproject_gdf(gdf, auto_crs=True) + reproject_gdf.to_file(output_path) + print("Reprojected CRS:", reproject_gdf.crs) + else: + print("No CRS information found.") \ No newline at end of file diff --git a/tests/test_vector.py b/tests/test_vector.py index 7ec21d5..f55d39d 100644 --- a/tests/test_vector.py +++ b/tests/test_vector.py @@ -4,7 +4,7 @@ import pytest import geopandas as gpd from shapely.geometry import Point -from geo_tools.core.vector import ( +from app.core.vector import ( add_area_column, clip_to_extent, dissolve_by, diff --git a/项目分析报告.md b/项目分析报告.md new file mode 100644 index 0000000..bb2a13f --- /dev/null +++ b/项目分析报告.md @@ -0,0 +1,260 @@ +# Geo-Tools 项目分析报告 + +## 1. 项目概览 + +Geo-Tools 是一个专业的地理信息数据处理工具库，基于 geopandas、shapely、fiona 等开源地理空间库，提供了统一、简洁的 API 接口，简化地理空间数据处理流程。 + +- **核心价值**：为地理空间数据处理提供标准化、高效的工具集 +- **应用场景**：科学研究、GIS 分析、土地利用规划、农业分析等 +- **技术特点**：统一数据读写接口、丰富的空间分析功能、坐标系转换支持 + +## 2. 目录结构 + +项目采用模块化设计，清晰分离核心功能、IO 操作、分析工具和辅助功能。主要目录结构如下： + +``` +geo_tools/ +├── app/ # 主应用目录 +│ ├── analysis/ # 空间分析模块 +│ │ ├── spatial_ops.py # 空间叠加与邻域分析 +│ │ └── stats.py # 空间统计工具 +│ ├── config/ # 配置模块 +│ ├── core/ # 核心地理处理层 +│ │ ├── geometry.py # 几何运算工具 +│ │ ├── projection.py # 坐标系转换 +│ │ ├── raster.py # 栅格数据处理 +│ │ └── vector.py # 矢量数据处理 +│ ├── io/ # 数据读写模块 +│ │ ├── readers.py # 数据读取 +│ │ └── writers.py # 数据写入 +│ └── utils/ # 工具模块 +│ ├── config.py # 配置工具 +│ ├── logger.py # 日志工具 +│ └── validators.py # 数据验证 +├── data/ # 数据目录 +│ └── sample/ # 示例数据 +├── scripts/ # 脚本目录 +│ └── 其他工具/ # 特定领域工具脚本 +├── tests/ # 测试目录 +└── pyproject.toml # 项目配置文件 +``` + +## 3. 核心功能模块 + +### 3.1 数据读写（IO） + +**功能描述**：提供统一的矢量数据读取接口，支持多种格式的数据读写。 + +**支持格式**： +- Shapefile (.shp) +- GeoJSON (.geojson / .json) +- GeoPackage (.gpkg) +- File Geodatabase (.gdb) +- KML / KMZ +- FlatGeobuf (.fgb) +- CSV（含 WKT 或 经纬度列） + +**核心方法**： +- `read_vector()`：统一的矢量数据读取入口 +- `read_gdb()`：读取 Esri File Geodatabase +- `read_gpkg()`：读取 GeoPackage 文件 +- `read_csv_points()`：从 CSV 读取点数据 +- `list_gdb_layers()`：列出 GDB 中所有图层 + +### 3.2 几何处理（Core） + +**功能描述**：基于 Shapely 提供几何对象的创建、操作和分析功能。 + +**核心功能**： +- **几何有效性**：验证和修复几何对象 +- **基础几何运算**：缓冲区、质心、边界框、凸包 +- **集合运算**：交集、并集、差集、对称差集 +- **空间关系判断**：包含、在内部、相交、距离计算 + +**核心方法**： +- `is_valid_geometry()`：判断几何对象是否有效 +- `fix_geometry()`：尝试修复无效几何 +- `buffer_geometry()`：对几何对象执行缓冲区运算 +- `intersect()/union()/difference()`：几何集合运算 +- `contains()/within()/intersects()`：空间关系判断 + +### 3.3 坐标系转换（Core） + +**功能描述**：基于 pyproj 提供坐标系信息查询和投影转换功能。 + +**核心功能**： +- CRS 信息查询 +- 坐标点批量转换 +- 矢量数据重投影 +- 自动推荐适合的投影 CRS + +**核心方法**： +- `get_crs_info()`：获取 CRS 的基本信息 +- `transform_coordinates()`：批量转换坐标点 +- `reproject_gdf()`：将 GeoDataFrame 重投影到目标坐标系 +- `suggest_projected_crs()`：根据经纬度自动推荐投影 CRS + +### 3.4 空间分析（Analysis） + +**功能描述**：提供空间叠加、邻域分析和空间统计功能。 + +**核心功能**： +- **空间叠加**：缓冲区叠加、图层叠加 +- **邻域分析**：最近邻分析 +- **位置选择**：按空间关系选择要素 +- **空间统计**：面积加权均值、属性汇总、点计数 + +**核心方法**： +- `buffer_and_overlay()`：缓冲区后执行叠置分析 +- `nearest_features()`：找到最近的要素 +- `select_by_location()`：按位置关系选择要素 +- `area_weighted_mean()`：计算面积加权均值 +- `summarize_attributes()`：属性统计汇总 +- `count_by_polygon()`：统计面要素内的点数量 + +## 4. 技术架构与依赖 + +### 4.1 核心依赖 + +| 依赖库 | 版本要求 | 用途 | +|-------|---------|------| +| geopandas | >=0.14 | 地理空间数据处理核心库 | +| shapely | >=2.0 | 几何对象操作 | +| fiona | >=1.9 | 矢量数据读写 | +| pyproj | >=3.6 | 坐标系转换 | +| pandas | >=2.0 | 数据处理 | +| numpy | >=1.24 | 数值计算 | +| pydantic | >=2.0 | 数据验证 | +| pydantic-settings | >=2.0 | 配置管理 | +| python-dotenv | >=1.0 | 环境变量管理 | + +### 4.2 可选依赖 + +| 依赖组 | 包含库 | 用途 | +|-------|-------|------| +| dev | pytest, black, isort, flake8, mypy | 测试和代码质量工具 | +| notebook | jupyter, matplotlib, contextily, folium | 交互式分析和可视化 | +| raster | rasterio, xarray, rio-cogeo | 栅格数据处理 | + +### 4.3 技术架构 + +项目采用分层架构设计： + +1. **底层依赖层**：基于 geopandas、shapely、fiona 等核心库 +2. **核心处理层**：提供几何运算、坐标系转换等基础功能 +3. **IO 层**：统一数据读写接口，支持多种格式 +4. **分析层**：提供空间分析和统计功能 +5. **工具层**：提供配置、日志、验证等辅助功能 +6. **应用层**：通过脚本和 API 接口供用户使用 + +## 5. 业务定位与应用场景 + +### 5.1 业务定位 + +Geo-Tools 定位为专业的地理信息数据处理工具库，主要服务于以下场景： + +- **科学研究**：为地理空间相关研究提供数据处理和分析工具 +- **GIS 专业人员**：简化日常工作流程，提高处理效率 +- **土地利用规划**：支持空间分析和统计功能 +- **农业分析**：如耕作层厚度栅格制作、面积加权均值计算等 + +### 5.2 应用场景 + +从项目中的脚本和功能来看，主要应用场景包括： + +1. **数据格式转换**：在不同地理数据格式之间转换 +2. **空间分析**：缓冲区分析、叠加分析、邻域分析 +3. **空间统计**：面积加权均值、属性汇总、点计数 +4. **坐标系管理**：不同坐标系之间的转换和管理 +5. **农业分析**：如耕作层厚度分析、区域面积统计 + +### 5.3 典型工作流 + +1. **数据读取**：从各种格式读取地理空间数据 +2. **数据校验**：验证几何有效性和坐标系信息 +3. **数据处理**：几何运算、坐标系转换 +4. **空间分析**：执行叠加分析、邻域分析等 +5. **统计计算**：计算面积加权均值、属性汇总等 +6. **结果输出**：将处理结果写入各种格式 + +## 6. 核心 API/类/函数 + +### 6.1 数据读写 API + +| 函数名 | 功能描述 | 参数说明 | 返回值 | +|-------|---------|---------|-------| +| `read_vector()` | 统一的矢量数据读取入口 | path: 数据路径
layer: 图层名
crs: 目标坐标系 | GeoDataFrame | +| `read_gdb()` | 读取 Esri File Geodatabase | gdb_path: GDB 路径
layer: 图层名 | GeoDataFrame | +| `list_gdb_layers()` | 列出 GDB 中所有图层 | gdb_path: GDB 路径 | list[str] | +| `write_vector()` | 统一的矢量数据写入 | gdf: GeoDataFrame
path: 输出路径 | None | + +### 6.2 几何处理 API + +| 函数名 | 功能描述 | 参数说明 | 返回值 | +|-------|---------|---------|-------| +| `is_valid_geometry()` | 判断几何对象是否有效 | geom: 几何对象 | bool | +| `fix_geometry()` | 尝试修复无效几何 | geom: 几何对象 | BaseGeometry | +| `buffer_geometry()` | 对几何对象执行缓冲区运算 | geom: 几何对象
distance: 缓冲距离 | BaseGeometry | +| `intersect()/union()/difference()` | 几何集合运算 | geom_a, geom_b: 几何对象 | BaseGeometry | + +### 6.3 坐标系转换 API + +| 函数名 | 功能描述 | 参数说明 | 返回值 | +|-------|---------|---------|-------| +| `get_crs_info()` | 获取 CRS 的基本信息 | crs_input: CRS 输入 | dict | +| `transform_coordinates()` | 批量转换坐标点 | xs, ys: 坐标序列
source_crs: 源 CRS
target_crs: 目标 CRS | (list[float], list[float]) | +| `reproject_gdf()` | 将 GeoDataFrame 重投影 | gdf: GeoDataFrame
target_crs: 目标 CRS | GeoDataFrame | +| `suggest_projected_crs()` | 自动推荐投影 CRS | lon, lat: 经纬度 | str | + +### 6.4 空间分析 API + +| 函数名 | 功能描述 | 参数说明 | 返回值 | +|-------|---------|---------|-------| +| `buffer_and_overlay()` | 缓冲区后执行叠置分析 | source: 源图层
distance: 缓冲距离
target: 目标图层 | GeoDataFrame | +| `nearest_features()` | 找到最近的要素 | source: 查询图层
target: 被查询图层
k: 最近邻数量 | GeoDataFrame | +| `select_by_location()` | 按位置关系选择要素 | source: 源图层
selector: 选择图层
predicate: 空间谓词 | GeoDataFrame | +| `area_weighted_mean()` | 计算面积加权均值 | gdf: GeoDataFrame
value_col: 值列名
group_col: 分组列名 | pd.Series 或 pd.DataFrame | + +## 7. 技术亮点与优势 + +1. **统一的数据读写接口**：支持多种格式，简化数据输入输出流程 +2. **丰富的几何运算功能**：基于 Shapely 提供全面的几何操作 +3. **智能的坐标系管理**：自动推荐适合的投影 CRS，简化坐标系转换 +4. **强大的空间分析能力**：支持缓冲区分析、叠加分析、邻域分析等 +5. **灵活的空间统计功能**：面积加权均值、属性汇总、点计数等 +6. **良好的错误处理**：提供友好的错误信息，便于调试 +7. **模块化设计**：清晰的代码结构，便于维护和扩展 +8. **完善的文档**：详细的函数文档和示例代码 + +## 8. 潜在应用领域 + +1. **土地利用规划**：分析土地利用类型、面积统计 +2. **农业分析**：耕作层厚度分析、农田面积计算 +3. **环境监测**：生态区域分析、污染源影响范围评估 +4. **城市规划**：基础设施布局、服务范围分析 +5. **交通规划**：路线分析、站点服务范围 +6. **灾害评估**：灾害影响范围分析、风险评估 +7. **市场分析**：商业网点覆盖范围、人口分布分析 + +## 9. 项目发展建议 + +1. **增强栅格数据处理能力**：目前栅格处理是可选依赖，可考虑增强其功能 +2. **添加更多空间分析算法**：如空间插值、网络分析等 +3. **开发可视化工具**：集成更多可视化功能，便于结果展示 +4. **提供更多领域特定工具**：针对不同领域开发专用工具 +5. **完善测试覆盖**：增加更多测试用例，提高代码质量 +6. **优化性能**：对大型数据集的处理性能进行优化 +7. **提供更多示例**：增加不同领域的应用示例 + +## 10. 总结 + +Geo-Tools 是一个功能强大、设计合理的地理信息数据处理工具库，基于成熟的开源地理空间库，提供了统一、简洁的 API 接口，简化了地理空间数据的处理流程。 + +该项目的核心价值在于： +- 提供了标准化的地理空间数据处理工具 +- 支持多种数据格式的读写 +- 提供丰富的空间分析和统计功能 +- 简化了坐标系管理和转换 +- 为不同领域的地理空间分析提供了基础工具 + +Geo-Tools 不仅是一个实用的工具库，也是学习地理空间数据处理的良好资源，通过模块化的设计和清晰的文档，为用户提供了一个易于使用和扩展的地理信息处理框架。 \ No newline at end of file