AcidMap/app/services/water_service.py

import os
import geopandas as gpd
import pandas as pd
from pyproj import Transformer
from shapely.geometry import Point
import rasterio
from typing import Optional, Dict, Any
from datetime import datetime
import numpy as np
import logging
import shutil
import sys

# 导入MappingUtils
from ..utils.mapping_utils import MappingUtils, csv_to_raster_workflow, dataframe_to_raster_workflow
# 导入土地利用服务
from .land_use_service import land_use_service
# 导入边界服务
from .admin_boundary_service import get_boundary_gdf_by_name
from ..database import SessionLocal
import geopandas as gpd
import tempfile

# 配置日志
logger = logging.getLogger(__name__)

# 设置全局字体
import matplotlib.pyplot as plt

plt.rcParams['font.family'] = 'Arial'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']

# 设定常用的colormap（与mapping_utils保持一致）
COLORMAPS = {
    "yellow_orange_brown": ['#FFFECE', '#FFF085', '#FEBA17', '#BE3D2A', '#74512D', '#4E1F00'],
    "blue_series": ['#F6F8D5', '#98D2C0', '#4F959D', '#205781', '#143D60', '#2A3335'],
    "yellow_green": ['#FFEFC8', '#F8ED8C', '#D3E671', '#89AC46', '#5F8B4C', '#355F2E'],
    "green_brown": ['#F0F1C5', '#BBD8A3', '#6F826A', '#BF9264', '#735557', '#604652'],
    "yellow_pink_purple": ['#FCFAEE', '#FBF3B9', '#FFDCCC', '#FDB7EA', '#B7B1F2', '#8D77AB'],
    "green_yellow_red_purple": ['#15B392', '#73EC8B', '#FFEB55', '#EE66A6', '#D91656', '#640D5F']
}

# 建立中文到英文的colormap映射
COLORMAP_MAPPING = {
    "黄-橙-棕": "yellow_orange_brown",
    "蓝色系": "blue_series",
    "淡黄-草绿": "yellow_green",
    "绿色-棕色": "green_brown",
    "黄-粉-紫": "yellow_pink_purple",
    "绿-黄-红-紫": "green_yellow_red_purple"
}


# 主路径配置
def get_base_dir():
    """获取基础目录路径"""
    if getattr(sys, 'frozen', False):
        # 打包后的可执行文件
        return os.path.dirname(sys.executable)
    else:
        # 脚本运行模式
        return os.path.dirname(os.path.abspath(__file__))


def get_boundary_gdf_from_database(area: str, level: str) -> Optional[gpd.GeoDataFrame]:
    """
    直接从数据库获取边界数据作为GeoDataFrame
    
    @description: 与其他模块保持一致的边界获取方法
    @param area: 地区名称
    @param level: 行政层级（county | city | province）
    @returns: 边界GeoDataFrame或None
    """
    try:
        with SessionLocal() as db:
            norm_area = area.strip()
            boundary_gdf = get_boundary_gdf_by_name(db, norm_area, level=level)
            if boundary_gdf is not None:
                logger.info(f"从数据库获取边界数据: {norm_area} ({level})")
            return boundary_gdf
                
    except Exception as e:
        logger.warning(f"从数据库获取边界数据失败: {str(e)}")
        
    return None





