1import os
   2import smash
   3
   4# import pyhdf5_handler
   5# from pyhdf5_handler.src import hdf5_handler
   6import numpy as np
   7import pandas as pd
   8import datetime
   9
  10# from smashbox.src import param
  11from smashbox.model import setup
  12from smashbox.model import mesh
  13from smashbox.plot import myplot
  14from smashbox.tools import geo_toolbox
  15from smashbox.tools import tools
  16from smashbox.model import atmos_data_connector
  17from smashbox.stats import mystats
  18from smashbox.read_inputdata import smashmodel
  19from smashbox.model import smash_model
  20
  21import pyhdf5_handler
  22import copy
  23import pickle
  24
  25
  26class model:
  27    """Main class model() which has a complete set of functions, class and attributes to
  28    build, run and manipulate the Smash model, compute statistical criterium and
  29    plot graphics."""
  30
  31    def __init__(self, name, myparam):
  32
  33        self._model_name = name
  34        """The attribute _model_name is contain the name of the current object"""
  35
  36        self._myparam = copy.deepcopy(myparam)
  37        """The attribute _myparam is a copy of the parent class smashbox.myparam.
  38        This class stores the main parameters used for building the Smash model"""
  39
  40        self.mysetup = setup.setup(self._myparam.param)
  41        """The attribute mysetup own the class setup.setup(). This class stores the Smash
  42        setup for the hydrological simulation and some helpers to manipulate these
  43        parameters."""
  44
  45        self.mymesh = mesh.mesh(self.mysetup)
  46        """The attribute mymesh own the class mesh.mesh(). This class stores the Smash
  47        mesh used for the hydrological simulation and some helpers to manipulate this mesh.
  48        """
  49
  50        self._fstates = None
  51        self._istates = None
  52        self._myatmos_data_connector = None
  53
  54        self.mysmashmodel = smash_model.smash_model()
  55        """The attribute mysmashmodel stores the Smash model object created with
  56        attributes mysetup and -mymesh."""
  57
  58        self.warmup_model = smash_model.smash_model()
  59        """The attribute warmup_model store a smash model used for warmup and compatible 
  60        with the model in attribute mysmashmodel"""
  61
  62        self.optimize_model = smash_model.smash_model()
  63        """The attribute optimize_model store a smash model used for optimization and compatible 
  64        with the model in attribute mysmashmodel"""
  65
  66        self.validation_model = smash_model.smash_model()
  67        """The attribute validation_model stores the Smash model object used for validation"""
  68
  69        self.myplot = myplot.myplot(self)
  70        """The attribute myplot own the class myplot.myplot(). This class stores plotting
  71        # capabilities on the smash model object."""
  72
  73    def generate_mesh(
  74        self,
  75        max_depth: float = 1.0,
  76        query: str | None = None,
  77        area_error_th: None | float = None,
  78        lacuna_threshold: None | float = None,
  79    ):
  80        """
  81        Generate the mesh of the Smash model
  82
  83        Parameters
  84        ----------
  85
  86        max_depth : `int`, default 1
  87            The maximum depth accepted by the algorithm to find the catchment outlet.
  88            A **max_depth** of 1 means that the algorithm will search among the
  89            combinations in
  90            (``row - 1``, ``row``, ``row + 1``; ``col - 1``, ``col``, ``col + 1``),
  91            the coordinates that minimize
  92            the relative error between the given catchment area and the modeled
  93            catchment area calculated from the
  94            flow directions file.
  95        :param query: Any pandas dataframe query as string: '(SURF>20) & (SURF<100)'. This query
  96        must be build using the field (column name) in the outlet database.
  97        https://pandas.pydata.org/docs/user_guide/indexing.html#the-query-method
  98        :type query: str
  99        area_error_th: float | None
 100            The tolerance error for the difference between the observed and simulated
 101             surface. The error is computed as follow:
 102                 Serror=abs(Ssim-Sobs)/Sobs
 103            All outlets where `Serror > area_error_th` will be automatically removed from
 104            the mesh.
 105        :type area_error_th: float
 106        :param lacuna_threshold: Lacuna threshold for the discharges in percent. All gauge where the proportion of lacuna exceed this threshold will be removed.
 107        :type: float | Nonetype
 108
 109        Examples
 110        --------
 111
 112        >>> es=smashbox.SmashBox()
 113        >>> sb.newmodel("RealCollobrier")
 114        >>> sb.RealCollobrier.generate_mesh(min_surf=5, max_surf=100)
 115
 116        """
 117        self.mymesh.generate_mesh(
 118            self._myparam.param,
 119            max_depth=max_depth,
 120            query=query,
 121            area_error_th=area_error_th,
 122            lacuna_threshold=lacuna_threshold,
 123        )
 124
 125    def model(
 126        self,
 127        setup: dict | None = None,
 128        mesh: dict | None = None,
 129        read_data=None,
 130    ):
 131        """
 132        Smash model object creation. This function wrap smash.Model(). Setup and mesh
 133        argument are optional since these dictionnary are hosted by the smashbox object.
 134
 135        Parameters
 136        ----------
 137        setup: dict | None
 138            The Smash setup (optionnal), if None the smashbox setup will be used
 139        mesh: dict | None
 140            The Smash mesh (optionnal), if None the smashbox mesh will be used
 141
 142        Examples
 143        --------
 144
 145        >>> es=smashbox.SmashBox()
 146        >>> sb.newmodel("RealCollobrier")
 147        >>> sb.RealCollobrier.generate_mesh(run=False)
 148        >>> sb.RealCollobrier.model()
 149
 150        """
 151        if setup is None:
 152            setup = self.mysetup.setup.copy()
 153
 154        if mesh is None:
 155            mesh = self.mymesh.mesh.copy()
 156
 157        if read_data is False:
 158            setup.update(
 159                {
 160                    "read_prcp": False,
 161                    "read_pet": False,
 162                    "read_snow": False,
 163                    "read_qobs": False,
 164                    "read_temp": False,
 165                }
 166            )
 167        if read_data is True:
 168            setup.update(
 169                {
 170                    "read_prcp": True,
 171                    "read_pet": True,
 172                    "read_snow": True,
 173                    "read_qobs": True,
 174                    "read_temp": True,
 175                }
 176            )
 177
 178        self.mysmashmodel.smash = self._model(setup=setup, mesh=mesh)
 179
 180        if read_data is True:
 181            self._gathering_atmosdata()
 182
 183        self._gathering_parameters(self.mysmashmodel.smash)
 184
 185    def _model(self, setup=None, mesh=None):
 186        """
 187        Smash model object creation. This function wrap smash.Model()
 188
 189        Parameters
 190        ----------
 191        setup: dict | None
 192            The Smash setup (optionnal), if None the smashbox setup will be used
 193        mesh: dict | None
 194            The Smash mesh (optionnal), if None the smashbox mesh will be used
 195
 196        Examples
 197        --------
 198
 199        >>> es=smashbox.SmashBox()
 200        >>> sb.newmodel("RealCollobrier")
 201        >>> sb.RealCollobrier.generate_mesh(run=False)
 202        >>> sb.RealCollobrier.model()
 203
 204        """
 205        if setup is None:
 206            print("</> The input setup is None ")
 207            return None
 208
 209        if mesh is None:
 210            print(
 211                "</> input mesh is None. use smashbox.model.model.generate_mesh()"
 212            )
 213            return None
 214
 215        if self._myparam.param.enhanced_smash_input_data:
 216            model = smashmodel.SmashModel(setup, mesh)
 217        else:
 218            model = smash.Model(setup, mesh)
 219
 220        return model
 221
 222    def _gathering_parameters(self, model):
 223
 224        if self.optimize_model.smash is not None:
 225            print(
 226                f"</> Getting parameters from previously optimized model ..."
 227            )
 228            model.rr_parameters = self.optimize_model.smash.rr_parameters
 229        else:
 230            print(
 231                f"</> Importing parameters from {self._myparam.param._smash_parameters} ..."
 232            )
 233            self.import_parameters(model=model)
 234            self._transform_parameters(
 235                model=model, dt_origin=self._myparam.param._smash_parameters_dt
 236            )
 237
 238    def _transform_parameters(
 239        self, model=None, dt_origin: None | float = None
 240    ):
 241        """
 242        Function to transform the parameters according the model time-step and the original
 243        time-step used for calibrated the parameters.
 244        :param dt_origin: Original time-step used for generate the calibrated parameter,
 245        defaults to None
 246        :type dt_origin: None | float, optional
 247
 248        """
 249        # if dt_origin is None:
 250        #     raise ValueError(
 251        #         " Argument dt_origin is None. This must be filled with the value "
 252        #         "of the timestep used to calibrate the parameters."
 253        #     )
 254
 255        # if not hasattr(self, "mysmashmodel.smash") or self.mysmashmodel.smash is None:
 256        #     raise ValueError(
 257        #         "</> mysmashmodel.smash attribute does not exist or is None. "
 258        #         "The model must be created first..."
 259        #     )
 260        if model is None:
 261            return
 262
 263        dt_target = model.setup.dt
 264
 265        if dt_origin is None:
 266            return
 267
 268        if dt_target == dt_origin:
 269            return
 270
 271        print(
 272            f"</> Tranforming parameters `ct`, `kexec`, `llr` calibrated with a time-step "
 273            f"of {dt_origin}s to the modeled time-step {dt_target}s."
 274        )
 275        parameters_tranfsorm = ["ct", "kexc", "llr"]
 276        power_tranform = [1.0 / 4.0, -1.0 / 8.0, 1.0]
 277
 278        for i, param in enumerate(parameters_tranfsorm):
 279            index_param = np.where(param == model.rr_parameters.keys)[0]
 280            if len(index_param) > 0:
 281                model.rr_parameters.values[:, :, index_param[0]] = (
 282                    model.rr_parameters.values[:, :, index_param[0]]
 283                    * (dt_origin / dt_target) ** power_tranform[i]
 284                )
 285
 286    def _gathering_atmosdata(self):
 287
 288        if self._myatmos_data_connector is not None:
 289
 290            print("</> Gathering atmos data ...")
 291
 292            if (
 293                self._myatmos_data_connector.input_ntimestep
 294                != self.mysmashmodel.smash.setup.ntime_step
 295            ):
 296                print(
 297                    "</> Warnings: Inconsistant ntime_step "
 298                    f"{self._myatmos_data_connector.input_ntimestep}"
 299                    f"!={self.mysmashmodel.setup.ntime_step}, "
 300                    "the Smash model will be rebuild."
 301                )
 302
 303                self._model()
 304                self._gathering_parameters(self.mysmashmodel.smash)
 305
 306            if self._myatmos_data_connector.smash_prcp is not None:
 307                self.mysmashmodel.smash.atmos_data.prcp = (
 308                    self._myatmos_data_connector.smash_prcp
 309                )
 310
 311            if self._myatmos_data_connector.smash_pet is not None:
 312                self.mysmashmodel.smash.atmos_data.pet = (
 313                    self._myatmos_data_connector.smash_pet
 314                )
 315
 316            if self.mysmashmodel.smash.setup.prcp_partitioning:
 317                print("</> Compute prcp partitionning ...")
 318                smash.fcore._mw_atmos_statistic.compute_prcp_partitioning(
 319                    model.setup, model.mesh, model._input_data
 320                )
 321
 322            if self.mysmashmodel.smash.setup.compute_mean_atmos:
 323                print("</> Computing mean atmospheric data")
 324                smash.fcore._mw_atmos_statistic.compute_mean_atmos(
 325                    self.mysmashmodel.smash.setup,
 326                    self.mysmashmodel.smash.mesh,
 327                    self.mysmashmodel.smash._input_data,
 328                )
 329
 330    @tools.autocast_args
 331    def model_warmup(self, warmup: int = 365):
 332        """
 333        Smash model warmup function. This function warm the curent model by creating a new
 334        on with attribute warmup_model.smash. The final states of warmup_model.smash are copied to the
 335        initial states of mysmashmodel.smash.
 336
 337        Parameters
 338        ----------
 339
 340        warmup: None | int
 341            a integer of the number of days used for warming the model.
 342
 343        """
 344        print("</> Warmup the Smash-model...")
 345
 346        if warmup is not None:
 347            try:
 348                timedelta = datetime.timedelta(days=warmup)
 349            except:
 350                raise ValueError(
 351                    f"warmup arg `{warmup}` is not a valid time delta."
 352                )
 353
 354            w_setup = copy.deepcopy(self.mysetup.setup)
 355
 356            w_setup["end_time"] = w_setup["start_time"]
 357            w_setup["start_time"] = datetime.datetime.strftime(
 358                datetime.datetime.fromisoformat(w_setup["start_time"])
 359                - timedelta,
 360                "%Y-%m-%d %H:%M",
 361            )
 362            w_setup["read_qobs"] = False
 363            w_setup["read_prcp"] = True
 364            w_setup["read_pet"] = True
 365
 366            if self.warmup_model.smash is None:
 367                self.warmup_model.smash = self._model(
 368                    setup=w_setup, mesh=self.mymesh.mesh
 369                )
 370
 371            self._gathering_parameters(self.warmup_model.smash)
 372            self.warmup_model.smash.forward_run()
 373
 374            # TODO FIX BUG IN SMASH when states > 1.0
 375            mask_sup = np.where(
 376                self.warmup_model.smash.rr_final_states.values > 1
 377            )
 378            self.warmup_model.smash.rr_final_states.values[mask_sup] = 0.9999
 379
 380            if (
 381                hasattr(self, "mysmashmodel")
 382                and self.mysmashmodel.smash is not None
 383            ):
 384
 385                self.mysmashmodel.smash.rr_initial_states = (
 386                    self.warmup_model.smash.rr_final_states
 387                )
 388
 389    @tools.autocast_args
 390    def validate(
 391        self,
 392        start_time=None,
 393        end_time=None,
 394        warmup: int = 365,
 395        cost_options: dict | None = None,
 396        common_options: dict | None = None,
 397        return_options: dict | None = None,
 398    ):
 399        """
 400        Smash model validation function. This function run the smash model on a other period
 401        where start_time and end_time may differ from the setup. Currently, it won't work with the
 402        atmos_data_connector wich only apply to the main model 'mysmashmodel'.
 403
 404        Parameters
 405        ----------
 406
 407        :param start_time: The start time of the optimization, format "YYYY-mm-dd HH:MM", defaults to None
 408        :type start_time: None | str, optional
 409        :param end_time: The end time of the optimization, format "YYYY-mm-dd HH:MM", defaults to None
 410        :type end_time: None | str, optional
 411        :param warmup: a integer of the number of days used for warming the model.
 412        :type warmup: int, optional, default 365
 413        """
 414        print("</> Validation of the Smash-model...")
 415
 416        v_setup = copy.deepcopy(self.mysetup.setup)
 417
 418        v_setup["start_time"] = start_time
 419        v_setup["end_time"] = end_time
 420        v_setup["read_qobs"] = True
 421        v_setup["read_prcp"] = True
 422        v_setup["read_pet"] = True
 423
 424        if self.validation_model.smash is None:
 425            self.validation_model.smash = self._model(
 426                setup=v_setup, mesh=self.mymesh.mesh
 427            )
 428            print(
 429                f"</> Getting parameters from the validation model `validation_model.smash` ..."
 430            )
 431
 432            self._gathering_parameters(self.validation_model.smash)
 433
 434        if warmup > 0:
 435
 436            try:
 437                timedelta = datetime.timedelta(days=warmup)
 438            except:
 439                raise ValueError(
 440                    f"warmup arg `{warmup}` is not a valid time delta."
 441                )
 442
 443            v_setup["end_time"] = v_setup["start_time"]
 444            v_setup["start_time"] = datetime.datetime.strftime(
 445                datetime.datetime.fromisoformat(v_setup["start_time"])
 446                - timedelta,
 447                "%Y-%m-%d %H:%M",
 448            )
 449            warming_model = self._model(setup=v_setup, mesh=self.mymesh.mesh)
 450            print(
 451                f"</> Getting parameters from the warming model model `warming_model` ..."
 452            )
 453            self._gathering_parameters(warming_model)
 454            warming_model.forward_run()
 455            self.validation_model.smash.rr_initial_states = (
 456                warming_model.rr_final_states.copy()
 457            )
 458
 459        self.validation_model.extra_smash_results = (
 460            self.validation_model.smash.forward_run(
 461                cost_options=cost_options,
 462                common_options=common_options,
 463                return_options=return_options,
 464            )
 465        )
 466
 467    def optimize(
 468        self,
 469        start_time: None | str = None,
 470        end_time: None | str = None,
 471        mapping: str = "uniform",
 472        optimizer: None | str = None,
 473        optimize_options: None | str = None,
 474        cost_options: None | str = None,
 475        common_options: None | str = None,
 476        return_options: None | str = None,
 477        callback=None,
 478    ):
 479        """
 480        Optimize the current model (with the current setup and mesh),
 481        store the model in the attribute optimize_model.smash and set the calibrated
 482        parameters to the model behind the attribute mysmashmodel.smash.
 483        Start_time and end_time can be specified here to change the period of the
 484        calibration compare to the current setup. All other arguments are
 485        equivalent to the smash.model.optimize function
 486        (see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize).
 487         One difference from Smash is that gauges with no data are
 488         automatically removed from the optimization without raising an error.
 489         This provide a convient way to calibrate quickly the parameters using the current
 490         mesh which may include gauges with data for calibration
 491         and location gauge for discharges computation.
 492
 493        :param start_time: The start time of the optimization, format "YYYY-mm-dd HH:MM", defaults to None
 494        :type start_time: None | str, optional
 495        :param end_time: The end time of the optimization, format "YYYY-mm-dd HH:MM", defaults to None
 496        :type end_time: None | str, optional
 497        :param mapping: Type of mapping, see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize, defaults to "uniform"
 498        :type mapping: str, optional
 499        :param optimizer: Name of optimizer, see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize, defaults to None
 500        :type optimizer: None | str, optional
 501        :param optimize_options: Dictionary containing optimization options for fine-tuning the optimization process, see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize, defaults to None
 502        :type optimize_options: None | str, optional
 503        :param cost_options: Dictionary containing computation cost options for simulated and observed responses. see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize, defaults to None
 504        :type cost_options: None | str, optional
 505        :param common_options: Dictionary containing common options with two elements: ncpu (int, default 1) and verbose (bool, default is False), defaults to None
 506        :type common_options: None | str, optional
 507        :param return_options: Dictionary containing return options to save additional simulation results, see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize, defaults to None
 508        :type return_options: None | str, optional
 509        :param callback: A callable called after each iteration with the signature callback(iopt: Optimize), see https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.optimize.html#smash.optimize, defaults to None
 510        :type callback: TYPE, optional
 511
 512        """
 513        o_setup = copy.deepcopy(self.mysetup.setup)
 514
 515        if start_time is not None and end_time is not None:
 516            o_setup["start_time"] = start_time
 517            o_setup["end_time"] = end_time
 518
 519        o_setup["read_qobs"] = True
 520        o_setup["read_prcp"] = True
 521        o_setup["read_pet"] = True
 522
 523        if self.optimize_model.smash is None:
 524            self.optimize_model.smash = self._model(
 525                setup=o_setup, mesh=self.mymesh.mesh
 526            )
 527            print(
 528                f"</> Getting parameters from the default model `mysmashmodel.smash` ..."
 529            )
 530
 531            for key in list(self.optimize_model.smash.rr_parameters.keys):
 532                pos = list(self.optimize_model.smash.rr_parameters.keys).index(
 533                    key
 534                )
 535                pos_w = list(self.mysmashmodel.smash.rr_parameters.keys).index(
 536                    key
 537                )
 538                self.optimize_model.smash.rr_parameters.values[:, :, pos] = (
 539                    self.mysmashmodel.smash.rr_parameters.values[:, :, pos_w]
 540                )
 541
 542        if self.warmup_model.smash is not None:
 543            for key in list(self.optimize_model.smash.rr_initial_states.keys):
 544                pos = list(
 545                    self.optimize_model.smash.rr_initial_states.keys
 546                ).index(key)
 547                pos_w = list(
 548                    self.warmup_model.smash.rr_initial_states.keys
 549                ).index(key)
 550                self.optimize_model.smash.rr_initial_states.values[
 551                    :, :, pos
 552                ] = self.warmup_model.smash.rr_final_states.values[:, :, pos_w]
 553
 554        default_cost_option = {
 555            "end_warmup": o_setup["start_time"],
 556            "gauge": "all",
 557        }
 558        if cost_options is not None:
 559            default_cost_option.update(cost_options)
 560
 561        # auto remove gauge if no observation
 562        if isinstance(default_cost_option["gauge"], str):
 563            if default_cost_option["gauge"] == "dws":
 564                gauge = np.empty(shape=0)
 565
 566                for i, pos in enumerate(
 567                    self.optimize_model.smash.mesh.gauge_pos
 568                ):
 569                    if (
 570                        self.optimize_model.smash.mesh.flwdst[tuple(pos)]
 571                        == 0.0
 572                    ):
 573                        gauge = np.append(
 574                            gauge, self.optimize_model.smash.mesh.code[i]
 575                        )
 576
 577            elif default_cost_option["gauge"] == "all":
 578                gauge = np.array(self.optimize_model.smash.mesh.code, ndmin=1)
 579            else:
 580                gauge = np.array(default_cost_option["gauge"], ndmin=1)
 581        elif isinstance(default_cost_option["gauge"], list):
 582            gauge = np.array(default_cost_option["gauge"], ndmin=1)
 583
 584        st = pd.Timestamp(self.optimize_model.smash.setup.start_time)
 585        et = pd.Timestamp(self.optimize_model.smash.setup.end_time)
 586        ew = pd.Timestamp(default_cost_option["end_warmup"])
 587        start_slice = int(
 588            (ew - st).total_seconds() / self.optimize_model.smash.setup.dt
 589        )
 590        # end_slice = start_slice+int((et - ew).total_seconds() / self.optimize_model.smash.setup.dt)
 591        end_slice = -1
 592        time_slice = slice(start_slice, end_slice)
 593
 594        del_gauge = []
 595        for i in range(len(gauge)):
 596            pos = np.where(self.optimize_model.smash.mesh.code == gauge[i])[0][
 597                0
 598            ]
 599            # print(i, pos)
 600            if np.all(
 601                self.optimize_model.smash.response_data.q[pos, time_slice] < 0
 602            ):
 603                del_gauge.append(i)
 604
 605                print(
 606                    f"No observed discharge available at gauge '{gauge[i]}' for the selected "
 607                    f"optimization period ['{ew}', '{et}']. This gauge is removed "
 608                    f"from the optimization."
 609                )
 610
 611        gauge = np.delete(gauge, del_gauge)
 612        default_cost_option.update({"gauge": gauge})
 613
 614        self.optimize_model.smash.optimize(
 615            mapping,
 616            optimizer,
 617            optimize_options,
 618            default_cost_option,
 619            common_options,
 620            return_options,
 621            callback,
 622        )
 623
 624        if (
 625            hasattr(self, "mysmashmodel")
 626            and self.mysmashmodel.smash is not None
 627        ):
 628            self.mysmashmodel.smash.rr_parameters = (
 629                self.optimize_model.smash.rr_parameters.copy()
 630            )
 631
 632    @tools.autocast_args
 633    def forward_run(
 634        self,
 635        warmup: int | None = None,
 636        invert_states: bool | None = False,
 637        cost_options: dict | None = None,
 638        common_options: dict | None = None,
 639        return_options: dict | None = None,
 640    ):
 641        """
 642        Smash model forward run.This function wrap smash.model.forward_run().
 643
 644        Parameters
 645        ----------
 646
 647        warmup: None | int
 648            a integer of the number of days used for warming the model.
 649        invert_states : bool = False
 650            invert states of the model, so that the final states are used
 651            for the initial states.
 652        cost_options : dict | None,
 653            Dictionary containing computation cost options for simulated and observed responses (https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.Model.forward_run.html).
 654        common_options: dict| None,
 655            Dictionary containing common options with two elements, ncpu and verbose (https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.Model.forward_run.html)
 656        return_options : dict | None,
 657            Dictionary containing return options to save additional simulation results. (https://smash.recover.inrae.fr/api_reference/principal_methods/smash/smash.Model.forward_run.html)
 658
 659        Examples
 660        --------
 661
 662        >>> es=smashbox.SmashBox()
 663        >>> sb.newmodel("RealCollobrier")
 664        >>> sb.RealCollobrier.generate_mesh(run=False)
 665        >>> sb.RealCollobrier.model()
 666        >>> sb.RealCollobrier.forward_run()
 667        """
 668
 669        if self.mysmashmodel.smash is None:
 670            self.model()
 671
 672        if warmup is not None:
 673            self.model_warmup(warmup=warmup)
 674
 675        if self.warmup_model.smash is not None:
 676            for key in list(self.mysmashmodel.smash.rr_initial_states.keys):
 677                pos = list(
 678                    self.mysmashmodel.smash.rr_initial_states.keys
 679                ).index(key)
 680                pos_w = list(
 681                    self.warmup_model.smash.rr_initial_states.keys
 682                ).index(key)
 683                self.mysmashmodel.smash.rr_initial_states.values[:, :, pos] = (
 684                    self.warmup_model.smash.rr_final_states.values[:, :, pos_w]
 685                )
 686
 687        if self.optimize_model.smash is not None:
 688            for key in list(self.mysmashmodel.smash.rr_parameters.keys):
 689                pos = list(self.mysmashmodel.smash.rr_parameters.keys).index(
 690                    key
 691                )
 692                pos_w = list(
 693                    self.optimize_model.smash.rr_parameters.keys
 694                ).index(key)
 695                self.mysmashmodel.smash.rr_parameters.values[:, :, pos] = (
 696                    self.optimize_model.smash.rr_parameters.values[:, :, pos_w]
 697                )
 698
 699        # needed here in order to replace the rainfall without rebuild the model
 700        self._gathering_atmosdata()
 701
 702        if invert_states:
 703            if self._fstates is not None and not np.all(
 704                self._fstates == -99.0
 705            ):
 706
 707                # TODO FIX BUG IN SMASH when states > 1.0
 708                mask_inf = np.where(
 709                    self.mysmashmodel.smash.mesh.active_cell == 1
 710                )
 711                mask_sup = np.where(self._fstates > 1)
 712                self._fstates[mask_sup] = 0.9999
 713
 714                self.mysmashmodel.smash.rr_initial_states.values[mask_inf] = (
 715                    self._fstates[mask_inf].copy()
 716                )
 717
 718            else:
 719                print(
 720                    "</> model final states does not exist."
 721                    " Invert states is not possible."
 722                    " Smash default initial states is used."
 723                )
 724
 725        self.mysmashmodel.extra_smash_results = (
 726            self.mysmashmodel.smash.forward_run(
 727                cost_options=cost_options,
 728                common_options=common_options,
 729                return_options=return_options,
 730            )
 731        )
 732
 733        # states backup
 734        self._fstates = self.mysmashmodel.smash.rr_final_states.values.copy()
 735        self._istates = self.mysmashmodel.smash.rr_initial_states.values.copy()
 736
 737    @tools.autocast_args
 738    def atmos_data_connector(
 739        self,
 740        input_prcp: np.ndarray | None = None,
 741        input_pet: np.ndarray | None = None,
 742        input_dt: float | None = None,
 743        input_res: tuple | list = (1000.0, 1000.0),
 744        input_start_time: str = "2050-01-01 01:00",
 745        input_bbox: dict | None = None,
 746        input_epsg: int = 2154,
 747        resampling_method="home_made_with_scipy_zoom",
 748    ):
 749        """
 750        Generate a connector for the external atmos_data comming from
 751        other model such as Graffas (spatial rainfall generator).
 752
 753        Parameters
 754        ----------
 755
 756        input_prcp : np.ndarray | None = None
 757            An np.ndarray with shape (nrow, ncol, npdt) wich stores the spatial rainfall
 758            which will be used by Smash.Ideally, the extend of this array match
 759            exactlly the extend of the Smash domain.
 760        input_pet : np.ndarray | None = None
 761            An np.ndarray with shape (nrow, ncol, npdt) wich stores the spatial
 762            evapotranspiration which will be used by Smash.
 763        input_dt : float = 3600.
 764            Time step in seconds of the input precipitation
 765        input_res : tuple | list = (1000., 1000.)
 766            The resolution of the input precipitation
 767        input_start_time : str = "2050-01-01 01:00"
 768            The date of the start_time
 769        input_bbox : dict | None = None
 770            The extend of the domain using a bbox.
 771            The convention used here is a dictionary like
 772            bbox={"left":xmin, "top": ymax, "right": xmax, "bottom": ymin}.
 773            If not provided, the extend of the smash domain will be used starting
 774            from xmin and ymin.
 775        input_epsg : int = 2154
 776            The epsg code of the coordinate system. If not provided,
 777            the coordinate system used in Smash will be used.
 778        resampling_method: str
 779            The method to use to resample and crop the input array. Default is
 780            'home_made_with_scipy_zoom' Choice are:
 781            ['rasterio_1', 'rasterio_2', 'home_made_with_scipy_zoom'].
 782            'home_made_with_scipy_zoom' is the fastest method. 'rasterio_1'
 783            is the slowest method. However, 'rasterio' method use much tested and
 784            reliable method to resample and crop the array.
 785
 786        Examples
 787        --------
 788
 789        >>> prcp_array=np.zeros(200,200,20)+1.
 790        >>> es=smashbox.SmashBox()
 791        >>> sb.newmodel("RealCollobrier")
 792        >>> sb.RealCollobrier.generate_mesh()
 793        >>> sb.RealCollobrier.atmos_data_connector(input_prcp=prcp_array)
 794
 795        """
 796        if input_prcp is None and input_pet is None:
 797            raise ValueError(
 798                "</> input_prcp and input_pet are None. input_prcp or input_pet must be"
 799                " numpy.ndarray with shape (nrow, ncol, ntimesteps)."
 800            )
 801
 802        if input_dt is None:
 803            print(
 804                "</> We suppose the input atmos data time-step equal to the model time-step."
 805            )
 806            input_dt = self.mysetup.setup["dt"]
 807
 808        if input_dt != self.mysetup.setup["dt"]:
 809            print(
 810                "</> Model and input atmos data have different time-step. "
 811                "Input atmos data will be resampled."
 812            )
 813            if input_prcp is not None:
 814                input_prcp = tools.time_resample_prcp_array(
 815                    input_prcp,
 816                    input_dt,
 817                    self.mysetup.setup["dt"],
 818                    t_axis=2,
 819                )
 820            if input_pet is not None:
 821                input_pet = tools.time_resample_prcp_array(
 822                    input_pet,
 823                    input_dt,
 824                    self.mysetup.setup["dt"],
 825                    t_axis=2,
 826                )
 827            input_dt = self.mysetup.setup["dt"]
 828
 829        # use conversion factor defined in setup
 830        if (
 831            "prcp_conversion_factor" in self.mysetup.setup
 832            and input_prcp is not None
 833        ):
 834            input_prcp = (
 835                input_prcp * self.mysetup.setup["prcp_conversion_factor"]
 836            )
 837        if (
 838            "pet_conversion_factor" in self.mysetup.setup
 839            and input_pet is not None
 840        ):
 841            input_pet = input_pet * self.mysetup.setup["pet_conversion_factor"]
 842
 843        self._myatmos_data_connector = (
 844            atmos_data_connector.atmos_data_connector(
 845                input_prcp=input_prcp,
 846                input_pet=input_pet,
 847                input_dt=input_dt,
 848                input_res=input_res,
 849                input_start_time=input_start_time,
 850                input_bbox=input_bbox,
 851                input_epsg=input_epsg,
 852            )
 853        )
 854
 855        smash_bbox = geo_toolbox.get_bbox_from_smash_mesh(self.mymesh.mesh)
 856
 857        if self.mymesh.mesh["nrow"] > input_prcp.shape[0]:
 858            print(
 859                f"</> Warning !! Input rainfall domain is smaller than the smash domain  {input_prcp.shape[0]}<{self.mymesh.mesh['nrow']}"
 860            )
 861        if self.mymesh.mesh["ncol"] > input_prcp.shape[1]:
 862            print(
 863                f"</> Warning !! Input rainfall domain is smaller than the smash domain  {input_prcp.shape[1]}<{self.mymesh.mesh['ncol']}"
 864            )
 865
 866        self._myatmos_data_connector.read_input_atmos_data(
 867            output_bbox=smash_bbox,
 868            output_epsg=self._myparam.param.epsg,
 869            output_res=(self.mymesh.mesh["xres"], self.mymesh.mesh["yres"]),
 870            output_shape=(self.mymesh.mesh["nrow"], self.mymesh.mesh["ncol"]),
 871            resampling_method=resampling_method,
 872        )
 873
 874        self._myatmos_data_connector.change_setup(self.mysetup)
 875
 876    def import_parameters(self, model=None):
 877        """
 878        Import Geotiff parameter in Smash. This function wrap
 879        smash.io.read_grid_parameters(). Path to the parameters is defined in
 880        self._myparam.param.smash_parameters.
 881
 882        Parameter
 883        ---------
 884
 885        model: smash.Model object
 886            A Smash model object. If None, the Smash model of smashbox will be used
 887
 888        Examples
 889        --------
 890
 891        >>> es=smashbox.SmashBox()
 892        >>> sb.newmodel("RealCollobrier")
 893        >>> sb.RealCollobrier.generate_mesh(run=False)
 894        >>> sb.RealCollobrier.model()
 895        >>> sb.RealCollobrier.import_parameters()
 896        >>> sb.RealCollobrier.forward_run()
 897        """
 898        if self._myparam.param.smash_parameters is None:
 899            print(
 900                "</> Warning: no calibrated Smash parameters is used, leaving it to"
 901                "  default."
 902            )
 903            return
 904
 905        if model is None:
 906            model = self.mysmashmodel.smash
 907
 908        smash.io.read_grid_parameters(
 909            model=model,
 910            path=self._myparam.param.smash_parameters,
 911        )
 912
 913    def export_parameters(self, path: os.PathLike = "./output_smash_param"):
 914        """
 915        Export Geotiff Smash parameter as Geotiff. This function wrap
 916        smash.io.save_grid_parameters().
 917
 918        Parameters
 919        ----------
 920
 921        path : os.PathLike = "./output_smash_param"
 922            path to a directory where the parameter will be saved.
 923
 924        Examples
 925        --------
 926
 927        >>> es=smashbox.SmashBox()
 928        >>> sb.newmodel("RealCollobrier")
 929        >>> sb.RealCollobrier.generate_mesh(run=False)
 930        >>> sb.RealCollobrier.model()
 931        >>> sb.RealCollobrier.import_parameters()
 932        >>> sb.RealCollobrier.export_parameters()
 933        """
 934        smash.io.export_parameters(self.mysmashmodel.smash, path)
 935
 936    def save_model_container_hdf5(
 937        self, path_to_hdf5: str | None = None, save_smash_model: bool = True
 938    ):
 939        """
 940
 941        :param path_to_hdf5: Path to the hdf5 file, defaults to None. If None,
 942        the fucntion will return a dictionnary containing all data of the self object.
 943        :type path_to_hdf5: str | None, optional
 944        :param save_smash_model: Savec the smash models objects or no, defaults
 945        to True. If False the Smash models objects are not saved.
 946        :type save_smash_model: str, optional
 947        :return: if a path_to_hdf5 is None, a dictionary containing all
 948        attribute of the input object self is returned.
 949        :rtype: dict | None
 950
 951        """
 952
 953        if path_to_hdf5 is not None:
 954
 955            structure = (
 956                pyhdf5_handler.src.object_handler.generate_object_structure(
 957                    self, include_method=False
 958                )
 959            )
 960
 961            if not save_smash_model:
 962                del structure["warmup_model"]["smash"]
 963                del structure["mysmashmodel"]["smash"]
 964                del structure["optimize_model"]["smash"]
 965                del structure["validation_model"]["smash"]
 966
 967            pyhdf5_handler.save_object_to_hdf5file(
 968                path_to_hdf5=path_to_hdf5,
 969                instance=self,
 970                keys_data=structure,
 971                location=f"./{self._model_name}/",
 972            )
 973
 974            # Add missing data
 975            structure = (
 976                pyhdf5_handler.src.object_handler.generate_object_structure(
 977                    self._myparam, include_method=False
 978                )
 979            )
 980
 981            pyhdf5_handler.save_object_to_hdf5file(
 982                path_to_hdf5=path_to_hdf5,
 983                instance=self._myparam,
 984                keys_data=structure,
 985                location=f"./{self._model_name}/",
 986            )
 987            pyhdf5_handler.save_dict_to_hdf5file(
 988                path_to_hdf5=path_to_hdf5,
 989                dictionary={"istates": self._istates},
 990                location=f"./{self._model_name}/",
 991            )
 992            pyhdf5_handler.save_dict_to_hdf5file(
 993                path_to_hdf5=path_to_hdf5,
 994                dictionary={"fstates": self._fstates},
 995                location=f"./{self._model_name}/",
 996            )
 997
 998        else:
 999            dict_results = (
1000                pyhdf5_handler.src.object_handler.read_object_as_dict(self)
1001            )
1002
1003            if not save_smash_model:
1004                del dict_results["warmup_model"]["smash"]
1005                del dict_results["mysmashmodel"]["smash"]
1006                del structure["optimize_model"]["smash"]
1007                del structure["validation_model"]["smash"]
1008
1009            return dict_results
1010
1011    def save_model_container(
1012        self,
1013        path: str | None = None,
1014        save_full_model=False,
1015    ):
1016        """
1017
1018        :param path: Path to the hdf5 file, defaults to None. If None,
1019        the fucntion will return a dictionnary containing all data of the self object.
1020        :type path: str | None, optional
1021
1022        """
1023
1024        if path is None:
1025            raise ValueError("Argument `path` is None.")
1026            return
1027
1028        if save_full_model is True:
1029            save_func = getattr(smash.io, "save_model")
1030        else:
1031            save_func = getattr(smash.io, "save_model_ddt")
1032
1033        if not os.path.exists(os.path.join(path, self._model_name)):
1034            os.makedirs(os.path.join(path, self._model_name))
1035
1036        structure = (
1037            pyhdf5_handler.src.object_handler.generate_object_structure(
1038                self, include_method=False
1039            )
1040        )
1041
1042        if "warmup_model" in structure:
1043            del structure["warmup_model"]["smash"]
1044        if "mysmashmodel" in structure:
1045            del structure["mysmashmodel"]["smash"]
1046        if "optimize_model" in structure:
1047            del structure["optimize_model"]["smash"]
1048        if "validation_model" in structure:
1049            del structure["validation_model"]["smash"]
1050
1051        pyhdf5_handler.save_object_to_hdf5file(
1052            path_to_hdf5=os.path.join(path, self._model_name, "smashbox.hdf5"),
1053            instance=self,
1054            keys_data=structure,
1055            location=f"./",
1056        )
1057
1058        structure = (
1059            pyhdf5_handler.src.object_handler.generate_object_structure(
1060                self._myparam, include_method=False
1061            )
1062        )
1063        pyhdf5_handler.save_object_to_hdf5file(
1064            path_to_hdf5=os.path.join(path, self._model_name, "smashbox.hdf5"),
1065            instance=self._myparam,
1066            keys_data=structure,
1067            location=f"./",
1068        )
1069
1070        pyhdf5_handler.save_dict_to_hdf5file(
1071            path_to_hdf5=os.path.join(path, self._model_name, "smashbox.hdf5"),
1072            dictionary={"istates": self._istates},
1073            location=f"./",
1074        )
1075        pyhdf5_handler.save_dict_to_hdf5file(
1076            path_to_hdf5=os.path.join(path, self._model_name, "smashbox.hdf5"),
1077            dictionary={"fstates": self._fstates},
1078            location=f"./",
1079        )
1080
1081        if self.mysmashmodel.smash is not None:
1082            save_func(
1083                self.mysmashmodel.smash,
1084                os.path.join(path, self._model_name, "mysmashmodel.hdf5"),
1085            )
1086        if self.warmup_model.smash is not None:
1087            save_func(
1088                self.warmup_model.smash,
1089                os.path.join(path, self._model_name, "warmup_model.hdf5"),
1090            )
1091        if self.optimize_model.smash is not None:
1092            save_func(
1093                self.optimize_model.smash,
1094                os.path.join(path, self._model_name, "optimize_model.hdf5"),
1095            )
1096        if self.validation_model.smash is not None:
1097            save_func(
1098                self.validation_model.smash,
1099                os.path.join(path, self._model_name, "validation_model.hdf5"),
1100            )
1101
1102    # def save_model_container_pickle(
1103    #     self,
1104    #     path: str | None = None,
1105    # ):
1106    #     """
1107
1108    #     :param path: Path to the hdf5 file, defaults to None. If None,
1109    #     the fucntion will return a dictionnary containing all data of the self object.
1110    #     :type path: str | None, optional
1111
1112    #     """
1113
1114    #     if path is not None:
1115    #         if not os.path.exists(os.path.join(path, self._model_name)):
1116    #             os.makedirs(os.path.join(path, self._model_name))
1117
1118    #         if self.mysmashmodel.smash is not None:
1119    #             smash.io.save_model(
1120    #                 self.mysmashmodel.smash,
1121    #                 os.path.join(path, self._model_name, "mysmashmodel.smash.hdf5"),
1122    #             )
1123    #         if self.warmup_model.smash is not None:
1124    #             smash.io.save_model(
1125    #                 self.warmup_model.smash,
1126    #                 os.path.join(path, self._model_name, "warmup_model.hdf5"),
1127    #             )
1128    #         if self.optimize_model.smash is not None:
1129    #             smash.io.save_model(
1130    #                 self.optimize_model.smash,
1131    #                 os.path.join(path, self._model_name, "optimize_model.hdf5"),
1132    #             )
1133
1134    #         structure = pyhdf5_handler.src.object_handler.generate_object_structure(
1135    #             self, include_method=False
1136    #         )
1137
1138    #         structure = [
1139    #             "mysetup",
1140    #             "mymesh",
1141    #             "mysmashmodel",
1142    #             "warmup_model",
1143    #             "optimize_model",
1144    #             "extra_smash_results",
1145    #             # "mystats", #input arg is self, and it is a weak ref (cannot pickle it)
1146    #             # "myplot", #input arg is self, and it is a weak ref (cannot pickle it)
1147    #             "_model_name",
1148    #             "_myparam",
1149    #             "_istates",
1150    #             "_fstates",
1151    #             "_myatmos_data_connector",
1152    #         ]
1153
1154    #         if "warmup_model" in structure:
1155    #             structure.remove("warmup_model")
1156    #         if "mysmashmodel" in structure:
1157    #             structure.remove("mysmashmodel")
1158    #         if "optimize_model" in structure:
1159    #             structure.remove("optimize_model")
1160
1161    #         for key in structure:
1162    #             print(key)
1163    #             with open(os.path.join(path, self._model_name, f"{key}.pkl"), "wb") as f:
1164    #                 obj = getattr(self, key)
1165    #                 if obj is not None:
1166    #                     pickle.dump(obj, f)
1167
1168    # def import_parameters(self)
1169    # if os.path.exists(self._myparam.param.smash_parameters) and self._myparam.param.smash_parameters.endswith(".hdf5"):
1170
1171    #     hdf5= pyhdf5_handler.open_hdf5(self._myparam.param.smash_parameters)
1172    #     hdf5_key=list(hdf5.keys())
1173    #     hdf5.close()
1174
1175    #     if "model_ddt" in hdf5_key:
1176    #         with_param = pyhdf5_handler.read_hdf5file_as_dict(self._myparam.param.smash_parameters,location="model_ddt/rr_parameters")
1177    #         with_mesh = pyhdf5_handler.get_hdf5file_item(self._myparam.param.smash_parameters,location="model_ddt",item="mesh")
1178
1179    #         smash_parameters.transfert_params_to_model(
1180    #             from_mesh=
1181    #             with_mesh,
1182    #             with_param=with_param,
1183    #             to_model=self.mysmashmodel.smash)
1184
1185    #     elif "model" in hdf5_key:
1186    #         with_param = pyhdf5_handler.get_hdf5file_item(self._myparam.param.smash_parameters,location="model",item="rr_parameters")
1187    #         with_mesh = pyhdf5_handler.get_hdf5file_item(self._myparam.param.smash_parameters,location="model",item="mesh")
1188
1189    #         smash_parameters.transfert_params_to_model(
1190    #             from_mesh=
1191    #             with_mesh,
1192    #             with_param=with_param,
1193    #             to_model=self.mysmashmodel.smash)
1194
1195    #     elif "rr_parameters" in hdf5.keys() and "mesh" in hdf5.keys():
1196    #         with_param = pyhdf5_handler.get_hdf5file_item(self._myparam.param.smash_parameters,location="./",item="rr_parameters")
1197    #         with_mesh = pyhdf5_handler.get_hdf5file_item(self._myparam.param.smash_parameters,location="./",item="mesh")
1198
1199    #         smash_parameters.transfert_params_to_model(
1200    #             from_mesh=with_mesh,
1201    #             with_param=with_param,
1202    #             to_model=self.mysmashmodel.smash)
1203
1204    #     else:
1205    #         raise ValueError("</> Error: cannot load the parameters ... unknown format...")
1206
1207    # elif os.path.isdir(self._myparam.param.smash_parameters):
1208
1209    #     smash_parameters.load_param_from_tiffformat(model=self.mysmashmodel.smash, path_to_parameters=self._myparam.param.smash_parameters)
1210
1211    # else:
1212    #     raise ValueError(f"</> Error: cannot load the parameters. Path or file '{self._myparam.param.smash_parameters}' does not exist ... ")
1213
1214    # def write_smashparam(self, hdf5file : str | os.PathLike = "output_smash_param.hdf5"):
1215
1216    #     rr_parameters=pyhdf5_handler.read_object_as_dict(self.mysmashmodel.smash.rr_parameters)
1217    #     pyhdf5_handler.save_dict_to_hdf5(hdf5file, rr_parameters)
1218
1219    # def export_parameters(self, path : os.path = "./output_smash_param")
1220
1221    # list_param=list(self.mysmashmodel.smash.rr_parameters.keys)
1222
1223    # for smparam in list_param:
1224
1225    #     array=self.mysmashmodel.smash.rr_parameters.values[:,:,list_param.index(smparam)]
1226
1227    #     smash_parameters.write_array_to_geotiff(os.path.join(path, smparam+".tif"),
1228    #                            array,
1229    #                            self.mysmashmodel.smash.mesh.xmin,
1230    #                            self.mysmashmodel.smash.mesh.ymax,
1231    #                            output_res = (self.mymesh.mesh["xres"],
1232    #                                          self.mymesh.mesh["yres"])
1233    #                            )