AcidMap/app/services/cd_flux_removal_service.py

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

import logging
import math
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.parameters import Parameters
from ..models.CropCd_output import CropCdOutputData


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) -> Dict[str, Any]:
        """
        根据地区计算籽粒移除Cd通量
        
        @param area: 地区名称
        @returns: 计算结果字典
        
        计算公式：籽粒移除(g/ha/a) = EXP(LnCropCd) * F11 * 0.5 * 15 / 1000
        """
        try:
            with SessionLocal() as db:
                # 查询指定地区的参数
                parameter = db.query(Parameters).filter(Parameters.area == area).first()
                
                if not parameter:
                    return {
                        "success": False,
                        "message": f"未找到地区 '{area}' 的参数数据",
                        "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) -> 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:
                # 查询指定地区的参数
                parameter = db.query(Parameters).filter(Parameters.area == area).first()
                
                if not parameter:
                    return {
                        "success": False,
                        "message": f"未找到地区 '{area}' 的参数数据",
                        "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
            }