# 土地数据处理函数
def process_land_data(land_type, coefficient_params=None, save_csv=True):
    """处理土地类型数据并生成清洗后的简化数据"""
    base_dir = get_base_dir()
    xls_file = os.path.join(base_dir, "..", "static", "water", "Data", "Irrigation_water_SamplingPoint.xlsx")

    logger.info(f"处理土地类型: {land_type}")
    logger.info(f"Excel文件: {xls_file}")

    # 从数据库获取土地利用数据，替代原来的SHP文件读取
    logger.info(f"从数据库获取 '{land_type}' 类型的土地利用数据...")
    land_centers_df = land_use_service.get_land_centers_for_processing(land_type)
    
    if land_centers_df is None or land_centers_df.empty:
        logger.warning(f"数据库中没有找到 '{land_type}' 类型的土地利用数据")
        return None, None, None

    logger.info(f"从数据库获取到 {len(land_centers_df)} 个 '{land_type}' 类型的土地数据")

    # 读取Excel采样点数据
    if not os.path.exists(xls_file):
        logger.error(f"采样点Excel文件不存在: {xls_file}")
        return None, None, None
        
    df_xls = pd.read_excel(xls_file)
    logger.info(f"读取到 {len(df_xls)} 个采样点数据")

    # 设置土地类型系数
    default_params = {
        "水田": (711, 0.524),
        "水浇地": (427, 0.599),
        "旱地": (200, 0.7)
    }

    params = coefficient_params or default_params
    param1, param2 = params.get(land_type, (200, 0.7))
    Num = param1 * param2
    logger.info(f"系数: {param1} * {param2} = {Num}")

    # 处理每个面要素，使用数据库中的中心点坐标
    cd_values = []
    centers = []
    
    for index, row in land_centers_df.iterrows():
        center_lon = row['center_lon']
        center_lat = row['center_lat']
        centers.append((center_lon, center_lat))

        # 计算到所有采样点的距离
        distances = df_xls.apply(
            lambda x: Point(center_lon, center_lat).distance(Point(x['经度'], x['纬度'])),
            axis=1
        )
        min_idx = distances.idxmin()
        nearest = df_xls.loc[min_idx]

        # 计算Cd含量值
        cd_value = nearest['Cd (ug/L)'] * Num
        cd_values.append(cd_value)

    # 创建简化数据DataFrame
    simplified_data = pd.DataFrame({
        'lon': [c[0] for c in centers],
        'lat': [c[1] for c in centers],
        'Prediction': cd_values
    })

    # 应用3σ原则检测异常值
    mean_value = simplified_data['Prediction'].mean()
    std_value = simplified_data['Prediction'].std()
    lower_bound = mean_value - 3 * std_value
    upper_bound = mean_value + 3 * std_value

    logger.info(f"Cd含量 - 平均值: {mean_value:.4f}, 标准差: {std_value:.4f}")
    logger.info(f"检测范围: 下限 = {lower_bound:.4f}, 上限 = {upper_bound:.4f}")

    # 识别异常值
    outliers = simplified_data[
        (simplified_data['Prediction'] < lower_bound) |
        (simplified_data['Prediction'] > upper_bound)
        ]
    logger.info(f"检测到异常值数量: {len(outliers)}")

    # 创建清洗后的数据
    cleaned_data = simplified_data[
        (simplified_data['Prediction'] >= lower_bound) &
        (simplified_data['Prediction'] <= upper_bound)
        ]
    logger.info(f"清洗后数据点数: {len(cleaned_data)}")

    # 可选：保存CSV文件用于兼容性和调试
    cleaned_csv_path = None
    if save_csv:
        # 创建输出目录
        data_dir = os.path.join(base_dir, "..", "static", "water", "Data")
        os.makedirs(data_dir, exist_ok=True)
        logger.info(f"数据目录: {data_dir}")
        
        cleaned_csv_path = os.path.join(data_dir, f"中心点经纬度与预测值&{land_type}_清洗.csv")
        cleaned_data.to_csv(cleaned_csv_path, index=False, encoding='utf-8-sig')
        logger.info(f"保存CSV: {cleaned_csv_path}")
    else:
        logger.info("跳过CSV文件生成（内存处理优化）")

    return cleaned_data, Num, cleaned_csv_path  # 返回DataFrame, 系数, 和CSV路径（如果生成了）


# 可视化函数（完全使用统一的MappingUtils接口）
def mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path, output_size=8):
    """
    生成栅格地图可视化
    
    @description: 使用统一的MappingUtils接口生成栅格地图
    @param shp_path: 边界shapefile路径
    @param tif_path: 栅格数据文件路径  
    @param color_map_name: 色彩方案名称（中文）
    @param title_name: 地图标题
    @param output_path: 输出图片路径
    @param output_size: 图片尺寸
    @returns: 输出图片路径或None（失败时）
    """
    try:
        logger.info(f"生成栅格地图: {title_name}")

        # 确保输出目录存在
        os.makedirs(os.path.dirname(output_path), exist_ok=True)

        # 创建MappingUtils实例
        mapper = MappingUtils()

        # 转换颜色方案名称
        colormap_key = COLORMAP_MAPPING.get(color_map_name, "green_yellow_red_purple")

        # 直接调用统一的绘图接口，无需额外包装
        generated_path = mapper.create_raster_map(
            shp_path=shp_path,
            tif_path=tif_path,
            output_path=os.path.splitext(output_path)[0],  # 去掉扩展名，MappingUtils会自动添加.jpg
            colormap=colormap_key,
            title=title_name,
            output_size=output_size,
            figsize=None,
            dpi=300,
            enable_interpolation=False,
            interpolation_method='linear'
        )

        # 如果生成的路径与期望路径不同，则重命名
        if generated_path != output_path and os.path.exists(generated_path):
            shutil.move(generated_path, output_path)
            logger.info(f"重命名图片: {generated_path} -> {output_path}")
            return output_path
        
        return generated_path

    except Exception as e:
        logger.error(f"栅格地图生成失败: {str(e)}")
        return None


