AcidMap/Water/Python_codes/Counter.py

import os
import pandas as pd
from pyproj import Transformer
from shapely.geometry import Point
import numpy as np

# 设置工作目录
os.chdir(r"D:\土壤-大创\Irrigation_Water\Raster")

# 读取耕地中心点CSV文件
csv_cultivated_land = r"..\Data\四县三用地&水田.csv"
df_land = pd.read_csv(csv_cultivated_land)

# 确保CSV文件包含必要的字段
required_fields = ['中心点经度', '中心点纬度', 'DLMC']
for field in required_fields:
    if field not in df_land.columns:
        raise ValueError(f"耕地CSV文件缺少必要的字段: {field}")

# 筛选出旱地
land_type = "旱地"  # 避免使用str作为变量名，因为它是Python内置函数
df_land = df_land[df_land['DLMC'] == land_type]

# 定义坐标系转换（假设输入已是WGS84，根据实际情况修改）
source_crs = "EPSG:4326"
transformer = Transformer.from_crs(source_crs, "EPSG:4326", always_xy=True)

# 读取采样点数据
xls_file = r"..\Data\Irrigation_water_SamplingPoint.xlsx"
df_xls = pd.read_excel(xls_file)

# 设置系数
A = 711 * 0.524
B = 427 * 0.599
C = 200 * 0.7

# 根据土地类型选择系数
coefficient = 0.0
if land_type == "水浇地":
    coefficient = B
elif land_type == "水田":
    coefficient = A
else:
    coefficient = C

# 存储结果的列表
data_rows = []

# 遍历耕地中心点
for index, row in df_land.iterrows():
    # 获取正确的坐标字段
    lon = row['中心点经度']
    lat = row['中心点纬度']
    
    # 坐标转换（如果需要）
    if source_crs != "EPSG:4326":
        lon, lat = transformer.transform(lon, lat)
    
    # 创建中心点几何对象
    center_point = Point(lon, lat)
    
    # 计算到所有采样点的距离
    distances = df_xls.apply(
        lambda x: center_point.distance(Point(x['经度'], x['纬度'])), 
        axis=1
    )
    
    # 找到最近的采样点
    nearest_index = distances.idxmin()
    nearest_sample = df_xls.loc[nearest_index]
    
    # 计算Cd值
    cd_value = nearest_sample['Cd (ug/L)'] * coefficient
    
    # 添加到结果列表
    data_rows.append({
        'DLMC': row['DLMC'],
        '中心点经度': lon,
        '中心点纬度': lat,
        'Cd (ug/L)': cd_value
    })

# 创建DataFrame并选择指定列
csv_data = pd.DataFrame(data_rows)

# 保存结果
csv_file_path = r"D:\土壤-大创\Water\Data\旱地_Cd.csv"
csv_data.to_csv(csv_file_path, index=False, encoding='utf-8-sig')

print(f"CSV文件已保存至: {csv_file_path}")