missum/ArcGis_Py
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

初始化

Browse Source
This commit is contained in:
missum
2026-04-22 12:27:49 +08:00
commit 4857cb6e45
73 changed files with 20927 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

0
tools/core/acid_stats/__init__.py Normal file
View File

610
tools/core/acid_stats/土地利用类型酸化统计表.py Normal file
View File

@@ -0,0 +1,610 @@
# -*- coding: utf-8 -*-
import os
import arcpy
import pandas as pd
import numpy as np
from collections import OrderedDict
from openpyxl import Workbook
from openpyxl.styles import Font
from openpyxl.utils import get_column_letter
from tools.config.arcgis_field_cal_code import codeblock_cal_shfj, codeblock_dltb_ejdl, codeblock_dltb_yjdl
from tools.core.utils.excel_utils import ExcelStyleUtils
yjdl_order = ["耕地", "园地", "林地", "草地", "其他"]
ejdl_order = ["水田", "旱地", "水浇地", "果园", "茶园", "橡胶园", "其他园地"]
# --- 2. 辅助函数 ---
# 等级计算
def get_acidification_degree(delta_ph):
"""根据ΔpH值判断酸化程度"""
if pd.isna(delta_ph) or delta_ph == 0:
return "-"
# 请根据您的实际分级标准调整这里的阈值
if delta_ph > 1.0:
return "重度酸化"
elif 0.5 < delta_ph <= 1.0:
return "中度酸化"
elif 0.3 < delta_ph <= 0.5:
return "轻度酸化"
elif -0.3 <= delta_ph <= 0.3:
return "未酸化"
else: # dPH < -0.3
return "碱化"
# --- 3. 数据处理与分析 均值---
def process_data_for_table5_3(gdb_path, mean_table_name, sample_table_name):
"""
【最终版 v2】: 增加对制图样点数的处理，以支持加权平均计算。
"""
print("【最终版 v2】开始处理数据...")
def clean_df(df, columns):
# ... (此函数不变)
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 处理样点数据，计算“样点均值” ---
print("--> 步骤1: 计算样点均值...")
sample_table_path = os.path.join(gdb_path, sample_table_name)
sample_fields = ['YJDL', 'EJDL', 'dPH']
df_samples = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, sample_fields,'dPH>0.3', skip_nulls=False))
df_samples = clean_df(df_samples, ['YJDL', 'EJDL'])
# 按 YJDL, EJDL 分组，计算 dPH 的均值
df_sample_means = df_samples.groupby(['YJDL', 'EJDL'])['dPH'].mean().reset_index()
df_sample_means.rename(columns={'dPH': '样点均值'}, inplace=True)
print("样点均值计算完成。")
# --- b. 处理制图数据，获取“制图均值”和“制图样点数” ---
print("--> 步骤2: 获取制图均值和样点数...")
mean_table_path = os.path.join(gdb_path, mean_table_name)
# **【核心修改】: 增加读取 COUNT 字段**
mean_fields = ['YJDL', 'EJDL', 'MEAN', 'COUNT']
df_map_data = pd.DataFrame(arcpy.da.TableToNumPyArray(mean_table_path, mean_fields, skip_nulls=False))
df_map_data = clean_df(df_map_data, ['YJDL', 'EJDL'])
df_map_data.rename(columns={'MEAN': '制图均值', 'COUNT': '制图样点数'}, inplace=True)
print("制图数据获取完成。")
# --- c. 合并数据 ---
print("--> 步骤3: 合并数据...")
df_skeleton = pd.concat([
df_sample_means[['YJDL', 'EJDL']],
df_map_data[['YJDL', 'EJDL']]
]).drop_duplicates().reset_index(drop=True)
df_final = pd.merge(df_skeleton, df_sample_means, on=['YJDL', 'EJDL'], how='left')
# **【核心修改】: 合并整个 df_map_data，而不仅仅是均值列**
df_final = pd.merge(df_final, df_map_data, on=['YJDL', 'EJDL'], how='left')
# --- d. 计算酸化程度 ---
print("--> 步骤4: 计算酸化程度...")
# **【核心修改】: 在计算酸化程度之前，先过滤掉不展示的行**
# 我们只对 dPH 在酸化范围内 ( > 0.3) 的数据感兴趣
# 但为了计算合计，我们需要保留所有数据，所以这一步只计算，不删除
df_final['酸化程度_样本'] = df_final['样点均值'].apply(get_acidification_degree)
df_final['酸化程度_制图'] = df_final['制图均值'].apply(get_acidification_degree)
# (可选) 按“一级地类”和“二级地类”排序
in_ejdl_order = ejdl_order + [x for x in df_final['EJDL'].unique() if x not in ejdl_order]
df_final["YJDL"] = pd.Categorical(df_final['YJDL'], categories=yjdl_order, ordered=True)
df_final["EJDL"] = pd.Categorical(df_final['EJDL'], categories=in_ejdl_order, ordered=True)
df_final.sort_values(['YJDL', 'EJDL'], inplace=True)
print("数据处理流程完成！")
return df_final
# --- 4. Excel 制表 均值---
def write_to_excel_table5_3(df, output_path):
"""
将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel。")
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同土地利用类型pH变化统计"
# --- b. 绘制表头 ---
ws.merge_cells('A1:B1'); ws['A1'] = '土地利用类型'
ws.merge_cells('C1:F1'); ws['C1'] = 'ΔpH'
ws['A2'] = '一级'
ws['B2'] = '二级'
ws['C2'] = '样点均值'
ws['D2'] = '酸化程度'
ws['E2'] = '制图均值'
ws['F2'] = '酸化程度'
# --- c. 填充数据 ---
current_row = 3
# **【核心修改】: 先对整个DataFrame进行过滤，只保留需要展示的行**
# 只有当“样点酸化程度”或“制图酸化程度”不为“未酸化”、“碱化”或“-”时，才展示该行
acid_levels_to_show = ["轻度酸化", "中度酸化", "重度酸化"]
df_to_write = df[
df['酸化程度_样本'].isin(acid_levels_to_show) |
df['酸化程度_制图'].isin(acid_levels_to_show)
].copy() # 使用 .copy() 避免 SettingWithCopyWarning
for yl, group_yl_df in df_to_write.groupby('YJDL', sort=False, observed=False):
print(f"正在写入一级地类: {yl}...")
yl_start_row = current_row
# 遍历该一级地类下的所有“二级地类”
for _, row_data in group_yl_df.iterrows():
ws.cell(row=current_row, column=2).value = row_data['EJDL']
# 填充样点数据
sample_mean = row_data.get('样点均值')
if pd.notna(sample_mean):
ws.cell(row=current_row, column=3).value = f"{sample_mean:.2f}" if sample_mean > 0.3 else "-"
ws.cell(row=current_row, column=4).value = row_data.get('酸化程度_样本', '-') if sample_mean > 0.3 else "-"
else:
ws.cell(row=current_row, column=3).value = "-"
ws.cell(row=current_row, column=4).value = "-"
# 填充制图数据
map_mean = row_data.get('制图均值')
if pd.notna(map_mean):
ws.cell(row=current_row, column=5).value = f"{map_mean:.2f}" if map_mean > 0.3 else "-"
ws.cell(row=current_row, column=6).value = row_data.get('酸化程度_制图', '-') if map_mean > 0.3 else "-"
else:
ws.cell(row=current_row, column=5).value = "-"
ws.cell(row=current_row, column=6).value = "-"
current_row += 1
# 计算并写入“合计”行
if ws.cell(row=current_row-1, column=2).value in ["林地", "草地", "其他"]:
ws.merge_cells(start_row=yl_start_row, start_column=1, end_row=yl_start_row, end_column=2)
ws.cell(row=yl_start_row, column=1).value = yl
continue
ws.cell(row=current_row, column=2).value = '合计'
# 计算合计行的均值 (均值的均值)
total_sample_mean = group_yl_df['样点均值'].mean()
if pd.notna(total_sample_mean):
ws.cell(row=current_row, column=3).value = f"{total_sample_mean:.2f}"
ws.cell(row=current_row, column=4).value = get_acidification_degree(total_sample_mean)
else:
ws.cell(row=current_row, column=3).value = "-"
ws.cell(row=current_row, column=4).value = "-"
# b. **【核心修正】: 计算合计行的“制图均值”(加权平均)**
# 准备加权平均的分子和分母
weighted_sum = 0
total_count = 0
# 遍历当前一级地类分组中的每一行
for _, row in group_yl_df.iterrows():
mean_val = row.get('制图均值')
count_val = row.get('制图样点数')
# 只有当均值和样点数都存在且有效时，才参与计算
if pd.notna(mean_val) and pd.notna(count_val) and count_val > 0:
weighted_sum += mean_val * count_val # Σ (mean * count)
total_count += count_val # Σ (count)
# 计算加权平均值
weighted_avg = (weighted_sum / total_count) if total_count > 0 else 0
if weighted_avg > 0:
ws.cell(row=current_row, column=5).value = f"{weighted_avg:.2f}"
ws.cell(row=current_row, column=6).value = get_acidification_degree(weighted_avg)
else:
ws.cell(row=current_row, column=5).value = "-"
ws.cell(row=current_row, column=6).value = "-"
# 合并“一级地类”单元格
if yl_start_row <= current_row:
ws.merge_cells(start_row=yl_start_row, start_column=1, end_row=current_row, end_column=1)
ws.cell(row=yl_start_row, column=1).value = yl
current_row += 1
# --- a. 定义样式 ---
header_font = Font(name='等线', size=11, bold=True)
# --- d. 应用样式和调整列宽 ---
max_col_letter = get_column_letter(ws.max_column)
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}{current_row-1}')
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}2', header_font)
print("正在自动调整列宽...")
# 自动调整列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
# --- 2. 数据处理与分析 (使用 Pandas) ---
def process_data_for_table5_4(gdb_path, area_table_name, sample_table_name, target_area_dict):
"""
【最终修正版 v2】: 先建立统一的层级结构，再分别合并统计结果。
"""
print("【最终修正版 v2】开始处理数据...")
def clean_df(df, columns):
# ... (此函数不变)
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 从两个表中提取并建立唯一的 (YJDL, EJDL) 层级结构 "骨架" ---
print("--> 步骤1: 建立统一的层级结构...")
sample_table_path = os.path.join(gdb_path, sample_table_name)
area_table_path = os.path.join(gdb_path, area_table_name)
df_samples_raw = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, ['YJDL', 'EJDL'], skip_nulls=False))
df_area_raw = pd.DataFrame(arcpy.da.TableToNumPyArray(area_table_path, ['YJDL', 'EJDL'], skip_nulls=False))
# 清理并合并两个表中的 (YJDL, EJDL) 组合
df_samples_raw = clean_df(df_samples_raw, ['YJDL', 'EJDL'])
df_area_raw = clean_df(df_area_raw, ['YJDL', 'EJDL'])
# 使用 concat 连接两个DataFrame，然后用 drop_duplicates 去除重复的组合
df_skeleton = pd.concat([df_samples_raw, df_area_raw]).drop_duplicates().reset_index(drop=True)
if df_skeleton.empty:
print("警告: 无法从源数据中建立任何有效的 (YJDL, EJDL) 层级结构。")
return pd.DataFrame(), {}
print(f"已建立包含 {len(df_skeleton)} 个唯一土壤类型的层级结构。")
# --- b. 独立统计样点数据 ---
print("--> 步骤2: 独立统计样点数据...")
df_samples = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, ['EJDL', 'YJDL', 'dPH'], skip_nulls=False))
df_samples = clean_df(df_samples, ['YJDL', 'EJDL'])
if not df_samples.empty:
# ... (统计逻辑不变)
bins = [-np.inf, -0.3, 0.3, 0.5, 1.0, np.inf]
labels = ["碱化", "未酸化", "轻度酸化", "中度酸化", "重度酸化"]
df_samples['SHFJ'] = pd.cut(df_samples['dPH'], bins=bins, labels=labels, right=True)
sample_counts = df_samples.groupby(['YJDL', 'EJDL', 'SHFJ'], observed=False).size().reset_index(name='样点数')
ts_total_samples = sample_counts.groupby(['YJDL', 'EJDL'])['样点数'].transform('sum')
sample_counts['样点占比'] = (sample_counts['样点数'] / ts_total_samples) * 100
df_sample_stats = sample_counts.pivot_table(
index=['YJDL', 'EJDL'], columns='SHFJ', values=['样点数', '样点占比'], fill_value=0, observed=False
).reset_index()
df_sample_stats.columns = [f'{col[0]}_{col[1]}'.strip('_') if col[1] else col[0] for col in df_sample_stats.columns]
# 将样点统计结果合并到骨架上
df_final = pd.merge(df_skeleton, df_sample_stats, on=['YJDL', 'EJDL'], how='left')
else:
df_final = df_skeleton.copy()
# --- c. 独立统计面积数据 ---
print("--> 步骤3: 独立统计面积数据...")
df_area = pd.DataFrame(arcpy.da.TableToNumPyArray(area_table_path, ['EJDL', 'YJDL', 'SHFJ', 'AREA'], skip_nulls=False))
df_area = clean_df(df_area, ['YJDL', 'EJDL'])
if not df_area.empty:
# 计算平差系数
landuse_types = {'耕地':'01', '园地':'02', '林地':'03', '草地':'04', '其他':'12'}
df_area['AREA_MU'] = df_area['AREA'] * 0.0015
yjdl_area = df_area.groupby(['YJDL'])['AREA_MU'].sum().reset_index()
yjdl_area.columns = ['YJDL', 'ORIGINAL_TOTAL_MU']
adjustment_factors = []
for _, row in yjdl_area.iterrows():
yjdl = row['YJDL']
original_total = row['ORIGINAL_TOTAL_MU']
target_total = target_area_dict.get(landuse_types[yjdl], original_total) # 如果没有指定，就用原始面积
adjustment_factor = target_total / original_total
adjustment_factors.append({
'YJDL': yjdl,
'原始总面积_亩': original_total,
'目标总面积_亩': target_total,
'平差系数': adjustment_factor
})
factor_df = pd.DataFrame(adjustment_factors)
# 4. 对每个二级地类应用平差系数
# 合并原始数据和平差系数
df_with_factors = df_area.merge(factor_df[['YJDL', '平差系数']], on='YJDL')
df_with_factors['制图面积_亩'] = df_with_factors['AREA_MU'] * df_with_factors['平差系数']
ts_total_area = df_with_factors.groupby(['YJDL', 'EJDL'])['制图面积_亩'].transform('sum')
df_with_factors['面积占比'] = (df_with_factors['制图面积_亩'] / ts_total_area) * 100
df_area_stats = df_with_factors.pivot_table(
index=['YJDL', 'EJDL'], columns='SHFJ', values=['制图面积_亩', '面积占比'], fill_value=0
).reset_index()
df_area_stats.columns = [f'{col[0]}_{col[1]}'.strip('_') if col[1] else col[0] for col in df_area_stats.columns]
# 将面积统计结果合并到 df_final 上
# 注意，这里我们合并到已经包含样点数据的 df_final 上
df_final = pd.merge(df_final, df_area_stats, on=['YJDL', 'EJDL'], how='left')
# --- d. 最后清理和构建映射 ---
df_final.fillna(0, inplace=True)
print("--> 步骤4: 自动构建层级结构...")
dynamic_soil_mapping = df_final.groupby('YJDL')['EJDL'].unique().apply(list).to_dict()
dynamic_soil_mapping = OrderedDict(sorted(dynamic_soil_mapping.items(),key=lambda item: yjdl_order.index(item[0])))
in_ejdl_order = ejdl_order + [x for x in df_final['EJDL'].unique() if x not in ejdl_order]
for yl in dynamic_soil_mapping:
# dynamic_soil_mapping[yl].sort()
dynamic_soil_mapping[yl] = sorted( dynamic_soil_mapping[yl], key=lambda x: in_ejdl_order.index(x))
print("数据处理流程完成！")
return df_final, dynamic_soil_mapping
# --- 3. Excel 制表 面积---
def write_to_excel_table5_4(df, soil_mapping, output_path):
"""
【最终修正版】: 将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel，将创建一个空的报告。")
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同类型土壤酸化程度统计"
# --- b. 绘制表头 (不变) ---
ws.merge_cells('A1:B1'); ws['A1'] = '土地利用类型'
ws['A2'] = '一级'
ws['B2'] = '二级'
acid_levels = ['轻度酸化', '中度酸化', '重度酸化']
all_possible_levels = ['碱化', '未酸化', '轻度酸化', '中度酸化', '重度酸化']
acid_level_headers = ['轻度酸化(0.3<ΔpH≤0.5)', '中度酸化(0.5<ΔpH≤1.0)', '重度酸化(ΔpH>1.0)']
col_start = 3
for header in acid_level_headers:
ws.merge_cells(start_row=1, start_column=col_start, end_row=1, end_column=col_start + 3)
ws.cell(row=1, column=col_start).value = header
ws.cell(row=2, column=col_start).value = '样点数/个'
ws.cell(row=2, column=col_start + 1).value = '占比/%'
ws.cell(row=2, column=col_start + 2).value = '制图面积/亩'
ws.cell(row=2, column=col_start + 3).value = '占比/%'
col_start += 4
# --- c. 填充数据 (完全重构的逻辑) ---
current_row = 3
# 使用 .groupby('YJDL', sort=False) 来保证我们之前设置的排序顺序
for yl, ts_list in soil_mapping.items():
# **【关键】** group_yl 是一个只包含当前一级地类数据的子DataFrame
# 我们可以安全地在这个子DataFrame上进行迭代和计算
print(f"正在写入一级地类: {yl}...")
yl_start_row = current_row
# 筛选出当前一级地类的所有数据
group_yl_df = df[df['YJDL'] == yl]
# 1. 遍历该一级地类下的所有“二级地类”并写入数据
for ts in ts_list:
ws.cell(row=current_row, column=2).value = ts
# 在子集中查找当前二级地类的数据行
row_data = group_yl_df[group_yl_df['EJDL'] == ts]
# --- 填充单元格的逻辑开始 ---
col_start = 3 # 从第 C 列开始填充
# 检查是否找到了该土属的数据
if not row_data.empty:
# 如果找到了数据 (row_data 不为空)，我们就获取这一行的数据
# .iloc[0] 获取第一行（也是唯一一行）的数据，作为一个 Series 对象
data_series = row_data.iloc[0]
# 遍历每一个酸化等级，填充对应的四列数据
for level in acid_levels:
# 1. 构建要从 data_series 中查找的列名
sample_col = f'样点数_{level}'
sample_pct_col = f'样点占比_{level}'
area_col = f'制图面积_亩_{level}'
area_pct_col = f'面积占比_{level}'
# 2. 从 data_series 中安全地获取值
# 使用 .get(key, default_value) 的好处是，如果列名不存在，它会返回默认值(0)，而不会报错
sample_val = data_series.get(sample_col, 0)
sample_pct_val = data_series.get(sample_pct_col, 0)
area_val = data_series.get(area_col, 0)
area_pct_val = data_series.get(area_pct_col, 0)
# 3. 将获取到的值填入单元格
# - 对于数值，我们判断它是否大于0。如果是，就填入数值；否则，填入 "-"
# - 对于样点数，我们将其转为整数
# - 对于占比和面积，我们保留两位小数
# 样点数/个
ws.cell(row=current_row, column=col_start).value = int(sample_val) if sample_val > 0 else "-"
# 占比/%
ws.cell(row=current_row, column=col_start + 1).value = f"{sample_pct_val:.2f}%" if sample_val > 0 else "-"
# 制图面积/万亩
ws.cell(row=current_row, column=col_start + 2).value = f"{area_val:.0f}" if area_val > 0 else "-"
# 占比/%
ws.cell(row=current_row, column=col_start + 3).value = f"{area_pct_val:.2f}%" if area_val > 0 else "-"
# 移动到下一个酸化等级的起始列
col_start += 4
else:
# 如果没有找到该土属的数据 (row_data 为空)
# 这意味着该土属在源数据中不存在任何样点或面积信息
# 我们将整行所有统计单元格都填充为 "-"
# acid_levels 列表包含3个等级，每个等级4列，总共12列
for _ in range(len(acid_levels) * 4):
ws.cell(row=current_row, column=col_start).value = "-"
col_start += 1
# --- 填充单元格的逻辑结束 ---
# 完成一行填充后，行号加1，为下一行做准备
current_row += 1
# 2. 计算并写入这个一级地类的“合计”行
if ws.cell(row=current_row-1, column=2).value in ["林地","草地", "其他"]:
ws.merge_cells(start_row=yl_start_row, start_column=1, end_row=yl_start_row, end_column=2)
ws.cell(row=yl_start_row, column=1).value = yl
continue
ws.cell(row=current_row, column=2).value = '合计'
# 计算总样点数和总面积，仅针对当前 group_yl
yl_grand_total_samples = 0
for lvl in all_possible_levels:
if f'样点数_{lvl}' in group_yl_df:
yl_grand_total_samples += group_yl_df[f'样点数_{lvl}'].sum()
yl_grand_total_area = 0
for lvl in all_possible_levels:
if f'制图面积_亩_{lvl}' in group_yl_df:
yl_grand_total_area += group_yl_df[f'制图面积_亩_{lvl}'].sum()
col_start = 3
for level in acid_levels:
sample_sum = group_yl_df.get(f'样点数_{level}', 0).sum()
col_name = f'制图面积_亩_{level}'
area_sum = group_yl_df[col_name].sum() if col_name in group_yl_df else 0
# area_sum = group_yl_df.get(f'制图面积_亩_{level}', 0).sum()
sample_perc = (sample_sum / yl_grand_total_samples * 100) if yl_grand_total_samples > 0 else 0
area_perc = (area_sum / yl_grand_total_area * 100) if yl_grand_total_area > 0 else 0
ws.cell(row=current_row, column=col_start).value = int(sample_sum) if sample_sum > 0 else "-"
ws.cell(row=current_row, column=col_start + 1).value = f"{sample_perc:.2f}%" if sample_sum > 0 else "-"
ws.cell(row=current_row, column=col_start + 2).value = f"{area_sum:.0f}" if area_sum > 0 else "-"
ws.cell(row=current_row, column=col_start + 3).value = f"{area_perc:.2f}%" if area_sum > 0 else "-"
col_start += 4
# 3. 合并“一级地类”单元格
if yl_start_row <= current_row:
ws.merge_cells(start_row=yl_start_row, start_column=1, end_row=current_row, end_column=1)
ws.cell(row=yl_start_row, column=1).value = yl
current_row += 1
# --- a. 定义样式 (不变) ---
header_font = Font(name='等线', size=11, bold=True)
# --- d. 应用样式和调整列宽 (最终健壮版) ---
max_col_letter = get_column_letter(ws.max_column)
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}{current_row-1}')
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}2', header_font)
print("正在自动调整列宽...")
# 调整列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
def main(gdb_path:str, ph_features:str,dltb_class_feature:str, shph_tif:str, output_path:str,target_areas_dict:dict):
try:
# --- 1. 用户配置 ---
# 输出配置
output_excel_path = os.path.join(output_path, "土地利用类型酸化统计表.xlsx") # 生成的Excel报告文件路径
# 设置工作空间和变量
arcpy.env.workspace = gdb_path
arcpy.env.overwriteOutput = True
sample_table_name = "历史样点PH信息_Table" # 图2: 样点信息表名
in_zone_feature = dltb_class_feature # 地类图斑
in_class_feature = ph_features # 已重分类好的酸化PH图层
in_value_raster = shph_tif # 赋值栅格,酸化PH栅格
out_table_area = r"土地利用类型_酸化面积表" # 输出的面积统计表名
out_table_mean = r"土地利用类型_酸化均值表" # 输出的均值表名
print("开始处理数据...")
if not arcpy.Exists(out_table_area):
# 判断输入表是否存在SHFJ字段
try:
if not arcpy.ListFields(in_zone_feature, "EJDL"):
arcpy.management.CalculateField(in_zone_feature, "EJDL", "calculate_ejdl(!DLBM!,!DLMC!)", "PYTHON3", codeblock_dltb_ejdl)
arcpy.management.CalculateField(in_zone_feature, "YJDL", "calculate_yjdl(!DLBM!)", "PYTHON3", codeblock_dltb_yjdl)
if not arcpy.ListFields(in_class_feature, "SHFJ"):
arcpy.management.CalculateField(in_class_feature, "SHFJ", "calculate_shfj(!gridcode!)", "PYTHON3", codeblock_cal_shfj)
except Exception as e:
print(f"计算SHFJ字段时发生错误: {e}")
# 拿到地类图斑的坐标系
desc = arcpy.Describe(in_zone_feature)
spatial_ref = desc.spatialReference
# 1.用arcpy.analysis.TabulateIntersection进行交集制表,面积使用地类图斑投影坐标系下面积
with arcpy.EnvManager(outputCoordinateSystem=spatial_ref):
arcpy.analysis.TabulateIntersection(
in_zone_feature,
["YJDL", "EJDL"],
in_class_feature,
out_table_area,
"SHFJ",
out_units="SQUARE_METERS",
)
if not arcpy.Exists(out_table_mean):
# 判断输入表是否存在YJDL_EJDL字段
if not arcpy.ListFields(in_zone_feature, "YJDL_EJDL"):
# 如果不存在，则添加该字段
arcpy.management.AddField(in_zone_feature, "YJDL_EJDL", "TEXT")
# 计算YJDL_EJDL字段的值
arcpy.management.CalculateField(in_zone_feature,"YJDL_EJDL","!YJDL! + '_' + !EJDL!","PYTHON3")
# 2.用arcpy.sa.ZonalStatisticsAsTable进行区域统计
mean_table = arcpy.sa.ZonalStatisticsAsTable(
in_zone_feature, "YJDL_EJDL", in_value_raster, out_table_mean, "DATA", "MEAN"
)
# 2.1 添加土壤类型字段并计算
arcpy.management.AddFields(
out_table_mean,
[["YJDL", "TEXT"],["EJDL", "TEXT"]],
)
arcpy.management.CalculateField(mean_table, "YJDL", "!YJDL_EJDL!.split('_')[0]", "PYTHON3")
arcpy.management.CalculateField(mean_table, "EJDL", "!YJDL_EJDL!.split('_')[1]", "PYTHON3")
# 生成表5.4的面积统计Excel报告
final_dataframe, soil_structure = process_data_for_table5_4(gdb_path, out_table_area, sample_table_name,target_areas_dict)
write_to_excel_table5_4(final_dataframe, soil_structure, output_excel_path)
# 生成表5.3的均值统计Excel报告
final_mean_dataframe = process_data_for_table5_3(gdb_path, out_table_mean, sample_table_name)
write_to_excel_table5_3(final_mean_dataframe, output_excel_path.replace(".xlsx", "_mean.xlsx"))
except Exception as e:
print(f"\n处理过程中发生严重错误: {e}")
import traceback
traceback.print_exc()
finally:
import gc
gc.collect()
# --- 4. 主程序入口 ---
# if __name__ == "__main__":
# main()

