Zdrojový kód pro modules.SARCA_lib

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#  Project: RadAgro
#  File: SARCA_lib.py
#
#  Author: Dr. Jakub Brom
#
#  Copyright (c) 2020. Dr. Jakub Brom, University of South Bohemia in
#  České Budějovice, Faculty of Agriculture.
#
#  This program is free software: you can redistribute it and/or modify
#      it under the terms of the GNU General Public License as published by
#      the Free Software Foundation, either version 3 of the License, or
#      (at your option) any later version.
#
#      This program is distributed in the hope that it will be useful,
#      but WITHOUT ANY WARRANTY; without even the implied warranty of
#      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#      GNU General Public License for more details.
#
#      You should have received a copy of the GNU General Public License
#      along with this program.  If not, see <https://www.gnu.org/licenses/>
#
#  Last changes: 26.11.20 16:58
#
#  Begin: 2020/10/31
#
#  Description: Calculation of crop growth features and radioactivity
#  contamination features


import numpy as np


[dokumentace]class SARCALib: """Library for calculation of crops growth parameters and radioactivity contamination of crops""" def __init__(self): return
[dokumentace] def calculateGrowthCoefs(self, dw_max, dw_min=0.1): """ Calculate default values of growth curve parameters - slope \ (m) and intercept (n). :param dw_max: Maximal dry mass of particular crop \ :math:`(t.ha^{-1})` :param dw_min: Minimal dry mass of particular crop \ :math:`(t.ha^{-1})`. Default value is 0.1 \ :math:`(t.ha^{-1})`. :return: Slope of growth curve :return: Intercept of growth curve """ intercept = np.log(abs(dw_max/(dw_min-0.001) - 1)) const = -np.log(abs(dw_max/(dw_max-0.001) - 1)) slope = intercept + const return slope, intercept
[dokumentace] def timeSeriesConst(self, t_accident, t0, t1): """ Calculation of relative position of the radioactive accident in the crop growth time series [0, 1]. :param t_accident: Date of radioactive accident as number of \ day in year :param t0: Date of sowing crop as number of day in year :param t1: Date of harvesting crop as number of day in year :return: Relative position of the radioactive accident in the \ crop growth time series """ ignore_zero = np.seterr(all="ignore") t0 = np.where(t0 >= t1, 92, t0) ts = np.where(t0 >= t1, 0, (t_accident - t0) / (t1 - t0)) ts[ts > 1.0] = 1.0 return ts
[dokumentace] def timeScale(self, ts, coef_m, coef_n): """ Calculates scaling value t_scale of crop growth time series. :param ts: Relative position of the radioactive accident in \ the crop growth time series :param coef_m: Scaling parameter :param coef_n: Scaling parameter :return: Scaling value of growing time series """ t_scale = coef_m * ts - coef_n return t_scale
[dokumentace] def dryMass(self, DW_max, t_accident, t0, t1, coef_m, coef_n): """ Calculation of actual biomass of the particular crop. :param DW_max: Maximal dry mass of the crops above ground \ biomass in the growing season :math:`(t.ha^{-1})` :param t_accident: Date of radioactive accident as number of \ day in year :param t0: Date of sowing crop as number of day in year :param t1: Date of harvesting crop as number of day in year :param coef_m: Scaling parameter :param coef_n: Scaling parameter :return: Dry mass of above ground biomass :math:`(t.ha^{-1})` """ # Time scaling ts = self.timeSeriesConst(t_accident, t0, t1) t_scale = self.timeScale(ts, coef_m, coef_n) # Dry mas of above ground biomass calculation dry_mass = np.where(ts != 0.0, DW_max / (1.0 + np.exp(-t_scale)), 0.0) return dry_mass
[dokumentace] def leafAreaIndex(self, dw_act, DW_max, LAI_max, R_min, t_accident, t0, t1): """ Calculation of actual LAI of the particular crop. :param dw_act: actual dry mass of the crops above ground \ biomass :math:`(t.ha^{-1})` :param DW_max: Maximal dry mass of the crops above ground \ biomass in the growing season :math:`(t.ha^{-1})` :param LAI_max: Maximal LAI of the crops in the growing season :param R_min: Minimal relative wight moisture of the biomass (%) :param t_accident: Date of radioactive accident as number of \ day in year :param t0: Date of sowing crop as number of day in year :param t1: Date of harvesting crop as number of day in year :return: Actual LAI of the crop. """ # Time scaling ts = self.timeSeriesConst(t_accident, t0, t1) # Dry mas of above ground biomass calculation LAI = np.where(ts <= 0.7, LAI_max ** 2 * dw_act / DW_max, (-3.6511 * LAI_max + 0.19993 * R_min - 6.66309) * ts ** 2 + (3.9841 * LAI_max - 0.14 * R_min + 4.6668) * ts) LAI = np.where(LAI > LAI_max, LAI_max, LAI) return LAI
[dokumentace] def interceptFactor(self, LAI, precip, DW, k=1.0, S=0.2): """ Interception factor for both dry and wet deposition of radionuclides. :param LAI: Leaf Area Index (unitless) :param k: Constant of radionuclide: I = 0.5, Sr and Ba = 2.0, \ Cs and another radionuclides = 1.0 :param precip: Precipitation amount (mm) for period of \ deposition (ca 24 hours after radiation accident). :param DW: Amount of fresh biomass :math:`(t.ha^{-1})` :param S: Mean thickness of water film at plant leaves (mm). \ Default S = 0.2 mm :return: Interception Factor (rel.) """ try: IF = LAI * k * S * (1.0 - np.exp((-np.log(2.0)) / (3.0 * S) * ( precip + 0.0001))) / (precip + 0.0001) IF = np.array(IF) IF[IF > 0.8] = 0.8 IF[DW < 0.1] = 0.0 except ArithmeticError: raise ArithmeticError("Interception factor has not been " "calculated") return IF
[dokumentace] def contBiomass(self, depo, IF): """ Radiaoctive contamination of biomass :math:`(Bq.m^{-2})` :param depo: Total radioactive deposition :math:`(Bq.m^{-2})` :param IF: Interception Factor (rel.) :return: Radioactive contamination of biomass :math:`(Bq.m^{-2})` """ try: cont_biomass = depo * IF cont_biomass = np.where(cont_biomass < 0.0, 0.0, cont_biomass) except ArithmeticError: raise ArithmeticError("Vegetation biomass radioactive " "contamination has not been " "calculated") return cont_biomass
[dokumentace] def contSoil(self, depo, IF): """ Radiaoctive contamination of soil :math:`(Bq.m^{-2})` :param depo: Total radioactive deposition :math:`(Bq.m^{-2})` :param IF: Interception Factor (rel.) :return: Radioactive contamination of soil :math:`(Bq.m^{-2})` """ try: cont_soil = depo * (1 - IF) cont_soil = np.where(cont_soil < 0.0, 0.0, cont_soil) except ArithmeticError: raise ArithmeticError("Soil radioactive contamination" " has not been calculated") return cont_soil
[dokumentace] def contMass(self, cont_biomass, fresh_biomass): """ Calculation radioactive contamination of fresh vegetation mass \ :math:`(Bq.kg^{-1})` :param cont_biomass: Radioactive deposition on biomass \ :math:`(Bq.m^{-2})` :param fresh_biomass: Amount of fresh biomass of vegetation \ :math:`(t.ha^{-1})` :return: Fresh vegetation mass radioactive contamination \ :math:`(Bq.kg^{-1})` """ ignore_zero = np.seterr(all="ignore") try: cont_mass = cont_biomass / (fresh_biomass * 0.1) cont_mass[cont_mass < 0.0] = 0.0 except ArithmeticError: raise ArithmeticError("Mass Radioactive " "contamination of fresh biomass has" " not been calculated") return cont_mass
[dokumentace] def referLevel(self, depo, ref_level1=5000, ref_level2=3000000): """ Mask of radioactive deposition for three reference levels (categories): \n *0: Low* \n *1: Middle* \n *2: High* \n :param depo: Total radioactive deposition :math:`(Bq.m^{-2})` :param ref_level1: Lower reference level treshold \ :math:`(Bq.m^{-2})`. Default value = 5000 :math:`Bq.m^{-2}` :param ref_level2: Upper reference level treshold \ :math:`(Bq.m^{-2})`. Default value = 3000000 :math:`Bq.m^{-2}` :return: Mask of radioactive deposition for three reference levels """ try: reference_groups = np.where(depo >= ref_level2, 2.0, 1.0) reference_groups = np.where(depo >= ref_level1, reference_groups, 0.0) except ArithmeticError: raise ArithmeticError("Mask for reference levels has not " "been calculated") return reference_groups