404 lines
16 KiB
Python
404 lines
16 KiB
Python
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) |