AcidMap/app/services/cd_flux_removal_service.py

"""
Cd通量移除计算服务
@description: 提供籽粒移除和秸秆移除的Cd通量计算功能
@author: AcidMap Team
@version: 1.0.0
"""

import logging
import math
import os
import pandas as pd
from datetime import datetime
from typing import Dict, Any, List, Optional
from sqlalchemy.orm import sessionmaker, Session
from sqlalchemy import create_engine, and_
from ..database import SessionLocal, engine
from ..models import FluxCdInputData
from ..models.parameters import Parameters
from ..models.CropCd_output import CropCdOutputData
from ..models.farmland import FarmlandData
from ..utils.mapping_utils import MappingUtils
from .admin_boundary_service import get_boundary_geojson_by_name, get_boundary_gdf_by_name
import geopandas as gpd
from shapely.geometry import shape, Point
import tempfile
import json


class CdFluxRemovalService:
    """
    Cd通量移除计算服务类
    
    @description: 提供基于CropCd_output_data和Parameters表数据的籽粒移除和秸秆移除Cd通量计算功能
    """
    
    def __init__(self):
        """
        初始化Cd通量移除服务
        """
        self.logger = logging.getLogger(__name__)
        # 严格匹配策略：不做名称变体或后缀映射
        
    def calculate_grain_removal_by_area(self, area: str, level: Optional[str] = None) -> Dict[str, Any]:
        """
        根据地区计算籽粒移除Cd通量
        
        @param area: 地区名称（仅用于地图边界，参数固定使用"韶关"）
        @returns: 计算结果字典
        
        计算公式：籽粒移除(g/ha/a) = EXP(LnCropCd) * F11 * 0.5 * 15 / 1000
        """
        try:
            with SessionLocal() as db:
                # 参数固定使用"韶关"，area参数仅用于地图边界
                parameter = db.query(Parameters).filter(Parameters.area == "韶关").first()
                
                if not parameter:
                    return {
                        "success": False,
                        "message": f"未找到韶关地区的参数数据",
                        "data": None
                    }
                
                # 查询CropCd输出数据
                crop_cd_outputs = db.query(CropCdOutputData).all()
                
                if not crop_cd_outputs:
                    return {
                        "success": False,
                        "message": f"未找到CropCd输出数据",
                        "data": None
                    }
                
                # 计算每个样点的籽粒移除Cd通量
                results = []
                for output in crop_cd_outputs:
                    crop_cd_value = math.exp(output.ln_crop_cd)  # EXP(LnCropCd)
                    grain_removal = crop_cd_value * parameter.f11 * 0.5 * 15 / 1000
                    
                    results.append({
                        "farmland_id": output.farmland_id,
                        "sample_id": output.sample_id,
                        "ln_crop_cd": output.ln_crop_cd,
                        "crop_cd_value": crop_cd_value,
                        "f11_yield": parameter.f11,
                        "grain_removal_flux": grain_removal
                    })
                
                # 计算统计信息
                flux_values = [r["grain_removal_flux"] for r in results]
                statistics = {
                    "total_samples": len(results),
                    "mean_flux": sum(flux_values) / len(flux_values),
                    "max_flux": max(flux_values),
                    "min_flux": min(flux_values)
                }
                
                return {
                    "success": True,
                    "message": f"地区 '{area}' 的籽粒移除Cd通量计算成功",
                    "data": {
                        "area": area,
                        "calculation_type": "grain_removal",
                        "formula": "EXP(LnCropCd) * F11 * 0.5 * 15 / 1000",
                        "unit": "g/ha/a",
                        "results": results,
                        "statistics": statistics
                    }
                }
                
        except Exception as e:
            self.logger.error(f"计算地区 '{area}' 的籽粒移除Cd通量失败: {str(e)}")
            return {
                "success": False,
                "message": f"计算失败: {str(e)}",
                "data": None
            }
    
    def calculate_straw_removal_by_area(self, area: str, level: Optional[str] = None) -> Dict[str, Any]:
        """
        根据地区计算秸秆移除Cd通量
        
        @param area: 地区名称（仅用于地图边界，参数固定使用"韶关"）
        @returns: 计算结果字典
        
        计算公式：秸秆移除(g/ha/a) = [EXP(LnCropCd)/(EXP(LnCropCd)*0.76-0.0034)] * F11 * 0.5 * 15 / 1000
        """
        try:
            with SessionLocal() as db:
                # 参数固定使用"韶关"，area参数仅用于地图边界
                parameter = db.query(Parameters).filter(Parameters.area == "韶关").first()
                
                if not parameter:
                    return {
                        "success": False,
                        "message": f"未找到韶关地区的参数数据",
                        "data": None
                    }
                
                # 查询CropCd输出数据
                crop_cd_outputs = db.query(CropCdOutputData).all()
                
                if not crop_cd_outputs:
                    return {
                        "success": False,
                        "message": f"未找到CropCd输出数据",
                        "data": None
                    }
                
                # 计算每个样点的秸秆移除Cd通量
                results = []
                for output in crop_cd_outputs:
                    crop_cd_value = math.exp(output.ln_crop_cd)  # EXP(LnCropCd)
                    
                    # 计算分母：EXP(LnCropCd)*0.76-0.0034
                    denominator = crop_cd_value * 0.76 - 0.0034
                    
                    # 检查分母是否为零或负数，避免除零错误
                    if denominator <= 0:
                        self.logger.warning(f"样点 {output.farmland_id}-{output.sample_id} 的分母值为 {denominator}，跳过计算")
                        continue
                    
                    # 计算秸秆移除Cd通量
                    straw_removal = (crop_cd_value / denominator) * parameter.f11 * 0.5 * 15 / 1000
                    
                    results.append({
                        "farmland_id": output.farmland_id,
                        "sample_id": output.sample_id,
                        "ln_crop_cd": output.ln_crop_cd,
                        "crop_cd_value": crop_cd_value,
                        "denominator": denominator,
                        "f11_yield": parameter.f11,
                        "straw_removal_flux": straw_removal
                    })
                
                if not results:
                    return {
                        "success": False,
                        "message": "所有样点的计算都因分母值无效而失败",
                        "data": None
                    }
                
                # 计算统计信息
                flux_values = [r["straw_removal_flux"] for r in results]
                statistics = {
                    "total_samples": len(results),
                    "mean_flux": sum(flux_values) / len(flux_values),
                    "max_flux": max(flux_values),
                    "min_flux": min(flux_values)
                }
                
                return {
                    "success": True,
                    "message": f"地区 '{area}' 的秸秆移除Cd通量计算成功",
                    "data": {
                        "area": area,
                        "calculation_type": "straw_removal",
                        "formula": "[EXP(LnCropCd)/(EXP(LnCropCd)*0.76-0.0034)] * F11 * 0.5 * 15 / 1000",
                        "unit": "g/ha/a",
                        "results": results,
                        "statistics": statistics
                    }
                }
                
        except Exception as e:
            self.logger.error(f"计算地区 '{area}' 的秸秆移除Cd通量失败: {str(e)}")
            return {
                "success": False,
                "message": f"计算失败: {str(e)}",
                "data": None
            }
    
    def export_results_to_csv(self, results_data: Dict[str, Any], output_dir: str = "app/static/cd_flux") -> str:
        """
        将计算结果导出为CSV文件
        
        @param results_data: 计算结果数据
        @param output_dir: 输出目录
        @returns: CSV文件路径
        """
        try:
            # 确保输出目录存在
            os.makedirs(output_dir, exist_ok=True)
            
            # 生成时间戳
            timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
            
            # 生成文件名
            calculation_type = results_data.get("calculation_type", "flux_removal")
            area = results_data.get("area", "unknown")
            filename = f"{calculation_type}_{area}_{timestamp}.csv"
            csv_path = os.path.join(output_dir, filename)
            
            # 转换为DataFrame
            results = results_data.get("results", [])
            if not results:
                raise ValueError("没有结果数据可导出")
            
            df = pd.DataFrame(results)
            
            # 保存CSV文件
            df.to_csv(csv_path, index=False, encoding='utf-8-sig')
            
            self.logger.info(f"✓ 成功导出结果到: {csv_path}")
            return csv_path
            
        except Exception as e:
            self.logger.error(f"导出CSV文件失败: {str(e)}")
            raise
    
    def get_coordinates_for_results(self, results_data: Dict[str, Any]) -> List[Dict[str, Any]]:
        """
        获取结果数据对应的坐标信息
        
        @param results_data: 计算结果数据
        @returns: 包含坐标的结果列表
        """
        try:
            results = results_data.get("results", [])
            if not results:
                return []

            # 提取成对键，避免 N 次数据库查询
            farmland_sample_pairs = [(r["farmland_id"], r["sample_id"]) for r in results]

            with SessionLocal() as db:
                # 使用 farmland_id 分片查询，避免复合 IN 导致的兼容性与参数数量问题
                wanted_pairs = set(farmland_sample_pairs)
                unique_farmland_ids = sorted({fid for fid, _ in wanted_pairs})

                def chunk_list(items: List[int], chunk_size: int = 500) -> List[List[int]]:
                    return [items[i:i + chunk_size] for i in range(0, len(items), chunk_size)]

                rows: List[FarmlandData] = []
                for id_chunk in chunk_list(unique_farmland_ids, 500):
                    rows.extend(
                        db.query(FarmlandData)
                          .filter(FarmlandData.farmland_id.in_(id_chunk))
                          .all()
                    )

                pair_to_farmland = {
                    (row.farmland_id, row.sample_id): row for row in rows
                }

                coordinates_results: List[Dict[str, Any]] = []
                for r in results:
                    key = (r["farmland_id"], r["sample_id"])
                    farmland = pair_to_farmland.get(key)
                    if farmland is None:
                        continue

                    coord_result = {
                        "farmland_id": r["farmland_id"],
                        "sample_id": r["sample_id"],
                        "longitude": farmland.lon,
                        "latitude": farmland.lan,
                        "flux_value": r.get("grain_removal_flux") or r.get("straw_removal_flux")
                    }
                    coord_result.update(r)
                    coordinates_results.append(coord_result)

                self.logger.info(f"✓ 成功获取 {len(coordinates_results)} 个样点的坐标信息（分片批量查询）")
                return coordinates_results
                
        except Exception as e:
            self.logger.error(f"获取坐标信息失败: {str(e)}")
            raise
    
    def create_flux_visualization(self, area: str, calculation_type: str,
                                results_with_coords: List[Dict[str, Any]],
                                level: str = None,
                                output_dir: str = "app/static/cd_flux",
                                template_raster: str = "app/static/cd_flux/meanTemp.tif",
                                boundary_shp: str = None,
                                colormap: str = "green_yellow_red_purple",
                                resolution_factor: float = 4.0,
                                enable_interpolation: bool = False,
                                cleanup_intermediate: bool = True) -> Dict[str, str]:
        """
        创建Cd通量移除可视化图表
        
        @param area: 地区名称
        @param calculation_type: 计算类型（grain_removal 或 straw_removal）
        @param results_with_coords: 包含坐标的结果数据
        @param output_dir: 输出目录
        @param template_raster: 模板栅格文件路径
        @param boundary_shp: 边界shapefile路径
        @param colormap: 色彩方案
        @param resolution_factor: 分辨率因子
        @param enable_interpolation: 是否启用空间插值
        @returns: 生成的图片文件路径字典
        """
        try:
            if not results_with_coords:
                raise ValueError("没有包含坐标的结果数据")
            
            # 确保输出目录存在
            os.makedirs(output_dir, exist_ok=True)
            generated_files: Dict[str, str] = {}
            
            # 生成时间戳
            timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
            
            # 创建CSV文件用于绘图
            csv_filename = f"{calculation_type}_{area}_temp_{timestamp}.csv"
            csv_path = os.path.join(output_dir, csv_filename)
            
            # 准备绘图数据
            plot_data = []
            for result in results_with_coords:
                plot_data.append({
                    "longitude": result["longitude"],
                    "latitude": result["latitude"],
                    "flux_value": result["flux_value"]
                })
            
            # 保存为CSV
            df = pd.DataFrame(plot_data)
            df.to_csv(csv_path, index=False, encoding='utf-8-sig')
            
            # 初始化绘图工具
            mapper = MappingUtils()
            
            # 生成输出文件路径
            map_output = os.path.join(output_dir, f"{calculation_type}_{area}_map_{timestamp}")
            histogram_output = os.path.join(output_dir, f"{calculation_type}_{area}_histogram_{timestamp}")
            
            # 检查模板文件是否存在
            if not os.path.exists(template_raster):
                self.logger.warning(f"模板栅格文件不存在: {template_raster}")
                template_raster = None
            
            # 动态获取边界数据（严格使用指定层级）
            if not level:
                raise ValueError("必须提供行政层级 level：county | city | province")
            
            # 优化：直接从数据库获取边界GeoDataFrame，避免创建临时shapefile文件
            # 这样可以减少磁盘I/O操作和临时文件管理的开销
            boundary_gdf = self._get_boundary_gdf_from_database(area, level)
            boundary_shp = None  # 不再需要临时边界文件
            
            if boundary_gdf is None:
                self.logger.warning(f"未找到地区 '{area}' 的边界数据，将不使用边界裁剪")
            else:
                # 在绘图前进行样点边界包含性统计
                try:
                    if boundary_gdf is not None and len(boundary_gdf) > 0:
                        boundary_union = boundary_gdf.unary_union
                        total_points = len(results_with_coords)
                        inside_count = 0
                        outside_points: List[Dict[str, Any]] = []
                        for r in results_with_coords:
                            pt = Point(float(r["longitude"]), float(r["latitude"]))
                            if boundary_union.contains(pt) or boundary_union.touches(pt):
                                inside_count += 1
                            else:
                                outside_points.append({
                                    "farmland_id": r.get("farmland_id"),
                                    "sample_id": r.get("sample_id"),
                                    "longitude": r.get("longitude"),
                                    "latitude": r.get("latitude"),
                                    "flux_value": r.get("flux_value")
                                })

                        outside_count = total_points - inside_count
                        inside_pct = (inside_count / total_points * 100.0) if total_points > 0 else 0.0

                        self.logger.info(
                            f"样点边界检查 - 总数: {total_points}, 边界内: {inside_count} ({inside_pct:.2f}%), 边界外: {outside_count}")
                        if outside_count > 0:
                            sample_preview = outside_points[:10]
                            self.logger.warning(
                                f"存在 {outside_count} 个样点位于边界之外，绘图时将被掩膜隐藏。示例(最多10条): {sample_preview}")

                        # 在日志中打印边界检查统计结果
                        self.logger.info(
                            f"边界检查统计 - 地区: {area}, 层级: {level}, 计算类型: {calculation_type}, "
                            f"总样点: {total_points}, 边界内: {inside_count} ({inside_pct:.2f}%), "
                            f"边界外: {outside_count}"
                        )
                        if outside_count > 0 and len(outside_points) > 0:
                            sample_preview = outside_points[:5]  # 只显示前5个样本
                            self.logger.info(f"边界外样点示例(前5个): {sample_preview}")
                    else:
                        generated_files = {}
                except Exception as check_err:
                    self.logger.warning(f"样点边界包含性检查失败: {str(check_err)}")
                    # 保持已有 generated_files，不覆盖
            
            # 创建shapefile
            shapefile_path = csv_path.replace('.csv', '_points.shp')
            mapper.csv_to_shapefile(csv_path, shapefile_path, 
                                  lon_col='longitude', lat_col='latitude', value_col='flux_value')
            
            # 合并已生成文件映射
            generated_files.update({"csv": csv_path, "shapefile": shapefile_path})
            
            # 如果有模板栅格文件，创建栅格地图
            if template_raster:
                try:
                    # 创建栅格
                    raster_path = csv_path.replace('.csv', '_raster.tif')
                    raster_path, stats = mapper.vector_to_raster(
                        shapefile_path, template_raster, raster_path, 'flux_value',
                        resolution_factor=resolution_factor, boundary_shp=boundary_shp, boundary_gdf=boundary_gdf,
                        interpolation_method='nearest', enable_interpolation=enable_interpolation
                    )
                    generated_files["raster"] = raster_path
                    
                    # 创建栅格地图 - 使用英文标题避免中文乱码
                    title_mapping = {
                        "grain_removal": "Grain Removal Cd Flux",
                        "straw_removal": "Straw Removal Cd Flux"
                    }
                    map_title = title_mapping.get(calculation_type, "Cd Flux Removal")
                    
                    map_file = mapper.create_raster_map(
                        boundary_shp if boundary_shp else None,
                        raster_path,
                        map_output,
                        colormap=colormap,
                        title=map_title,
                        output_size=12,
                        dpi=300,
                        resolution_factor=4.0,
                        enable_interpolation=False,
                        interpolation_method='nearest',
                        boundary_gdf=boundary_gdf
                    )
                    generated_files["map"] = map_file
                    
                    # 创建直方图 - 使用英文标题避免中文乱码
                    histogram_title_mapping = {
                        "grain_removal": "Grain Removal Cd Flux Distribution",
                        "straw_removal": "Straw Removal Cd Flux Distribution"
                    }
                    histogram_title = histogram_title_mapping.get(calculation_type, "Cd Flux Distribution")
                    
                    histogram_file = mapper.create_histogram(
                        raster_path,
                        f"{histogram_output}.jpg",
                        title=histogram_title,
                        xlabel='Cd Flux (g/ha/a)',
                        ylabel='Frequency Density'
                    )
                    generated_files["histogram"] = histogram_file
                    
                except Exception as viz_error:
                    self.logger.warning(f"栅格可视化创建失败: {str(viz_error)}")
                    # 即使栅格可视化失败，也返回已生成的文件
            
            # 清理中间文件（默认开启，仅保留最终可视化）
            if cleanup_intermediate:
                try:
                    # 由于不再创建临时边界文件，所以传递None
                    self._cleanup_intermediate_files(generated_files, None)
                except Exception as cleanup_err:
                    self.logger.warning(f"中间文件清理失败: {str(cleanup_err)}")

            self.logger.info(f"✓ 成功创建 {calculation_type} 可视化，生成文件: {list(generated_files.keys())}")
            return generated_files
            
        except Exception as e:
            self.logger.error(f"创建可视化失败: {str(e)}")
            raise

    def _cleanup_intermediate_files(self, generated_files: Dict[str, str], boundary_shp: Optional[str]) -> None:
        """
        清理中间文件：CSV、Shapefile 及其配套文件、栅格TIFF；若边界为临时目录，则一并删除
        """
        import shutil
        import tempfile

        def _safe_remove(path: str) -> None:
            try:
                if path and os.path.exists(path) and os.path.isfile(path):
                    os.remove(path)
            except Exception:
                pass

        # 删除 CSV
        _safe_remove(generated_files.get("csv"))

        # 删除栅格
        _safe_remove(generated_files.get("raster"))

        # 删除 Shapefile 全家桶
        shp_path = generated_files.get("shapefile")
        if shp_path:
            base, _ = os.path.splitext(shp_path)
            for ext in (".shp", ".shx", ".dbf", ".prj", ".cpg"):
                _safe_remove(base + ext)

        # 如果边界文件来自系统临时目录，删除其所在目录
        if boundary_shp:
            temp_root = tempfile.gettempdir()
            try:
                if os.path.commonprefix([os.path.abspath(boundary_shp), temp_root]) == temp_root:
                    temp_dir = os.path.dirname(os.path.abspath(boundary_shp))
                    if os.path.isdir(temp_dir):
                        shutil.rmtree(temp_dir, ignore_errors=True)
            except Exception:
                pass
    
    def _get_boundary_file_for_area(self, area: str, level: str) -> Optional[str]:
        """
        为指定地区获取边界文件
        
        @param area: 地区名称
        @returns: 边界文件路径或None
        """
        try:
            # 仅从数据库严格获取边界（按指定层级精确匹配）
            boundary_path = self._create_boundary_from_database(area, level)
            if boundary_path:
                return boundary_path
            
            # 如果都没有找到，记录警告但不使用默认文件
            self.logger.warning(f"未找到地区 '{area}' 的专用边界文件，也无法从数据库获取")
            return None
            
        except Exception as e:
            self.logger.error(f"获取边界文件失败: {str(e)}")
            return None
    
    def _get_boundary_gdf_from_database(self, area: str, level: str) -> Optional[gpd.GeoDataFrame]:
        """
        直接从数据库获取边界数据作为GeoDataFrame
        
        @param area: 地区名称
        @param level: 行政层级
        @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:
                    self.logger.info(f"从数据库获取边界数据: {norm_area} ({level})")
                return boundary_gdf
                    
        except Exception as e:
            self.logger.warning(f"从数据库获取边界数据失败: {str(e)}")
            
        return None

    def _create_boundary_from_database(self, area: str, level: str) -> Optional[str]:
        """
        从数据库获取边界数据并创建临时shapefile
        
        @param area: 地区名称
        @returns: 临时边界文件路径或None
        @deprecated: 建议使用 _get_boundary_gdf_from_database 方法直接获取 GeoDataFrame
        """
        try:
            with SessionLocal() as db:
                norm_area = area.strip()
                boundary_geojson = get_boundary_geojson_by_name(db, norm_area, level=level)
                if boundary_geojson:
                    geometry_obj = shape(boundary_geojson["geometry"])
                    gdf = gpd.GeoDataFrame([boundary_geojson["properties"]], geometry=[geometry_obj], crs="EPSG:4326")
                    temp_dir = tempfile.mkdtemp()
                    boundary_path = os.path.join(temp_dir, f"{norm_area}_boundary.shp")
                    gdf.to_file(boundary_path, driver="ESRI Shapefile")
                    self.logger.info(f"从数据库创建边界文件: {boundary_path}")
                    return boundary_path
                    
        except Exception as e:
            self.logger.warning(f"从数据库创建边界文件失败: {str(e)}")
            
        return None


    def create_surface_runoff_visualization(self, area: str, level: str,
                                            output_dir: str = "app/static/surface_runoff",
                                            template_raster: str = "app/static/cd_flux/meanTemp.tif",
                                            boundary_shp: str = None,
                                            colormap: str = "green_yellow_red_purple",
                                            resolution_factor: float = 4.0,
                                            enable_interpolation: bool = False,
                                            cleanup_intermediate: bool = True) -> Dict[str, str]:
        """
        创建地表径流Cd通量可视化图表

        @param area: 地区名称
        @param level: 行政层级
        @param output_dir: 输出目录
        @param template_raster: 模板栅格文件路径
        @param boundary_shp: 边界shapefile路径
        @param colormap: 色彩方案
        @param resolution_factor: 分辨率因子
        @param enable_interpolation: 是否启用空间插值
        @param cleanup_intermediate: 是否清理中间文件
        @returns: 生成的图片文件路径字典
        """
        try:
            # 从数据库获取地表径流数据
            results_with_coords = []
            with SessionLocal() as db:
                # 查询地表径流数据（DB_Cd）及对应坐标
                runoff_data = db.query(
                    FarmlandData.farmland_id,
                    FarmlandData.sample_id,
                    FarmlandData.lon,
                    FarmlandData.lan,
                    FluxCdInputData.surface_runoff
                ).join(
                    FluxCdInputData,
                    (FluxCdInputData.farmland_id == FarmlandData.farmland_id) &
                    (FluxCdInputData.sample_id == FarmlandData.sample_id)
                ).all()

                if not runoff_data:
                    self.logger.warning("未找到地表径流数据")
                    return {}

                for data in runoff_data:
                    results_with_coords.append({
                        "farmland_id": data.farmland_id,
                        "sample_id": data.sample_id,
                        "longitude": data.lon,
                        "latitude": data.lan,
                        "flux_value": data.surface_runoff
                    })

            # 确保输出目录存在
            os.makedirs(output_dir, exist_ok=True)
            generated_files: Dict[str, str] = {}

            # 生成时间戳
            timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")

            # 创建CSV文件用于绘图
            csv_filename = f"surface_runoff_{area}_temp_{timestamp}.csv"
            csv_path = os.path.join(output_dir, csv_filename)

            # 准备绘图数据
            plot_data = []
            for result in results_with_coords:
                plot_data.append({
                    "longitude": result["longitude"],
                    "latitude": result["latitude"],
                    "flux_value": result["flux_value"]
                })

            # 保存为CSV
            df = pd.DataFrame(plot_data)
            df.to_csv(csv_path, index=False, encoding='utf-8-sig')

            # 初始化绘图工具
            mapper = MappingUtils()

            # 生成输出文件路径
            map_output = os.path.join(output_dir, f"surface_runoff_{area}_map_{timestamp}")
            histogram_output = os.path.join(output_dir, f"surface_runoff_{area}_histogram_{timestamp}")

            # 检查模板文件是否存在
            if not os.path.exists(template_raster):
                self.logger.warning(f"模板栅格文件不存在: {template_raster}")
                template_raster = None

            # 动态获取边界数据（严格使用指定层级）
            if not level:
                raise ValueError("必须提供行政层级 level：county | city | province")

            # 直接从数据库获取边界GeoDataFrame
            boundary_gdf = self._get_boundary_gdf_from_database(area, level)
            boundary_shp = None  # 不再需要临时边界文件

            if boundary_gdf is None:
                self.logger.warning(f"未找到地区 '{area}' 的边界数据，将不使用边界裁剪")
            else:
                # 在绘图前进行样点边界包含性统计
                try:
                    if boundary_gdf is not None and len(boundary_gdf) > 0:
                        boundary_union = boundary_gdf.unary_union
                        total_points = len(results_with_coords)
                        inside_count = 0
                        outside_points: List[Dict[str, Any]] = []
                        for r in results_with_coords:
                            pt = Point(float(r["longitude"]), float(r["latitude"]))
                            if boundary_union.contains(pt) or boundary_union.touches(pt):
                                inside_count += 1
                            else:
                                outside_points.append({
                                    "farmland_id": r.get("farmland_id"),
                                    "sample_id": r.get("sample_id"),
                                    "longitude": r.get("longitude"),
                                    "latitude": r.get("latitude"),
                                    "flux_value": r.get("flux_value")
                                })

                        outside_count = total_points - inside_count
                        inside_pct = (inside_count / total_points * 100.0) if total_points > 0 else 0.0

                        self.logger.info(
                            f"样点边界检查 - 总数: {total_points}, 边界内: {inside_count} ({inside_pct:.2f}%), 边界外: {outside_count}")
                        if outside_count > 0:
                            sample_preview = outside_points[:10]
                            self.logger.warning(
                                f"存在 {outside_count} 个样点位于边界之外，绘图时将被掩膜隐藏。示例(最多10条): {sample_preview}")
                except Exception as check_err:
                    self.logger.warning(f"样点边界包含性检查失败: {str(check_err)}")

            # 创建shapefile
            shapefile_path = csv_path.replace('.csv', '_points.shp')
            mapper.csv_to_shapefile(csv_path, shapefile_path,
                                    lon_col='longitude', lat_col='latitude', value_col='flux_value')

            # 合并已生成文件映射
            generated_files.update({"csv": csv_path, "shapefile": shapefile_path})

            # 如果有模板栅格文件，创建栅格地图
            if template_raster:
                try:
                    # 创建栅格
                    raster_path = csv_path.replace('.csv', '_raster.tif')
                    raster_path, stats = mapper.vector_to_raster(
                        shapefile_path, template_raster, raster_path, 'flux_value',
                        resolution_factor=resolution_factor, boundary_shp=boundary_shp, boundary_gdf=boundary_gdf,
                        interpolation_method='nearest', enable_interpolation=enable_interpolation
                    )
                    generated_files["raster"] = raster_path

                    # 创建栅格地图 - 使用英文标题避免中文乱码
                    map_title = "Surface Runoff Cd Flux"

                    map_file = mapper.create_raster_map(
                        boundary_shp if boundary_shp else None,
                        raster_path,
                        map_output,
                        colormap=colormap,
                        title=map_title,
                        output_size=12,
                        dpi=300,
                        resolution_factor=4.0,
                        enable_interpolation=False,
                        interpolation_method='nearest',
                        boundary_gdf=boundary_gdf
                    )
                    generated_files["map"] = map_file

                    # 创建直方图 - 使用英文标题避免中文乱码
                    histogram_title = "Surface Runoff Cd Flux Distribution"

                    histogram_file = mapper.create_histogram(
                        raster_path,
                        f"{histogram_output}.jpg",
                        title=histogram_title,
                        xlabel='Cd Flux (g/ha/a)',
                        ylabel='Frequency Density'
                    )
                    generated_files["histogram"] = histogram_file

                except Exception as viz_error:
                    self.logger.warning(f"栅格可视化创建失败: {str(viz_error)}")
                    # 即使栅格可视化失败，也返回已生成的文件

            # 清理中间文件（默认开启，仅保留最终可视化）
            if cleanup_intermediate:
                try:
                    self._cleanup_intermediate_files(generated_files, None)
                except Exception as cleanup_err:
                    self.logger.warning(f"中间文件清理失败: {str(cleanup_err)}")

            self.logger.info(f"✓ 成功创建地表径流Cd通量可视化，生成文件: {list(generated_files.keys())}")
            return generated_files

        except Exception as e:
            self.logger.error(f"创建地表径流可视化失败: {str(e)}")
            raise


    def get_groundwater_leaching_data(self, area: str) -> List[Dict[str, Any]]:
        """
        获取地下渗流（地下水渗漏）通量数据及坐标

        @param area: 地区名称
        @returns: 包含坐标和通量值的结果列表
        """
        try:
            with SessionLocal() as db:
                # 查询地下渗流数据及对应坐标
                leaching_data = db.query(
                    FluxCdInputData.groundwater_leaching,
                    FarmlandData.farmland_id,
                    FarmlandData.sample_id,
                    FarmlandData.lon,
                    FarmlandData.lan
                ).join(
                    FarmlandData,
                    and_(
                        FluxCdInputData.farmland_id == FarmlandData.farmland_id,
                        FluxCdInputData.sample_id == FarmlandData.sample_id
                    )
                ).all()

                if not leaching_data:
                    self.logger.warning(f"未找到地区 '{area}' 的地下渗流数据")
                    return []

                results = []
                for data in leaching_data:
                    results.append({
                        "farmland_id": data.farmland_id,
                        "sample_id": data.sample_id,
                        "longitude": data.lon,
                        "latitude": data.lan,
                        "flux_value": data.groundwater_leaching
                    })

                self.logger.info(f"成功获取 {len(results)} 个样点的地下渗流数据")
                return results

        except Exception as e:
            self.logger.error(f"获取地下渗流数据失败: {str(e)}")
            raise


    def create_groundwater_leaching_visualization(self, area: str, level: str,
                                                  output_dir: str = "app/static/groundwater_leaching",
                                                  template_raster: str = "app/static/cd_flux/meanTemp.tif",
                                                  colormap: str = "blues",
                                                  resolution_factor: float = 4.0,
                                                  enable_interpolation: bool = False,
                                                  cleanup_intermediate: bool = True) -> Dict[str, str]:
        """
        创建地下渗流（地下水渗漏）通量可视化图表

        @param area: 地区名称
        @param level: 行政层级
        @param output_dir: 输出目录
        @param template_raster: 模板栅格文件路径
        @param colormap: 色彩方案
        @param resolution_factor: 分辨率因子
        @param enable_interpolation: 是否启用空间插值
        @param cleanup_intermediate: 是否清理中间文件
        @returns: 生成的图片文件路径字典
        """
        try:
            # 获取地下渗流数据
            leaching_results = self.get_groundwater_leaching_data(area)

            if not leaching_results:
                return {
                    "success": False,
                    "message": f"没有找到地区 '{area}' 的地下渗流数据",
                    "data": None
                }

            # 确保输出目录存在
            os.makedirs(output_dir, exist_ok=True)
            generated_files: Dict[str, str] = {}

            # 生成时间戳
            timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")

            # 创建CSV文件用于绘图
            csv_filename = f"groundwater_leaching_{area}_temp_{timestamp}.csv"
            csv_path = os.path.join(output_dir, csv_filename)

            # 准备绘图数据
            plot_data = []
            for result in leaching_results:
                plot_data.append({
                    "longitude": result["longitude"],
                    "latitude": result["latitude"],
                    "flux_value": result["flux_value"]
                })

            # 保存为CSV
            df = pd.DataFrame(plot_data)
            df.to_csv(csv_path, index=False, encoding='utf-8-sig')

            # 初始化绘图工具
            mapper = MappingUtils()

            # 生成输出文件路径
            map_output = os.path.join(output_dir, f"groundwater_leaching_{area}_map_{timestamp}")
            histogram_output = os.path.join(output_dir, f"groundwater_leaching_{area}_histogram_{timestamp}")

            # 检查模板文件是否存在
            if not os.path.exists(template_raster):
                self.logger.warning(f"模板栅格文件不存在: {template_raster}")
                template_raster = None

            # 直接从数据库获取边界GeoDataFrame
            boundary_gdf = self._get_boundary_gdf_from_database(area, level)
            boundary_shp = None  # 不再需要临时边界文件

            # 在绘图前进行样点边界包含性统计
            if boundary_gdf is not None and len(boundary_gdf) > 0:
                boundary_union = boundary_gdf.unary_union
                total_points = len(leaching_results)
                inside_count = 0
                for r in leaching_results:
                    pt = Point(float(r["longitude"]), float(r["latitude"]))
                    if boundary_union.contains(pt) or boundary_union.touches(pt):
                        inside_count += 1

                outside_count = total_points - inside_count
                inside_pct = (inside_count / total_points * 100.0) if total_points > 0 else 0.0

                self.logger.info(
                    f"样点边界检查 - 总数: {total_points}, 边界内: {inside_count} ({inside_pct:.2f}%), 边界外: {outside_count}")

            # 创建shapefile
            shapefile_path = csv_path.replace('.csv', '_points.shp')
            mapper.csv_to_shapefile(csv_path, shapefile_path,
                                    lon_col='longitude', lat_col='latitude', value_col='flux_value')

            # 合并已生成文件映射
            generated_files.update({"csv": csv_path, "shapefile": shapefile_path})

            # 如果有模板栅格文件，创建栅格地图
            if template_raster:
                try:
                    # 创建栅格
                    raster_path = csv_path.replace('.csv', '_raster.tif')
                    raster_path, stats = mapper.vector_to_raster(
                        shapefile_path, template_raster, raster_path, 'flux_value',
                        resolution_factor=resolution_factor, boundary_gdf=boundary_gdf,
                        interpolation_method='nearest', enable_interpolation=enable_interpolation
                    )
                    generated_files["raster"] = raster_path

                    # 创建栅格地图 - 使用特定标题
                    map_title = "Groundwater Leaching Cd Flux"

                    map_file = mapper.create_raster_map(
                        boundary_shp,
                        raster_path,
                        map_output,
                        colormap=colormap,
                        title=map_title,
                        output_size=12,
                        dpi=300,
                        resolution_factor=4.0,
                        enable_interpolation=False,
                        interpolation_method='nearest',
                        boundary_gdf=boundary_gdf
                    )
                    generated_files["map"] = map_file

                    # 创建直方图 - 使用特定标题
                    histogram_title = "Groundwater Leaching Cd Flux Distribution"

                    histogram_file = mapper.create_histogram(
                        raster_path,
                        f"{histogram_output}.jpg",
                        title=histogram_title,
                        xlabel='Cd Flux (g/ha/a)',
                        ylabel='Frequency Density'
                    )
                    generated_files["histogram"] = histogram_file

                except Exception as viz_error:
                    self.logger.warning(f"栅格可视化创建失败: {str(viz_error)}")

            # 清理中间文件
            if cleanup_intermediate:
                try:
                    self._cleanup_intermediate_files(generated_files, None)
                except Exception as cleanup_err:
                    self.logger.warning(f"中间文件清理失败: {str(cleanup_err)}")

            self.logger.info(f"✓ 成功创建地下渗流通量可视化，生成文件: {list(generated_files.keys())}")
            return generated_files

        except Exception as e:
            self.logger.error(f"创建地下渗流可视化失败: {str(e)}")
            return {
                "success": False,
                "message": f"可视化创建失败: {str(e)}",
                "data": None
            }