def plot_tif_histogram(file_path, output_path, figsize=(8, 8),
                       xlabel='Cd(ug/L)', ylabel='frequency density',
                       title='Irrigation water input flux'):
    """
    生成TIFF文件的直方图
    
    @description: 使用统一的MappingUtils接口生成直方图
    @param file_path: TIFF文件路径
    @param output_path: 输出图片路径
    @param figsize: 图片尺寸元组
    @param xlabel: X轴标签
    @param ylabel: Y轴标签  
    @param title: 图表标题
    @returns: 输出图片路径或None（失败时）
    """
    try:
        logger.info(f"生成直方图: {file_path}")

        # 创建MappingUtils实例
        mapper = MappingUtils()

        # 直接调用统一的绘图接口，无需额外包装
        generated_path = mapper.create_histogram(
            file_path=file_path,
            save_path=output_path,
            figsize=figsize,
            xlabel=xlabel,
            ylabel=ylabel,
            title=title,
            bins=100,
            dpi=300
        )

        return generated_path

    except Exception as e:
        logger.error(f"直方图生成失败: {str(e)}")
        return None


def process_land_to_visualization(land_type, coefficient_params=None,
                                  color_map_name="绿-黄-红-紫",
                                  output_size=8,
                                  area: Optional[str] = None,
                                  level: Optional[str] = None,
                                  enable_interpolation: Optional[bool] = True,
                                  interpolation_method: Optional[str] = "linear",
                                  resolution_factor: Optional[float] = 4.0,
                                  save_csv: Optional[bool] = True):
    """
    完整的土地数据处理可视化流程（使用统一的MappingUtils接口，支持动态边界和插值控制）
    
    @description: 水样数据处理与可视化的完整工作流程，已优化为完全使用统一的绘图接口，
                  支持通过area和level参数动态控制地图边界，支持插值控制参数
    
    工作流程:
    1. 生成清洗后CSV - 从数据库获取土地利用数据并计算Cd含量
    2. 获取动态边界 - 根据area和level参数从数据库获取边界数据
    3. 使用csv_to_raster_workflow转换为GeoTIFF - 调用统一的栅格转换工具，支持插值控制
    4. 生成栅格地图 - 直接调用MappingUtils.create_raster_map()，支持动态边界
    5. 生成直方图 - 直接调用MappingUtils.create_histogram()

    @param land_type: 土地类型（水田、旱地、水浇地）
    @param coefficient_params: 可选的系数参数字典
    @param color_map_name: 色彩方案名称（中文）
    @param output_size: 输出图片尺寸
    @param area: 可选的地区名称，如"乐昌市"，用于动态控制地图边界
    @param level: 可选的行政层级，county | city | province
    @param enable_interpolation: 是否启用空间插值，默认True
    @param interpolation_method: 插值方法，nearest | linear | cubic，默认linear
    @param resolution_factor: 分辨率因子，默认4.0，越大分辨率越高
    @param save_csv: 是否生成CSV文件，默认True
    @returns: 包含所有生成文件路径的元组
    """
    base_dir = get_base_dir()
    logger.info(f"开始处理: {land_type}")

    # 1. 生成清洗后的数据（内存处理，可选择是否保存CSV）
    cleaned_data, used_coeff, cleaned_csv_path = process_land_data(
        land_type, 
        coefficient_params, 
        save_csv=save_csv  # 根据参数决定是否生成CSV文件
    )
    if cleaned_data is None:
        logger.error(f"处理土地数据失败: {land_type}")
        return None, None, None, None, None, None

    # 2. 获取动态边界数据（提前获取用于栅格转换）
    boundary_gdf = None
    if area and level:
        logger.info(f"尝试获取动态边界: area='{area}', level='{level}'")
        boundary_gdf = get_boundary_gdf_from_database(area, level)
        if boundary_gdf is not None:
            logger.info(f"成功获取动态边界用于栅格转换: {area} ({level}) - 包含 {len(boundary_gdf)} 个几何体")
            logger.info(f"边界数据范围: {boundary_gdf.bounds.to_dict('records')[0] if len(boundary_gdf) > 0 else 'N/A'}")
        else:
            logger.warning(f"未能获取动态边界，栅格转换将使用完整范围: {area} ({level})")
            logger.warning("可能的原因: 1) 数据库中没有该地区数据 2) 地区名称拼写错误 3) level参数不正确")
    else:
        logger.info("未指定area和level参数，尝试获取默认边界（乐昌市）")
        boundary_gdf = get_boundary_gdf_from_database("乐昌市", "county")
        if boundary_gdf is not None:
            logger.info(f"成功获取默认边界（乐昌市）用于栅格转换 - 包含 {len(boundary_gdf)} 个几何体")
        else:
            logger.info("未能获取默认边界，栅格转换将使用数据点范围")
            
    # 3. 使用dataframe_to_raster_workflow转换为GeoTIFF（内存处理，支持动态边界）
    raster_dir = os.path.join(base_dir, "..", "static", "water", "Raster")
    template_tif = os.path.join(raster_dir, "meanTemp.tif")
    output_dir = raster_dir  # 直接使用raster目录作为输出目录

    # 调用dataframe_to_raster_workflow，支持动态边界和插值控制（内存处理优化）
    workflow_result = dataframe_to_raster_workflow(
        df=cleaned_data,  # 直接使用DataFrame
        template_tif=template_tif,
        output_dir=output_dir,
        boundary_gdf=boundary_gdf,  # 直接使用GeoDataFrame
        resolution_factor=resolution_factor,  # 可配置的分辨率因子
        interpolation_method=interpolation_method,  # 可配置的插值方法
        field_name='Prediction',
        lon_col=0,  # DataFrame中经度列索引
        lat_col=1,  # DataFrame中纬度列索引
        value_col=2,  # DataFrame中数值列索引
        enable_interpolation=enable_interpolation  # 可配置的插值开关
    )

    # 获取输出的栅格文件路径
    output_tif = workflow_result['raster']
    logger.info(f"生成栅格文件: {output_tif}")

    # 4. 生成栅格地图（直接使用统一的MappingUtils接口，支持动态边界）
    map_output = os.path.join(raster_dir, f"{land_type}_Cd含量地图.jpg")
    
    try:
        # 创建MappingUtils实例并直接调用
        mapper = MappingUtils()
        
        # 转换颜色方案名称
        colormap_key = COLORMAP_MAPPING.get(color_map_name, "green_yellow_red_purple")
        
        # 统一使用boundary_gdf，不再依赖Shapefile
        map_result = mapper.create_raster_map(
            shp_path=None,  # 不使用Shapefile
            tif_path=output_tif,
            output_path=os.path.splitext(map_output)[0],  # 去掉扩展名
            colormap=colormap_key,
            title="Irrigation water input flux",
            output_size=output_size,
            dpi=300,
            enable_interpolation=enable_interpolation,  # 可配置的插值开关
            interpolation_method=interpolation_method,  # 可配置的插值方法
            boundary_gdf=boundary_gdf  # 使用从数据库获取的边界数据，如果为None则使用数据点范围
        )
        
        # 如果生成的路径与期望不同，则重命名
        if map_result != map_output and os.path.exists(map_result):
            shutil.move(map_result, map_output)
            logger.info(f"栅格地图已生成: {map_output}")
            
    except Exception as e:
        logger.error(f"栅格地图生成失败: {str(e)}")

    # 5. 生成直方图（直接使用统一的MappingUtils接口）
    hist_output = os.path.join(raster_dir, f"{land_type}_Cd含量直方图.jpg")
    
    try:
        # 直接调用统一绘图接口生成直方图
        hist_result = mapper.create_histogram(
            file_path=output_tif,
            save_path=hist_output,
            figsize=(8, 8),
            xlabel='Cd(ug/L)',
            ylabel='frequency density',
            title='Irrigation water input flux',
            bins=100,
            dpi=300
        )
        
        if hist_result:
            logger.info(f"直方图已生成: {hist_output}")
            
    except Exception as e:
        logger.error(f"直方图生成失败: {str(e)}")

    # 使用实际的CSV路径（如果生成了）或生成兼容性路径
    if cleaned_csv_path:
        cleaned_csv = cleaned_csv_path
    else:
        # 为了兼容性，生成路径（即使文件不存在）
        data_dir = os.path.join(base_dir, "..", "static", "water", "Data")
        cleaned_csv = os.path.join(data_dir, f"中心点经纬度与预测值&{land_type}_清洗.csv")
    
    return cleaned_csv, workflow_result['shapefile'], output_tif, map_output, hist_output, used_coeff


