AcidMap/Water/Python_codes/Calculate.py

import os
import sys
import pandas as pd
from shapely.geometry import Point
import numpy as np
from pathlib import Path

# 添加项目根目录到Python路径
project_root = Path(__file__).parent.parent.parent
sys.path.append(str(project_root))

# 导入数据库服务
from app.services.land_use_service import land_use_service

# 设置要处理的土地类型
str = "旱地"

# 从数据库获取土地利用数据，替代原来的SHP文件读取
print(f"从数据库获取 '{str}' 类型的土地利用数据...")
land_centers_df = land_use_service.get_land_centers_for_processing(str)

if land_centers_df is None or land_centers_df.empty:
    print(f"数据库中没有找到 '{str}' 类型的土地利用数据")
    exit(1)

print(f"从数据库获取到 {len(land_centers_df)} 个 '{str}' 类型的土地数据")

# 读取 xls 文件
xls_file = r"..\Data\Irrigation_water_SamplingPoint.xlsx"
df_xls = pd.read_excel(xls_file)

# 不再需要为GeoDataFrame添加字段，直接处理数据即可

# 根据面要素类型设置系数
A = 711 * 0.524
B = 427 * 0.599
C = 200 * 0.7

Num = 0.0
if str == "水浇地":
    Num = B
elif str == "水田":
    Num = A
else:
    Num = C

# 创建一个空列表来存储数据行
data_rows = []

# 遍历面要素，使用数据库中的中心点坐标
for index, row in land_centers_df.iterrows():
    # 直接使用数据库中已经转换好的WGS84坐标
    center_lon = row['center_lon']
    center_lat = row['center_lat']
    center_transformed = Point(center_lon, center_lat)
    
    # 计算面要素中心与采样点的距离
    distances = df_xls.apply(lambda x: center_transformed.distance(Point(x['经度'], x['纬度'])), axis=1)
    
    # 找到距离最近的采样点
    min_distance_index = distances.idxmin()
    nearest_sample = df_xls.loc[min_distance_index]
    
    # 计算最近采样点的重金属含量，并乘以系数Num
    cd_value = nearest_sample['Cd (ug/L)'] * Num
    hg_value = nearest_sample['Hg (ug/L)'] * Num
    cr_value = nearest_sample['Cr (ug/L)'] * Num
    pb_value = nearest_sample['Pb (ug/L)'] * Num
    as_value = nearest_sample['As (ug/L)'] * Num
    
    # 将数据行添加到列表中
    data_rows.append({
        '面要素ID': row['gid'],  # 使用数据库中的gid
        'DLMC': row['dlmc'],    # 使用数据库中的dlmc字段
        '中心点经度': center_lon,
        '中心点纬度': center_lat,
        '最近采样点ID': nearest_sample.name,
        '采样点经度': nearest_sample['经度'],
        '采样点纬度': nearest_sample['纬度'],
        'Cd (ug/L)': cd_value,
        'Hg (ug/L)': hg_value,
        'Cr (ug/L)': cr_value,
        'Pb (ug/L)': pb_value,
        'As (ug/L)': as_value
    })

# 从列表创建DataFrame
csv_data = pd.DataFrame(data_rows)

# 保存结果为CSV文件
output_dir = Path(__file__).parent.parent / "Data"
output_dir.mkdir(exist_ok=True)
csv_file_path = output_dir / f"四县三用地&{str}_DB.csv"
csv_data.to_csv(csv_file_path, index=False, encoding='utf-8-sig')

print(f"CSV文件已保存至: {csv_file_path}")
print(f"总共处理了 {len(csv_data)} 个 {str} 面要素")