629
tools/core/acid_stats/土壤类型图酸化统计表.py Normal file
View File

@@ -0,0 +1,629 @@
# -*- coding: utf-8 -*-
import os
import arcpy
import pandas as pd
import numpy as np
from openpyxl import Workbook
from openpyxl.styles import Font
from openpyxl.utils import get_column_letter
from tools.config.arcgis_field_cal_code import codeblock_cal_shfj
from tools.core.utils.excel_utils import ExcelStyleUtils
from tools.config.custom_sort import yl_order, ts_order
# --- 2. 辅助函数 ---
# 获取要素类各酸化等级面积
def get_acid_area_by_group(target_area_df):
try:
# 转为numpy数组供pandas统计使用
df = target_area_df.copy()
area_by_group = df.groupby("SHFJ")["AREA_MU"].sum()
for key in area_by_group.keys():
area_by_group[key] = area_by_group[key]
return area_by_group.to_dict()
except Exception as e:
print(f"计算面积时出错: {str(e)}")
return None
def apply_adjustment_by_each_level(df, target_area_dict):
"""
对DataFrame中的面积数据按每一个酸化等级独立进行平差。
参数:
df (pd.DataFrame): 包含面积统计的DataFrame。
target_area_dict (dict): 每个酸化等级的目标总面积字典。
例如: {'轻度酸化': 10000.0, '中度酸化': 8000.0, ...}
"""
print("\n开始按每个酸化等级独立进行平差...")
df_adjusted = df.copy()
for level, target_area in target_area_dict.items():
col_name = f'制图面积_亩_{level}'
adjusted_col_name = f'平差后面积_亩_{level}'
if col_name not in df.columns:
print(f"警告: 未找到列 '{col_name}'，跳过该等级平差。")
if adjusted_col_name not in df_adjusted.columns:
df_adjusted[adjusted_col_name] = 0 # 创建一个空列
continue
# a. 计算该等级的实际总面积
actual_area = df_adjusted[col_name].sum()
if actual_area > 0:
# b. 计算误差
error = target_area - actual_area
print(f"等级 '{level}': 目标面积={target_area:.2f}, 实际面积={actual_area:.2f}, 误差={error:.2f}")
# c. 按比例分配误差
adjustment = error * (df_adjusted[col_name] / actual_area)
df_adjusted[adjusted_col_name] = df_adjusted[col_name] + adjustment
df_adjusted[adjusted_col_name] = df_adjusted[adjusted_col_name].clip(lower=0)
else:
df_adjusted[adjusted_col_name] = df_adjusted[col_name]
print("按每个酸化等级独立平差完成。")
return df_adjusted
# 获取酸化程度
def get_acidification_degree(delta_ph):
"""根据ΔpH值判断酸化程度"""
if pd.isna(delta_ph) or delta_ph == 0:
return "-"
# 请根据您的实际分级标准调整这里的阈值
if delta_ph > 1.0:
return "重度酸化"
elif 0.5 < delta_ph <= 1.0:
return "中度酸化"
elif 0.3 < delta_ph <= 0.5:
return "轻度酸化"
elif -0.3 <= delta_ph <= 0.3:
return "未酸化"
else: # dPH < -0.3
return "碱化"
# --- 3. 数据处理与分析 均值表---
def process_data_for_table5_5(gdb_path, mean_table_name, sample_table_name):
"""
【最终版 v2】: 增加对制图样点数的处理，以支持加权平均计算。
"""
print("【最终版 v2】开始处理数据...")
def clean_df(df, columns):
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 处理样点数据，计算“样点均值” ---
print("--> 步骤1: 计算样点均值...")
sample_table_path = os.path.join(gdb_path, sample_table_name)
sample_fields = ['YL', 'TS', 'dPH']
df_samples = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, sample_fields, 'dPH>0.3', skip_nulls=False))
df_samples = clean_df(df_samples, ['YL', 'TS'])
# 按 YL, TS 分组，计算 dPH 的均值
df_sample_means = df_samples.groupby(['YL', 'TS'])['dPH'].mean().reset_index()
df_sample_means.rename(columns={'dPH': '样点均值'}, inplace=True)
print("样点均值计算完成。")
# --- b. 处理制图数据，获取“制图均值”和“制图样点数” ---
print("--> 步骤2: 获取制图均值和样点数...")
mean_table_path = os.path.join(gdb_path, mean_table_name)
mean_fields = ['YL', 'TS', 'MEAN', 'COUNT']
df_map_data = pd.DataFrame(arcpy.da.TableToNumPyArray(mean_table_path, mean_fields, skip_nulls=False))
df_map_data = clean_df(df_map_data, ['YL', 'TS'])
df_map_data.rename(columns={'MEAN': '制图均值', 'COUNT': '制图样点数'}, inplace=True)
print("制图数据获取完成。")
# --- c. 合并数据 ---
print("--> 步骤3: 合并数据...")
df_skeleton = pd.concat([
df_sample_means[['YL', 'TS']],
df_map_data[['YL', 'TS']]
]).drop_duplicates().reset_index(drop=True)
df_final = pd.merge(df_skeleton, df_sample_means, on=['YL', 'TS'], how='left')
# **【核心修改】: 合并整个 df_map_data，而不仅仅是均值列**
df_final = pd.merge(df_final, df_map_data, on=['YL', 'TS'], how='left')
# --- d. 计算酸化程度 ---
print("--> 步骤4: 计算酸化程度...")
# **【核心修改】: 在计算酸化程度之前，先过滤掉不展示的行**
# 我们只对 dPH 在酸化范围内 ( > 0.3) 的数据感兴趣
# 但为了计算合计，我们需要保留所有数据，所以这一步只计算，不删除
df_final['酸化程度_样本'] = df_final['样点均值'].apply(get_acidification_degree)
df_final['酸化程度_制图'] = df_final['制图均值'].apply(get_acidification_degree)
# (可选) 按“亚类”和“土属”排序
in_yl_order = yl_order + [x for x in df_final['YL'].unique() if x not in yl_order]
in_ts_order = ts_order + [x for x in df_final['TS'].unique() if x not in ts_order]
df_final["YL"] = pd.Categorical(df_final['YL'], categories=in_yl_order, ordered=True)
df_final["TS"] = pd.Categorical(df_final['TS'], categories=in_ts_order, ordered=True)
df_final.sort_values(['YL', 'TS'], inplace=True)
print("数据处理流程完成！")
return df_final
# --- 4. Excel 制表 均值表---
def write_to_excel_table5_5(df, output_path):
"""
将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel。")
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同类型土壤pH变化统计"
# --- b. 绘制表头 ---
ws.merge_cells('A1:A2'); ws['A1'] = '亚类'
ws.merge_cells('B1:B2'); ws['B1'] = '土属'
ws.merge_cells('C1:F1'); ws['C1'] = 'ΔpH'
ws['C2'] = '样点均值'
ws['D2'] = '酸化程度'
ws['E2'] = '制图均值'
ws['F2'] = '酸化程度'
# --- c. 填充数据 ---
current_row = 3
# **【核心修改】: 先对整个DataFrame进行过滤，只保留需要展示的行**
# 只有当“样点酸化程度”或“制图酸化程度”不为“未酸化”、“碱化”或“-”时，才展示该行
acid_levels_to_show = ["轻度酸化", "中度酸化", "重度酸化"]
df_to_write = df[
df['酸化程度_样本'].isin(acid_levels_to_show) |
df['酸化程度_制图'].isin(acid_levels_to_show)
].copy() # 使用 .copy() 避免 SettingWithCopyWarning
for yl, group_yl_df in df_to_write.groupby('YL', observed=True, sort=False):
print(f"正在写入亚类: {yl}...")
yl_start_row = current_row
# 遍历该亚类下的所有“土属”
for _, row_data in group_yl_df.iterrows():
ws.cell(row=current_row, column=2).value = row_data['TS']
# 填充样点数据
sample_mean = row_data.get('样点均值')
if pd.notna(sample_mean):
ws.cell(row=current_row, column=3).value = f"{sample_mean:.2f}" if sample_mean > 0.3 else "-"
ws.cell(row=current_row, column=4).value = row_data.get('酸化程度_样本', '-') if sample_mean > 0.3 else "-"
else:
ws.cell(row=current_row, column=3).value = "-"
ws.cell(row=current_row, column=4).value = "-"
# 填充制图数据
map_mean = row_data.get('制图均值')
if pd.notna(map_mean):
ws.cell(row=current_row, column=5).value = f"{map_mean:.2f}" if map_mean > 0.3 else "-"
ws.cell(row=current_row, column=6).value = row_data.get('酸化程度_制图', '-') if map_mean > 0.3 else "-"
else:
ws.cell(row=current_row, column=5).value = "-"
ws.cell(row=current_row, column=6).value = "-"
current_row += 1
# 计算并写入“合计”行
ws.cell(row=current_row, column=2).value = '合计'
# 计算合计行的均值 (均值的均值)
total_sample_mean = group_yl_df['样点均值'].mean()
if pd.notna(total_sample_mean):
ws.cell(row=current_row, column=3).value = f"{total_sample_mean:.2f}"
ws.cell(row=current_row, column=4).value = get_acidification_degree(total_sample_mean)
else:
ws.cell(row=current_row, column=3).value = "-"
ws.cell(row=current_row, column=4).value = "-"
# b. **【核心修正】: 计算合计行的“制图均值”(加权平均)**
# 准备加权平均的分子和分母
weighted_sum = 0
total_count = 0
# 遍历当前亚类分组中的每一行
for _, row in group_yl_df.iterrows():
mean_val = row.get('制图均值')
count_val = row.get('制图样点数')
# 只有当均值和样点数都存在且有效时，才参与计算
if pd.notna(mean_val) and pd.notna(count_val) and count_val > 0:
weighted_sum += mean_val * count_val # Σ (mean * count)
total_count += count_val # Σ (count)
# 计算加权平均值
weighted_avg = (weighted_sum / total_count) if total_count > 0 else 0
if weighted_avg > 0:
ws.cell(row=current_row, column=5).value = f"{weighted_avg:.2f}"
ws.cell(row=current_row, column=6).value = get_acidification_degree(weighted_avg)
else:
ws.cell(row=current_row, column=5).value = "-"
ws.cell(row=current_row, column=6).value = "-"
# 合并“亚类”单元格
if yl_start_row <= current_row:
ws.merge_cells(start_row=yl_start_row, start_column=1, end_row=current_row, end_column=1)
ws.cell(row=yl_start_row, column=1).value = yl
current_row += 1
# --- a. 定义样式 ---
header_font = Font(name='等线', size=11, bold=True)
# --- d. 应用样式和调整列宽 (最终健壮版) ---
max_col_letter = get_column_letter(ws.max_column)
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}{current_row-1}')
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}2', header_font)
print("正在自动调整列宽...")
# 自动调整列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
# --- 2. 数据处理与分析 (面积统计表) ---
def process_data_final(gdb_path, area_table_name, sample_table_name):
"""
【最终修正版 v2】: 先建立统一的层级结构，再分别合并统计结果。
"""
print("【最终修正版 v2】开始处理数据...")
def clean_df(df, columns):
# ... (此函数不变)
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 从两个表中提取并建立唯一的 (YL, TS) 层级结构 "骨架" ---
print("--> 步骤1: 建立统一的层级结构...")
sample_table_path = os.path.join(gdb_path, sample_table_name)
area_table_path = os.path.join(gdb_path, area_table_name)
df_samples_raw = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, ['YL', 'TS'], skip_nulls=False))
df_area_raw = pd.DataFrame(arcpy.da.TableToNumPyArray(area_table_path, ['YL', 'TS'], skip_nulls=False))
# 清理并合并两个表中的 (YL, TS) 组合
df_samples_raw = clean_df(df_samples_raw, ['YL', 'TS'])
df_area_raw = clean_df(df_area_raw, ['YL', 'TS'])
# 使用 concat 连接两个DataFrame，然后用 drop_duplicates 去除重复的组合
df_skeleton = pd.concat([df_samples_raw, df_area_raw]).drop_duplicates().reset_index(drop=True)
if df_skeleton.empty:
print("警告: 无法从源数据中建立任何有效的 (YL, TS) 层级结构。")
return pd.DataFrame(), {}
print(f"已建立包含 {len(df_skeleton)} 个唯一土壤类型的层级结构。")
# --- b. 独立统计样点数据 ---
print("--> 步骤2: 独立统计样点数据...")
df_samples = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, ['TS', 'YL', 'dPH'], skip_nulls=False))
df_samples = clean_df(df_samples, ['YL', 'TS'])
if not df_samples.empty:
bins = [-np.inf, -0.3, 0.3, 0.5, 1.0, np.inf]
labels = ["碱化", "未酸化", "轻度酸化", "中度酸化", "重度酸化"]
df_samples['SHFJ'] = pd.cut(df_samples['dPH'], bins=bins, labels=labels, right=True)
sample_counts = df_samples.groupby(['YL', 'TS', 'SHFJ'], observed=False).size().reset_index(name='样点数')
ts_total_samples = sample_counts.groupby(['YL', 'TS'])['样点数'].transform('sum')
sample_counts['样点占比'] = (sample_counts['样点数'] / ts_total_samples) * 100
df_sample_stats = sample_counts.pivot_table(
index=['YL', 'TS'], columns='SHFJ', values=['样点数', '样点占比'], fill_value=0, observed=False
).reset_index()
df_sample_stats.columns = [f'{col[0]}_{col[1]}'.strip('_') if col[1] else col[0] for col in df_sample_stats.columns]
# 将样点统计结果合并到骨架上
df_final = pd.merge(df_skeleton, df_sample_stats, on=['YL', 'TS'], how='left')
else:
df_final = df_skeleton.copy()
# --- c. 独立统计面积数据 ---
print("--> 步骤3: 独立统计面积数据...")
df_area = pd.DataFrame(arcpy.da.TableToNumPyArray(area_table_path, ['TS', 'YL', 'SHFJ', 'AREA'], skip_nulls=False))
df_area = clean_df(df_area, ['YL', 'TS'])
if not df_area.empty:
df_area['制图面积_亩'] = df_area['AREA'] * 0.0015
ts_total_area = df_area.groupby(['YL', 'TS'])['制图面积_亩'].transform('sum')
df_area['面积占比'] = (df_area['制图面积_亩'] / ts_total_area) * 100
df_area_stats = df_area.pivot_table(
index=['YL', 'TS'], columns='SHFJ', values=['制图面积_亩', '面积占比'], fill_value=0
).reset_index()
df_area_stats.columns = [f'{col[0]}_{col[1]}'.strip('_') if col[1] else col[0] for col in df_area_stats.columns]
# 将面积统计结果合并到 df_final 上
# 注意，这里我们合并到已经包含样点数据的 df_final 上
df_final = pd.merge(df_final, df_area_stats, on=['YL', 'TS'], how='left')
# --- d. 最后清理和构建映射 ---
df_final.fillna(0, inplace=True)
print("--> 步骤4: 自动构建层级结构...")
in_yl_order = yl_order + [x for x in df_final['YL'].unique() if x not in yl_order]
in_ts_order = ts_order + [x for x in df_final['TS'].unique() if x not in ts_order]
df_final["YL"] = pd.Categorical(df_final['YL'], categories=in_yl_order, ordered=True)
df_final["TS"] = pd.Categorical(df_final['TS'], categories=in_ts_order, ordered=True)
df_final.sort_values(['YL', 'TS'], inplace=True)
dynamic_soil_mapping = df_final.groupby('YL', observed=True)['TS'].unique().apply(list).to_dict()
# for yl in dynamic_soil_mapping:
# dynamic_soil_mapping[yl].sort()
print("数据处理流程完成！")
return df_final, dynamic_soil_mapping
# --- 3. Excel 制表 面积统计表 ---
def write_to_excel(df, soil_mapping, output_path):
"""
【最终修正版】: 将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel，将创建一个空的报告。")
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同类型土壤酸化程度统计"
# --- b. 绘制表头 (不变) ---
ws.merge_cells('A1:A2'); ws['A1'] = '亚类'
ws.merge_cells('B1:B2'); ws['B1'] = '土属'
acid_levels = ['轻度酸化', '中度酸化', '重度酸化']
all_possible_levels = ['碱化', '未酸化', '轻度酸化', '中度酸化', '重度酸化']
acid_level_headers = ['轻度酸化(0.3<ΔpH≤0.5)', '中度酸化(0.5<ΔpH≤1.0)', '重度酸化(ΔpH>1.0)']
col_start = 3
for header in acid_level_headers:
ws.merge_cells(start_row=1, start_column=col_start, end_row=1, end_column=col_start + 3)
ws.cell(row=1, column=col_start).value = header
ws.cell(row=2, column=col_start).value = '样点数/个'
ws.cell(row=2, column=col_start + 1).value = '占比/%'
ws.cell(row=2, column=col_start + 2).value = '制图面积/亩'
ws.cell(row=2, column=col_start + 3).value = '占比/%'
col_start += 4
# --- c. 填充数据 (完全重构的逻辑) ---
current_row = 3
# 使用 .groupby('YL', sort=False) 来保证我们之前设置的排序顺序
for yl, ts_list in soil_mapping.items():
# **【关键】** group_yl 是一个只包含当前亚类数据的子DataFrame
# 我们可以安全地在这个子DataFrame上进行迭代和计算
print(f"正在写入亚类: {yl}...")
yl_start_row = current_row
# 筛选出当前亚类的所有数据
group_yl_df = df[df['YL'] == yl]
# 1. 遍历该亚类下的所有“土属”并写入数据
for ts in ts_list:
ws.cell(row=current_row, column=2).value = ts
# 在子集中查找当前土属的数据行
row_data = group_yl_df[group_yl_df['TS'] == ts]
# --- 填充单元格的逻辑开始 ---
col_start = 3 # 从第 C 列开始填充
# 检查是否找到了该土属的数据
if not row_data.empty:
# 如果找到了数据 (row_data 不为空)，我们就获取这一行的数据
# .iloc[0] 获取第一行（也是唯一一行）的数据，作为一个 Series 对象
data_series = row_data.iloc[0]
# 遍历每一个酸化等级，填充对应的四列数据
for level in acid_levels:
# 1. 构建要从 data_series 中查找的列名
sample_col = f'样点数_{level}'
sample_pct_col = f'样点占比_{level}'
area_col = f'平差后面积_亩_{level}'
area_pct_col = f'面积占比_{level}'
# 2. 从 data_series 中安全地获取值
# 使用 .get(key, default_value) 的好处是，如果列名不存在，它会返回默认值(0)，而不会报错
sample_val = data_series.get(sample_col, 0)
sample_pct_val = data_series.get(sample_pct_col, 0)
area_val = data_series.get(area_col, 0)
area_pct_val = data_series.get(area_pct_col, 0)
# 3. 将获取到的值填入单元格
# - 对于数值，我们判断它是否大于0。如果是，就填入数值；否则，填入 "-"
# - 对于样点数，我们将其转为整数
# - 对于占比和面积，我们保留两位小数
# 样点数/个
ws.cell(row=current_row, column=col_start).value = int(sample_val) if sample_val > 0 else "-"
# 占比/%
ws.cell(row=current_row, column=col_start + 1).value = f"{sample_pct_val:.2f}%" if sample_val > 0 else "-"
# 制图面积/亩
ws.cell(row=current_row, column=col_start + 2).number_format = "0.00"
ws.cell(row=current_row, column=col_start + 2).value = f"{area_val:.0f}" if area_val > 0 else "-"
# 占比/%
ws.cell(row=current_row, column=col_start + 3).value = f"{area_pct_val:.2f}%" if area_val > 0 else "-"
# 移动到下一个酸化等级的起始列
col_start += 4
else:
# 如果没有找到该土属的数据 (row_data 为空)
# 这意味着该土属在源数据中不存在任何样点或面积信息
# 我们将整行所有统计单元格都填充为 "-"
# acid_levels 列表包含3个等级，每个等级4列，总共12列
for _ in range(len(acid_levels) * 4):
ws.cell(row=current_row, column=col_start).value = "-"
col_start += 1
# --- 填充单元格的逻辑结束 ---
# 完成一行填充后，行号加1，为下一行做准备
current_row += 1
# 2. 计算并写入这个亚类的“合计”行
ws.cell(row=current_row, column=2).value = '合计'
# 计算总样点数和总面积，仅针对当前 group_yl
yl_grand_total_samples = 0
for lvl in all_possible_levels:
if f'样点数_{lvl}' in group_yl_df:
yl_grand_total_samples += group_yl_df[f'样点数_{lvl}'].sum()
yl_grand_total_area = 0
for lvl in all_possible_levels:
if f'制图面积_亩_{lvl}' in group_yl_df:
yl_grand_total_area += group_yl_df[f'制图面积_亩_{lvl}'].sum()
col_start = 3
for level in acid_levels:
sample_sum = group_yl_df.get(f'样点数_{level}', 0).sum()
col_name = f'制图面积_亩_{level}'
area_sum = group_yl_df[col_name].sum() if col_name in group_yl_df else 0
# area_sum = group_yl_df.get(f'平差后面积_亩_{level}', 0).sum()
sample_perc = (sample_sum / yl_grand_total_samples * 100) if yl_grand_total_samples > 0 else 0
area_perc = (area_sum / yl_grand_total_area * 100) if yl_grand_total_area > 0 else 0
ws.cell(row=current_row, column=col_start).value = int(sample_sum) if sample_sum > 0 else "-"
ws.cell(row=current_row, column=col_start + 1).value = f"{sample_perc:.2f}%" if sample_sum > 0 else "-"
ws.cell(row=current_row, column=col_start + 2).value = f"{area_sum:.0f}" if area_sum > 0 else "-"
ws.cell(row=current_row, column=col_start + 3).value = f"{area_perc:.2f}%" if area_sum > 0 else "-"
col_start += 4
# 3. 合并“亚类”单元格
if yl_start_row <= current_row:
ws.merge_cells(start_row=yl_start_row, start_column=1, end_row=current_row, end_column=1)
ws.cell(row=yl_start_row, column=1).value = yl
current_row += 1
# --- a. 定义样式 (不变) ---
header_font = Font(name='等线', size=11, bold=True)
# d. 应用样式和调整列宽
max_col = 2 + len(acid_levels) * 4
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:{get_column_letter(max_col)}{current_row-1}')
ExcelStyleUtils.set_style(ws, f'A1:{get_column_letter(max_col)}2', header_font)
# 调整列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
def main(gdb_path, trlx_polygon, sh_ph_polygon, ph_raster, output_path, target_areas_df):
try:
# --- 1. 用户配置 ---
sample_table_name = "历史样点PH信息_Table" # 图2: 样点信息表名
# 输出配置
output_excel_path = os.path.join(output_path, "土壤类型酸化统计表.xlsx") # 生成的Excel报告文件路径
# 设置工作空间和变量
arcpy.env.workspace = gdb_path
arcpy.env.overwriteOutput = True
in_zone_feature = trlx_polygon # 土壤类型图
# in_class_feature = sh_ph_polygon # 已重分类好的酸化PH图层
in_class_feature = "最小面积统计单元"
in_value_raster = ph_raster # 酸化PH栅格
dltb_ph_statstable = "土地利用类型_酸化面积表" # 土壤类型_酸化面积表（gdb table）
out_table_area = r"土壤类型_酸化面积表" # 输出的交集表名
out_table_mean = r"土壤类型_酸化均值表" # 输出的均值表名
print("开始处理数据...")
if not arcpy.Exists(out_table_area):
# 判断输入表是否存在SHFJ字段
try:
arcpy.management.CalculateField(in_class_feature, "SHFJ", "calculate_shfj(!gridcode!)", "PYTHON3", codeblock_cal_shfj)
except Exception as e:
print(f"计算SHFJ字段时发生错误: {e}")
# 1.用arcpy.analysis.TabulateIntersection进行交集制表
arcpy.analysis.TabulateIntersection(
in_zone_feature,
["TS", "YL"],
in_class_feature,
out_table_area,
"SHFJ",
out_units="SQUARE_METERS",
)
if not arcpy.Exists(out_table_mean):
# 判断输入表是否存在YL_TS字段
if not arcpy.ListFields(in_zone_feature, "YL_TS"):
# 如果不存在，则添加该字段
arcpy.management.AddField(in_zone_feature, "YL_TS", "TEXT")
# 计算YL_TS字段的值
arcpy.management.CalculateField(in_zone_feature,"YL_TS","!YL! + '_' + !TS!","PYTHON3")
# 2.用arcpy.sa.ZonalStatisticsAsTable进行区域统计
mean_table = arcpy.sa.ZonalStatisticsAsTable(
in_zone_feature, "YL_TS", in_value_raster, out_table_mean, "DATA", "MEAN"
)
# 2.1 添加土壤类型字段并计算
arcpy.management.AddFields(
out_table_mean,
[["YL", "TEXT"],["TS", "TEXT"]],
)
arcpy.management.CalculateField(mean_table, "YL", "!YL_TS!.split('_')[0]", "PYTHON3")
arcpy.management.CalculateField(mean_table, "TS", "!YL_TS!.split('_')[1]", "PYTHON3")
# 生成表5.4的面积统计Excel报告
final_dataframe, soil_structure = process_data_final(gdb_path, out_table_area, sample_table_name)
# 统计地类图斑酸化总面积亩
each_acid_area = get_acid_area_by_group(target_areas_df)
print(f"容县土壤类型图斑总 acid 总面积（亩）：{each_acid_area}")
# 执行平差计算
if each_acid_area:
adjusted_dataframe = apply_adjustment_by_each_level(final_dataframe, each_acid_area)
print("使用平差值进行修正！")
write_to_excel(adjusted_dataframe, soil_structure, output_excel_path)
else:
print("未使用平差值进行修正！")
write_to_excel(final_dataframe, soil_structure, output_excel_path)
# 生成表5.4的均值统计Excel报告
final_mean_dataframe = process_data_for_table5_5(gdb_path, out_table_mean, sample_table_name)
write_to_excel_table5_5(final_mean_dataframe, output_excel_path.replace(".xlsx", "_mean.xlsx"))
# adjusted_dataframe.to_csv(output_excel_path.replace(".xlsx", "_adjusted.csv"), index=False)
except Exception as e:
print(f"\n处理过程中发生严重错误: {e}")
import traceback
traceback.print_exc()
finally:
import gc
gc.collect()
# --- 4. 主程序入口 ---
# if __name__ == "__main__":
# main()