def get_land_statistics(land_type: str) -> Optional[Dict[str, Any]]:
    """
    获取指定土地类型的Cd预测结果统计信息

    @param {str} land_type - 土地类型（水田、旱地或水浇地）
    @returns {Optional[Dict[str, Any]]} 统计信息，如果没有数据则返回None
    """
    try:
        logger.info(f"获取土地类型统计信息: {land_type}")

        # 获取基础目录
        base_dir = get_base_dir()
        data_dir = os.path.join(base_dir, "..", "static", "water", "Data")

        # 构建数据文件路径
        data_file = os.path.join(data_dir, f"中心点经纬度与预测值&{land_type}_清洗.csv")
        logger.info(f"数据文件路径: {data_file}")

        if not os.path.exists(data_file):
            logger.warning(f"未找到土地类型数据文件: {data_file}")
            return None

        # 读取预测数据
        df = pd.read_csv(data_file)
        logger.info(f"成功读取数据文件，包含 {len(df)} 行数据")

        # 检查必要的列是否存在
        if 'Prediction' not in df.columns:
            logger.warning("数据文件中缺少'Prediction'列")
            return None

        predictions = df['Prediction']

        # 计算基础统计信息
        stats = {
            "土地类型": land_type,
            "数据更新时间": datetime.fromtimestamp(os.path.getmtime(data_file)).strftime('%Y-%m-%d %H:%M:%S'),
            "数据点总数": len(predictions),
            "均值": float(predictions.mean()),
            "中位数": float(predictions.median()),
            "标准差": float(predictions.std()),
            "最小值": float(predictions.min()),
            "最大值": float(predictions.max()),
            "25%分位数": float(predictions.quantile(0.25)),
            "75%分位数": float(predictions.quantile(0.75)),
            "偏度": float(predictions.skew()),
            "峰度": float(predictions.kurtosis())
        }

        return stats

    except Exception as e:
        logger.error(f"获取土地类型统计信息失败: {str(e)}")
        return None


