933 lines
35 KiB
Python
933 lines
35 KiB
Python
# 生成完整的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()
|