167
tools/core/acid_stats/空间连接.py Normal file
View File

@@ -0,0 +1,167 @@
# -*- coding: utf-8 -*-
import sys
import arcpy
from pathlib import Path
sys.path.append(str(Path(__file__).parent))
from tools.config.arcgis_field_cal_code import codeblock_dltb_ejdl, codeblock_dltb_yjdl
def export_to_points(ph_points, dltb_features, trlx_features, xzq_features, assign_raster, workspace):
# --- 1. 设置工作空间和变量 ---
# 请根据您的实际情况修改以下路径
arcpy.env.workspace = workspace
arcpy.env.overwriteOutput = True
# 输入的要素类
input_features = ph_points # 历史样点PH数据
join_features_list = [trlx_features,xzq_features,dltb_features] # 连接图层 (规划分区)
# 输出的要素类
final_output_fc = "历史样点PH信息_Table"
# --- 3. 主处理逻辑 ---
try:
print("开始处理赋值样点PH信息...")
target_features = f"in_memory/temp_sample_raster"
# 将栅格数据提取至历史PH样点
arcpy.sa.ExtractValuesToPoints(
in_point_features=input_features,
in_raster=assign_raster,
out_point_features=target_features,
interpolate_values="NONE",
add_attributes="VALUE_ONLY"
)
print("开始计算地类一二级类别...")
# 计算地类图斑一级、二级类别
try:
arcpy.management.CalculateField(dltb_features, "EJDL", "calculate_ejdl(!DLBM!,!DLMC!)", "PYTHON3", codeblock_dltb_ejdl)
arcpy.management.CalculateField(dltb_features, "YJDL", "calculate_yjdl(!DLBM!)", "PYTHON3", codeblock_dltb_yjdl)
arcpy.management.CalculateField(dltb_features, "YJDLBM", "!DLBM![:2]", "PYTHON3")
raster_path = Path(assign_raster)
# if "二普" in raster_path.stem or "测土" in raster_path.stem:
arcpy.management.CalculateField(target_features, "dPH", "!RASTERVALU!-!PH!", "PYTHON3", field_type="DOUBLE")
# else:
# arcpy.management.CalculateField(target_features, "dPH", "!PH!-!RASTERVALU!", "PYTHON3", field_type="DOUBLE")
except Exception as e:
print(e)
# --- 2. 定义要保留的字段 ---
# 这是一个非常清晰的配置方式：指定每个图层要保留的字段列表
fields_to_keep = {
target_features: ["PH", "RASTERVALU", "dPH"],
trlx_features: ["YL", "TS"],
xzq_features: ["XZQMC"],
dltb_features: ["YJDL", "EJDL"]
}
print("开始配置字段映射...")
# 初始化当前的目标图层，最开始是原始的目标图层
current_target = target_features
# 存储所有中间生成的临时文件，以便最后清理
temp_outputs = []
temp_outputs.append(target_features)
# 获取目标图层的所有字段，以便在后续迭代中保留
retained_fields = fields_to_keep.get(target_features, [])
# 迭代处理每一个连接图层
for i, join_features in enumerate(join_features_list):
print(f"\n--- 开始处理第 {i+1} 个连接图层: {join_features} ---")
# 检查连接图层是否存在
if not arcpy.Exists(join_features):
print(f"警告: 连接图层 '{join_features}' 不存在，将跳过此连接。")
continue
# --- 配置 FieldMappings ---
field_mappings = arcpy.FieldMappings()
# a. 保留已经存在于 current_target 中的字段
# 这些字段是在之前的迭代中保留下来的
for field_name in retained_fields:
try:
field_map = arcpy.FieldMap()
field_map.addInputField(current_target, field_name)
field_mappings.addFieldMap(field_map)
except Exception:
# 如果字段在之前的某个步骤中未能成功添加，这里会捕获异常
print(f"注意: 在图层 '{current_target}' 中未找到字段 '{field_name}'，可能在之前的步骤中已被跳过。")
# b. 从当前的 join_features 中添加新字段
fields_from_current_join = fields_to_keep.get(join_features, [])
for field_name in fields_from_current_join:
try:
field_map = arcpy.FieldMap()
field_map.addInputField(join_features, field_name)
field_map.mergeRule = "First" # 对所有连接字段使用 "First" 规则
field_mappings.addFieldMap(field_map)
except Exception as e:
print(f"警告: 添加字段 '{field_name}' (来自 '{join_features}') 时出错，将跳过。错误信息: {e}")
# 如果本次迭代没有有效的字段映射，则跳过
if field_mappings.fieldCount == 0:
print(f"警告: 对于连接图层 '{join_features}' 没有有效的字段可以添加，跳过此连接。")
continue
# 定义本次连接的临时输出名
# 使用 in_memory 工作空间可以提高性能
temp_output = f"in_memory/temp_join_{i}"
temp_outputs.append(temp_output)
print(f"执行空间连接: '{current_target}' + '{join_features}' -> '{temp_output}'")
# 执行空间连接
arcpy.analysis.SpatialJoin(
target_features=current_target,
join_features=join_features,
out_feature_class=temp_output,
join_operation="JOIN_ONE_TO_ONE",
join_type="KEEP_ALL",
field_mapping=field_mappings,
match_option="INTERSECT"
)
# 更新 current_target 为本次操作的输出，以便下一次迭代使用
current_target = temp_output
# 更新已保留字段列表，为下一次迭代做准备
retained_fields.extend(fields_from_current_join)
print(f"连接成功。目前已保留的字段: {retained_fields}")
# --- 4. 保存最终结果并清理 ---
# 将最后一个临时输出复制或重命名为最终结果
if arcpy.Exists(current_target):
print(f"\n所有连接完成。将最终结果 '{current_target}' 保存为 '{final_output_fc}'...")
# arcpy.management.CopyFeatures(current_target, final_output_fc)
arcpy.conversion.ExportTable(current_target, final_output_fc)
print("最终结果已保存。")
# 验证输出字段
output_fields = [f.name for f in arcpy.ListFields(final_output_fc)]
print(f"最终输出的字段为: {output_fields}")
else:
print("警告: 没有任何连接操作成功执行，未生成最终输出。")
except arcpy.ExecuteError:
print("\n--- ArcPy 执行错误 ---")
print(arcpy.GetMessages(2))
except Exception as e:
print(f"\n--- 发生未预料的错误 ---")
print(e)
finally:
# 清理所有中间生成的临时文件
print("\n开始清理临时文件...")
for temp_file in temp_outputs:
if arcpy.Exists(temp_file):
arcpy.management.Delete(temp_file)
print(f"已删除临时文件: {temp_file}")
print("清理完成。")