# 主函数
def main():
    """主函数用于测试完整的处理流程"""
    logger.info("=" * 50)
    logger.info("土地数据处理与可视化系统")
    logger.info("=" * 50)

    try:
        # 处理水田数据
        logger.info("\n===== 处理水田数据 =====")
        results = process_land_to_visualization("水田")
        if results and all(results):
            cleaned_csv, shapefile, tif_file, map_file, hist_file, used_coeff = results
            logger.info(f"清洗后CSV: {cleaned_csv}")
            logger.info(f"Shapefile: {shapefile}")
            logger.info(f"GeoTIFF文件: {tif_file}")
            logger.info(f"地图图片: {map_file}")
            logger.info(f"直方图: {hist_file}")
            logger.info(f"使用系数: {used_coeff}")
        else:
            logger.error("水田数据处理失败")

        # 处理旱地数据（使用自定义参数）
        logger.info("\n===== 处理旱地数据 =====")
        custom_params = {"旱地": (220, 0.65)}
        results = process_land_to_visualization(
            "旱地",
            coefficient_params=custom_params,
            color_map_name="蓝色系"
        )
        if results and all(results):
            cleaned_csv, shapefile, tif_file, map_file, hist_file, used_coeff = results
            logger.info(f"清洗后CSV: {cleaned_csv}")
            logger.info(f"Shapefile: {shapefile}")
            logger.info(f"GeoTIFF文件: {tif_file}")
            logger.info(f"地图图片: {map_file}")
            logger.info(f"直方图: {hist_file}")
            logger.info(f"使用系数: {used_coeff}")
        else:
            logger.error("旱地数据处理失败")

    except Exception as e:
        logger.error(f"处理过程中发生错误: {str(e)}", exc_info=True)
    finally:
        logger.info("处理完成")


if __name__ == "__main__":
    main()