AcidMap/Water/Python_codes/Calculate.py

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

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

# 读取 Shp 文件
shp_file = r"..\Raster\四县三种用电.shp"
gdf_shp = gpd.read_file(shp_file)

str = "旱地"
# 筛选出字段 DLMC 为指定类型的面要素
gdf_shp = gdf_shp[gdf_shp['DLMC'] == str]

# 定义坐标系转换对象（从原始坐标系转换为WGS84）
transformer = Transformer.from_crs(gdf_shp.crs, "EPSG:4326", always_xy=True)

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

# 新建重金属含量字段，指定数据类型为双精度浮点数
gdf_shp['Cd'] = pd.Series(dtype=np.float64)
gdf_shp['Hg'] = pd.Series(dtype=np.float64)
gdf_shp['Cr'] = pd.Series(dtype=np.float64)
gdf_shp['Pb'] = pd.Series(dtype=np.float64)
gdf_shp['As'] = pd.Series(dtype=np.float64)

# 根据面要素类型设置系数
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 gdf_shp.iterrows():
    # 获取面要素的中心并转换坐标
    center_original = row['geometry'].centroid
    center_lon, center_lat = transformer.transform(center_original.x, center_original.y)
    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
    
    # 更新面要素的重金属含量
    gdf_shp.at[index, 'Cd'] = cd_value
    gdf_shp.at[index, 'Hg'] = hg_value
    gdf_shp.at[index, 'Cr'] = cr_value
    gdf_shp.at[index, 'Pb'] = pb_value
    gdf_shp.at[index, 'As'] = as_value
    
    # 将数据行添加到列表中
    data_rows.append({
        '面要素ID': index,
        'DLMC': row['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文件 - 修复路径中的转义字符问题
csv_file_path = os.path.join('D:\土壤-大创\Water\Data\四县三用地&旱地.csv')
csv_data.to_csv(csv_file_path, index=False, encoding='utf-8-sig')

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