641
tools/core/acid_stats/行政区划酸化统计表.py Normal file
View File

@@ -0,0 +1,641 @@
# -*- coding: utf-8 -*-
import os
import arcpy
import pandas as pd
import numpy as np
from openpyxl import Workbook
from openpyxl.styles import Font
from openpyxl.utils import get_column_letter
from tools.config.arcgis_field_cal_code import codeblock_cal_shfj
from tools.core.utils.excel_utils import ExcelStyleUtils
yjdl_order = ["耕地", "园地", "林地", "草地", "其他"]
ejdl_order = ["水田", "旱地", "水浇地", "果园", "茶园", "橡胶园", "其他园地"]
# --- 2. 辅助函数 ---
# 等级计算
def get_acidification_degree(delta_ph):
"""根据ΔpH值判断酸化程度"""
if pd.isna(delta_ph) or delta_ph == 0:
return "-"
# 请根据您的实际分级标准调整这里的阈值
if delta_ph > 1.0:
return "重度酸化"
elif 0.5 < delta_ph <= 1.0:
return "中度酸化"
elif 0.3 < delta_ph <= 0.5:
return "轻度酸化"
elif 0.1 < delta_ph <= 0.3:
return "弱酸化"
else: # dPH < -0.3
return "其他"
# --- 3. 数据处理与分析 均值---
def process_data_for_table5_7(gdb_path, mean_table_name, sample_table_name):
"""
【最终版 v2】: 增加对制图样点数的处理，以支持加权平均计算。
"""
print("开始处理数据...")
def clean_df(df, columns):
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 处理样点数据，计算“样点均值” ---
print("--> 步骤1: 计算样点均值...")
sample_table_path = os.path.join(gdb_path, sample_table_name)
sample_fields = ['XZQMC','YJDL','EJDL', 'dPH']
df_samples = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, sample_fields, 'dPH>0.3', skip_nulls=False))
df_samples = clean_df(df_samples, ['XZQMC','YJDL', 'EJDL'])
# 按 YJDL, EJDL 分组，计算 dPH 的均值
df_sample_means = df_samples.groupby(['XZQMC'])['dPH'].mean().reset_index()
df_sample_means.rename(columns={'dPH': '样点均值'}, inplace=True)
print("样点均值计算完成。")
# --- b. 处理制图数据，获取“制图均值”和“制图样点数” ---
print("--> 步骤2: 获取制图均值和样点数...")
mean_table_path = os.path.join(gdb_path, mean_table_name)
mean_fields = ['XZQMC', 'MEAN', 'COUNT']
df_map_data = pd.DataFrame(arcpy.da.TableToNumPyArray(mean_table_path, mean_fields, skip_nulls=False))
df_map_data = clean_df(df_map_data, ['XZQMC'])
df_map_data.rename(columns={'MEAN': '制图均值', 'COUNT': '制图样点数'}, inplace=True)
print("制图数据获取完成。")
# --- c. 合并数据 ---
print("--> 步骤3: 合并数据...")
df_skeleton = pd.concat([
df_sample_means[['XZQMC']],
df_map_data[['XZQMC']]
]).drop_duplicates().reset_index(drop=True)
df_final = pd.merge(df_skeleton, df_sample_means, on=['XZQMC'], how='left')
# **【核心修改】: 合并整个 df_map_data，而不仅仅是均值列**
df_final = pd.merge(df_final, df_map_data, on=['XZQMC'], how='left')
# --- d. 计算酸化程度 ---
print("--> 步骤4: 计算酸化程度...")
# **【核心修改】: 在计算酸化程度之前，先过滤掉不展示的行**
# 我们只对 dPH 在酸化范围内 ( > 0.3) 的数据感兴趣
# 但为了计算合计，我们需要保留所有数据，所以这一步只计算，不删除
df_final['酸化程度_样本'] = df_final['样点均值'].apply(get_acidification_degree)
df_final['酸化程度_制图'] = df_final['制图均值'].apply(get_acidification_degree)
df_final.sort_values(['XZQMC'], inplace=True)
print("数据处理流程完成！")
return df_final
# --- 4. Excel 制表 均值---
def write_to_excel_table5_7(df, output_path):
"""
将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel。")
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同土地利用类型pH变化统计"
# --- b. 绘制表头 ---
ws.merge_cells('A1:A2'); ws['A1'] = '乡镇/街道'
ws.merge_cells('B1:E1'); ws['B1'] = 'ΔpH'
ws['B2'] = '样点均值'
ws['C2'] = '酸化程度'
ws['D2'] = '制图均值'
ws['E2'] = '酸化程度'
# --- c. 填充数据 ---
current_row = 3
# **【核心修改】: 先对整个DataFrame进行过滤，只保留需要展示的行**
acid_levels_to_show = ["弱酸化", "轻度酸化", "中度酸化", "重度酸化", "其他"]
df_to_write = df[
df['酸化程度_样本'].isin(acid_levels_to_show) | df['酸化程度_制图'].isin(acid_levels_to_show)
].copy() # 使用 .copy() 避免 SettingWithCopyWarning
for _, row_data in df_to_write.iterrows():
print(f"正在写入一级地类...")
# 写入数据”
ws.cell(row=current_row, column=1).value = row_data['XZQMC']
# 填充样点数据
sample_mean = row_data.get('样点均值')
if pd.notna(sample_mean):
ws.cell(row=current_row, column=2).value = f"{sample_mean:.2f}" if sample_mean > 0.3 else "-"
ws.cell(row=current_row, column=3).value = row_data.get('酸化程度_样本', '-') if sample_mean > 0.3 else "-"
else:
ws.cell(row=current_row, column=2).value = "-"
ws.cell(row=current_row, column=3).value = "-"
# 填充制图数据
map_mean = row_data.get('制图均值')
if pd.notna(map_mean):
ws.cell(row=current_row, column=4).value = f"{map_mean:.2f}" if map_mean > 0.3 else "-"
ws.cell(row=current_row, column=5).value = row_data.get('酸化程度_制图', '-') if map_mean > 0.3 else "-"
else:
ws.cell(row=current_row, column=4).value = "-"
ws.cell(row=current_row, column=5).value = "-"
current_row += 1
# --- a. 定义样式 ---
header_font = Font(name='等线', size=11, bold=True)
# --- d. 应用样式和调整列宽 ---
max_col_letter = get_column_letter(ws.max_column)
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}{current_row-1}')
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}2', header_font)
print("正在自动调整列宽...")
# 设置列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
# --- 2. 数据处理与分析 面积 各乡镇---
def process_data_for_table5_4(gdb_path, area_table_name, target_area_dict):
"""
【最终修正版 v2】: 先建立统一的层级结构，再分别合并统计结果。
"""
print("【最终修正版 v2】开始处理数据...")
def clean_df(df, columns):
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 独立统计面积数据 ---
print("--> 步骤1: 独立统计面积数据...")
area_table_path = os.path.join(gdb_path, area_table_name)
df_area = pd.DataFrame(arcpy.da.TableToNumPyArray(area_table_path, ['XZQMC', 'SHFJ', 'AREA'], skip_nulls=False))
df_area = clean_df(df_area, ['XZQMC'])
df_final = pd.DataFrame()
if not df_area.empty:
# 计算平差系数
target_shfj_areas = target_area_dict.groupby(['SHFJ'])['AREA_MU'].sum().reset_index()
original_shfj_areas = df_area.groupby(['SHFJ'])['AREA'].sum().reset_index()
original_shfj_areas['AREA_MU'] = original_shfj_areas['AREA'] * 0.0015
adjustment_factors = []
for index, row in original_shfj_areas.iterrows():
shfj = row['SHFJ']
area_mu = row['AREA_MU']
adjustment_factor = target_shfj_areas[target_shfj_areas['SHFJ'] == shfj]['AREA_MU'].values[0] / area_mu
adjustment_factors.append({
'SHFJ': shfj,
'平差系数':adjustment_factor
})
factor_df = pd.DataFrame(adjustment_factors)
df_sh_area = df_area.merge(factor_df[['SHFJ', '平差系数']], on='SHFJ')
df_sh_area['制图面积_亩'] = df_sh_area['AREA'] * 0.0015 * df_sh_area['平差系数']
ts_total_area = df_sh_area.groupby(['XZQMC'])['制图面积_亩'].transform('sum')
df_sh_area['面积占比'] = (df_sh_area['制图面积_亩'] / ts_total_area) * 100
df_area_stats = df_sh_area.pivot_table(
index=['XZQMC'], columns='SHFJ', values=['制图面积_亩', '面积占比'], fill_value=0
).reset_index()
df_area_stats.columns = [f'{col[0]}_{col[1]}'.strip('_') if col[1] else col[0] for col in df_area_stats.columns]
df_final = df_area_stats
print("--> 步骤2: 计算酸化面积合计...")
# 定义属于酸化类别的面积列
acidic_area_cols = [
'制图面积_亩_轻度酸化',
'制图面积_亩_中度酸化',
'制图面积_亩_重度酸化'
]
# 确保这些列存在于DataFrame中，不存在的列用0代替
for col in acidic_area_cols:
if col not in df_final.columns:
df_final[col] = 0
# 将这三列相加，得到合计值
df_final['酸化面积合计_亩'] = df_final[acidic_area_cols].sum(axis=1)
# --- d. 最后清理和构建映射 ---
df_final.fillna(0, inplace=True)
print("数据处理流程完成！")
return df_final
# --- 3. Excel 制表 面积---
def write_to_excel_table5_4(df, output_path):
"""
【最终修正版】: 将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel，将创建一个空的报告。")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同乡镇酸化面积统计"
ws['A1'] = "没有有效的统计数据。"
wb.save(output_path)
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "不同乡镇酸化面积统计"
# --- b. 绘制表头 (不变) ---
ws.merge_cells('A1:A2'); ws['A1'] = '乡镇/街道'
acid_levels = ['弱酸化', '轻度酸化', '中度酸化', '重度酸化', '其他']
# acid_level_headers = ['0.1<ΔpH≤0.3', '0.3<ΔpH≤0.5', '0.5<ΔpH≤1.0', 'ΔpH>1.0', '其他']
# all_possible_levels = ['碱化', '未酸化', '轻度酸化', '中度酸化', '重度酸化']
acid_level_headers = ['弱酸化(0.1<ΔpH≤0.3)','轻度酸化(0.3<ΔpH≤0.5)', '中度酸化(0.5<ΔpH≤1.0)', '重度酸化(ΔpH>1.0)', '其他(未酸化)']
col_start = 2
for header in acid_level_headers:
ws.merge_cells(start_row=1, start_column=col_start, end_row=1, end_column=col_start + 1)
ws.cell(row=1, column=col_start).value = header
ws.cell(row=2, column=col_start).value = '面积/亩'
ws.cell(row=2, column=col_start + 1).value = '占比/%'
col_start += 2
# 增加合计列的表头**
total_col = col_start # 记录合计列的列号
ws.merge_cells(start_row=1, start_column=total_col, end_row=2, end_column=total_col)
ws.cell(row=1, column=total_col).value = '酸化面积合计'
# --- c. 填充数据 (完全重构的逻辑) ---
current_row = 3
# **【核心修改】: 不再需要 group_yl_df，直接遍历整个 df**
# 假设 df 已经按 XZQMC 排序（如果需要的话）
df_sorted = df.sort_values('XZQMC').reset_index(drop=True)
for index, row_data in df_sorted.iterrows():
ws.cell(row=current_row, column=1).value = row_data['XZQMC']
col_start = 2
for level in acid_levels:
area_col = f'制图面积_亩_{level}'
area_pct_col = f'面积占比_{level}'
area_val = row_data.get(area_col, 0)
area_pct_val = row_data.get(area_pct_col, 0)
ws.cell(row=current_row, column=col_start).value = f"{area_val:.0f}" if area_val > 0 else "-"
ws.cell(row=current_row, column=col_start + 1).value = f"{area_pct_val:.2f}%" if area_val > 0 else "-"
col_start += 2
# **【核心修改】: 填充酸化面积合计列的值**
total_area_val = row_data.get('酸化面积合计_亩', 0)
ws.cell(row=current_row, column=total_col).value = f"{total_area_val:.0f}" if total_area_val > 0 else "-"
current_row += 1
# **(可选) 增加一个所有乡镇的“总合计”行**
# print("--> 计算并写入总合计行...")
# ws.cell(row=current_row, column=1).value = '总合计'
# col_start = 2
# for level in acid_levels:
# area_col = f'制图面积_亩_{level}'
# area_sum = df_sorted.get(area_col, 0).sum()
# # 总合计行的占比是相对于所有乡镇的总面积
# grand_total_area = df_sorted[[f'制图面积_亩_{lvl}' for lvl in all_possible_levels if f'制图面积_亩_{lvl}' in df_sorted]].sum().sum()
# area_perc = (area_sum / grand_total_area * 100) if grand_total_area > 0 else 0
# ws.cell(row=current_row, column=col_start).value = f"{area_sum:.2f}" if area_sum > 0 else "-"
# ws.cell(row=current_row, column=col_start + 1).value = f"{area_perc:.2f}" if area_sum > 0 else "-"
# col_start += 2
# grand_total_acidic_area = df_sorted['酸化面积合计_亩'].sum()
# ws.cell(row=current_row, column=total_col).value = f"{grand_total_acidic_area:.2f}" if grand_total_acidic_area > 0 else "-"
# current_row += 1
# --- a. 定义样式 (不变) ---
header_font = Font(name='等线', size=11, bold=True)
# --- d. 应用样式和调整列宽 (最终健壮版) ---
max_col_letter = get_column_letter(ws.max_column)
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}{current_row-1}')
ExcelStyleUtils.set_style(ws, f'A1:{max_col_letter}2', header_font)
print("正在自动调整列宽...")
# 设置列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
# 步骤5.3: 生成表5.3 - 总表数据处理
def process_data_for_table5_2(gdb_path, area_table_name, sample_table_name, target_area_dict:pd.DataFrame):
def clean_df(df, columns):
for col in columns:
df[col] = df[col].astype(str).str.strip()
df.replace(['<Null>', 'None', '', '<空>'], np.nan, inplace=True)
df.dropna(subset=columns, inplace=True)
return df
# --- a. 从两个表中提取并建立唯一的 (YJDL, EJDL) 层级结构 "骨架" ---
print("--> 步骤1: 建立统一的层级结构...")
sample_table_path = os.path.join(gdb_path, sample_table_name)
area_table_path = os.path.join(gdb_path, area_table_name)
# --- b. 独立统计样点数据 ---
print("--> 步骤2: 独立统计样点数据...")
df_samples = pd.DataFrame(arcpy.da.TableToNumPyArray(sample_table_path, ['XZQMC', 'dPH'], skip_nulls=False))
df_samples = clean_df(df_samples, ['XZQMC'])
if not df_samples.empty:
bins = [-np.inf, 0.1, 0.3, 0.5, 1.0, np.inf]
labels = ["其他", "弱酸化", "轻度酸化", "中度酸化", "重度酸化"]
df_samples['SHFJ'] = pd.cut(df_samples['dPH'], bins=bins, labels=labels, right=True)
sample_counts = df_samples.groupby(['SHFJ'], observed=False).size().reset_index(name='样点数')
sample_counts = sample_counts.merge(df_samples.groupby(['SHFJ'], observed=False)['dPH'].mean(), on='SHFJ')
ts_total_samples = sample_counts['样点数'].sum()
sample_counts['样点占比'] = (sample_counts['样点数'] / ts_total_samples) * 100
# print(sample_counts)
# --- c. 独立统计面积数据 ---
print("--> 步骤3: 独立统计面积数据...")
df_area = pd.DataFrame(arcpy.da.TableToNumPyArray(area_table_path, ['XZQMC', 'SHFJ', 'AREA'], skip_nulls=False))
df_area = clean_df(df_area, ['XZQMC'])
if not df_area.empty:
# 计算平差系数
target_shfj_areas = target_area_dict.groupby(['SHFJ'])['AREA_MU'].sum().reset_index()
original_shfj_areas = df_area.groupby(['SHFJ'])['AREA'].sum().reset_index()
original_shfj_areas['AREA_MU'] = original_shfj_areas['AREA'] * 0.0015
adjustment_factors = []
for index, row in original_shfj_areas.iterrows():
shfj = row['SHFJ']
area_mu = row['AREA_MU']
adjustment_factor = target_shfj_areas[target_shfj_areas['SHFJ'] == shfj]['AREA_MU'].values[0] / area_mu
adjustment_factors.append({
'SHFJ': shfj,
'平差系数':adjustment_factor
})
factor_df = pd.DataFrame(adjustment_factors)
df_sh_area = df_area.merge(factor_df[['SHFJ', '平差系数']], on='SHFJ')
df_sh_area['制图面积_亩'] = df_sh_area['AREA'] * 0.0015 * df_sh_area['平差系数']
df_area_counts = df_sh_area.groupby(['SHFJ'], observed=False)[['制图面积_亩']].sum()
ts_total_area = df_area_counts['制图面积_亩'].sum()
df_area_counts['面积占比'] = (df_area_counts['制图面积_亩'] / ts_total_area) * 100
df_final = pd.merge(sample_counts, df_area_counts, on=['SHFJ'], how='left')
# # --- d. 最后清理和构建映射 ---
df_final.fillna(0, inplace=True)
return df_final
# --- 3. Excel 制表 总表---
def write_to_excel_table5_2(df, df_mean, output_path):
"""
【最终修正版】: 将处理好的数据写入格式化的 Excel 文件。
"""
if df.empty:
print("警告: 没有数据可以写入 Excel，将创建一个空的报告。")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws['A1'] = "没有有效的统计数据。"
wb.save(output_path)
return
print(f"开始生成 Excel 报告到 '{output_path}'...")
wb = Workbook()
ws = wb.create_sheet("Mysheet", 0)
ws.title = "行政区酸化程度等级分布及占比"
# --- b. 绘制表头 (不变) ---
ws.merge_cells('A1:B1'); ws['A1'] = '酸化程度'
ws.merge_cells('C1:D1'); ws['C1'] = '样点统计'
ws.merge_cells('E1:F1'); ws['E1'] = '制图统计'
ws.merge_cells('A8:B8'); ws['A8'] = '总计'
ws.merge_cells('A9:B9'); ws['A9'] = '全县酸化样点ΔpH 均值'
ws.merge_cells('A10:B10'); ws['A10'] = '全县酸化制图ΔpH 均值'
ws['A2'] = '分级'; ws['B2'] = '值域'
ws['C2'] = '数量/个'; ws['D2'] = '占比'
ws['E2'] = '面积/亩'; ws['F2'] = '占比'
acid_levels = ['弱酸化', '轻度酸化', '中度酸化', '重度酸化', '其他']
acid_level_headers = ['0.1<ΔpH≤0.3', '0.3<ΔpH≤0.5', '0.5<ΔpH≤1.0', 'ΔpH>1.0', '未酸化']
# --- c. 填充数据 ---
current_row = 3
# 1. 遍历该一级地类下的所有“二级地类”并写入数据
for index,level in enumerate(acid_levels):
ws.cell(row=current_row, column=1).value = level
ws.cell(row=current_row, column=2).value = acid_level_headers[index]
# 在子集中查找当前二级地类的数据行
row_data = df[df['SHFJ'] == level]
# --- 填充单元格的逻辑开始 ---
col_start = 3 # 从第 C 列开始填充
# 检查是否找到了该土属的数据
if not row_data.empty:
data_series = row_data.iloc[0]
# 1. 构建要从 data_series 中查找的列名
sample_col = f'样点数'
sample_pct_col = f'样点占比'
area_col = f'制图面积_亩'
area_pct_col = f'面积占比'
# 2. 从 data_series 中安全地获取值
sample_val = data_series.get(sample_col, 0)
sample_pct_val = data_series.get(sample_pct_col, 0)
area_val = data_series.get(area_col, 0)
area_pct_val = data_series.get(area_pct_col, 0)
# 3. 将获取到的值填入单元格
ws.cell(row=current_row, column=col_start).value = f"{sample_val:.0f}" if sample_val > 0 else "-"
# 占比/%
ws.cell(row=current_row, column=col_start + 1).value = f"{sample_pct_val:.2f}%" if sample_val > 0 else "-"
# 制图面积/亩
ws.cell(row=current_row, column=col_start + 2).value = f"{area_val:.0f}" if area_val > 0 else "-"
# 占比/%
ws.cell(row=current_row, column=col_start + 3).value = f"{area_pct_val:.2f}%" if area_val > 0 else "-"
# 移动到下一个酸化等级的起始列
col_start += 2
else:
for _ in range(4):
ws.cell(row=current_row, column=col_start).value = "-"
col_start += 1
current_row += 1
# 合计单元格填充
mask = df["SHFJ"].isin(acid_levels)
df_acid = df[mask]
weighted_avg = (df_acid["dPH"] * df_acid["样点数"]).sum() / df_acid["样点数"].sum()
mean_msk = df_mean["酸化程度_制图"].isin(acid_levels)
df_mean_acid = df_mean[mean_msk]
weighted_mean = (df_mean_acid["制图均值"] * df_mean_acid["制图样点数"]).sum() / df_mean_acid["制图样点数"].sum()
ws.merge_cells('C9:F9')
ws.merge_cells('C10:F10')
ws['C8'] = df[df['SHFJ'].isin(acid_levels)]['样点数'].sum()
ws['D8'] = f"{df[df['SHFJ'].isin(acid_levels)]['样点占比'].sum():.2f}%"
ws['E8'] = f"{df[df['SHFJ'].isin(acid_levels)]['制图面积_亩'].sum():.0f}"
ws['F8'] = f"{df[df['SHFJ'].isin(acid_levels)]['面积占比'].sum():.2f}%"
ws['C9'] = f"{weighted_avg:.2f}" # type: ignore
ws['C10'] = f"{weighted_mean:.2f}"
# --- a. 定义样式 (不变) ---
header_font = Font(name='宋体', size=11)
# --- d. 应用样式和调整列宽 (最终健壮版) ---
if current_row > 1: # 确保有数据才应用样式
ExcelStyleUtils.set_style(ws, f'A1:F10')
ExcelStyleUtils.set_style(ws, f'A1:F2', header_font)
print("正在自动调整列宽...")
# 设置列宽
ExcelStyleUtils.auto_adjust_column_width(ws)
# --- e. 保存文件 ---
wb.save(output_path)
print("Excel 报告生成成功！")
def main(gdb_path, xzq_features, ph_features, dltb_features, sh_ph_tif, output_path,target_areas_dict:dict):
try:
# --- 1. 用户配置 ---
# 输出配置
output_excel_path = os.path.join(output_path,"乡镇街道酸化统计表.xlsx") # 生成的Excel报告文件路径
# 设置工作空间和变量
arcpy.env.workspace = gdb_path
arcpy.env.overwriteOutput = True
sample_table_name = "历史样点PH信息_Table" # 图2: 样点信息表名
in_zone_feature = xzq_features # 规划分区图层
in_class_feature = ph_features # 已重分类好的酸化PH图层
dltb_class_feature = dltb_features
in_value_raster = sh_ph_tif # 赋值栅格
out_feature_class = "最小面积统计单元"
out_table_area = r"行政区划_酸化面积表" # 输出的交集表名
out_table_mean = r"行政区划_酸化均值表" # 输出的均值表名
print("开始处理数据...")
if not arcpy.Exists(out_feature_class):
# 判断输入表是否存在SHFJ字段
try:
arcpy.management.CalculateField(in_class_feature, "SHFJ", "calculate_shfj(!gridcode!)", "PYTHON3", codeblock_cal_shfj)
except Exception as e:
print(f"计算SHFJ字段时发生错误: {e}")
arcpy.analysis.Intersect(
in_features=[dltb_class_feature, in_class_feature],
out_feature_class=out_feature_class,
join_attributes="ALL",
output_type="INPUT"
)
if not arcpy.Exists(out_table_area):
# 1.用arcpy.analysis.TabulateIntersection进行交集制表
arcpy.analysis.TabulateIntersection(
in_zone_feature,
["XZQMC"],
out_feature_class,
out_table_area,
"SHFJ",
out_units="SQUARE_METERS",
)
if not arcpy.Exists(out_table_mean):
# 2.用arcpy.sa.ZonalStatisticsAsTable进行区域统计
arcpy.sa.ZonalStatisticsAsTable(
in_zone_feature, "XZQMC", in_value_raster, out_table_mean, "DATA", "MEAN"
)
# 计算按地类平差后的各酸化等级面积
if arcpy.Exists(out_feature_class):
df = pd.DataFrame(arcpy.da.TableToNumPyArray(out_feature_class, ["YJDL", "SHFJ", "Shape_Area"]))
df_area = df.groupby(["YJDL", "SHFJ"]).agg({"Shape_Area": "sum"}).reset_index()
yjdl_area = df_area.groupby(['YJDL'])['Shape_Area'].sum().reset_index()
landuse_types = {'耕地':'01', '园地':'02', '林地':'03', '草地':'04', '其他':'12'}
adjustment_factors = []
for _, row in yjdl_area.iterrows():
yjdl = row['YJDL']
original_total = row['Shape_Area'] * 0.0015
target_total = target_areas_dict.get(landuse_types[yjdl], original_total)
adjustment_factor = target_total / original_total
adjustment_factors.append({
'YJDL': yjdl,
'平差系数': adjustment_factor
})
factor_df = pd.DataFrame(adjustment_factors)
df_with_factors = df_area.merge(factor_df[['YJDL', '平差系数']], on='YJDL')
df_with_factors['AREA_MU'] = df_with_factors['Shape_Area'] * df_with_factors['平差系数'] * 0.0015
# print(df_with_factors)
# 生成表5.4的面积统计Excel报告
final_area_dataframe = process_data_for_table5_4(gdb_path, out_table_area, df_with_factors)
write_to_excel_table5_4(final_area_dataframe, output_excel_path)
# 生成表5.3的均值统计Excel报告
final_mean_dataframe = process_data_for_table5_7(gdb_path, out_table_mean, sample_table_name)
write_to_excel_table5_7(final_mean_dataframe, output_excel_path.replace(".xlsx", "_mean.xlsx"))
# 生成总表5.2的统计Excel报告
final_dataframe = process_data_for_table5_2(gdb_path, out_table_area, sample_table_name, df_with_factors)
write_to_excel_table5_2(final_dataframe, final_mean_dataframe, output_excel_path.replace(".xlsx", "_total.xlsx"))
return df_with_factors
except Exception as e:
print(f"\n处理过程中发生严重错误: {e}")
import traceback
traceback.print_exc()
finally:
import gc
gc.collect()
# --- 4. 主程序入口 ---
# if __name__ == "__main__":
# main()