smashbox.test_average_stats

  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4Created on Tue Nov  4 14:11:45 2025
  5
  6@author: maxime
  7"""
  8
  9if __name__ == "__main__":
 10    ################ Test rapoort ################
 11    import smashbox
 12    import numpy as np
 13
 14    bbox = {
 15        "left": 875000.0,
 16        "bottom": 6228000.0,
 17        "right": 1001000.0,
 18        "top": 6320000.0,
 19    }
 20
 21    sb = smashbox.SmashBox()
 22
 23    sb.myparam.set_param("bbox", bbox)
 24
 25    sb.newmodel("rex")
 26
 27    sb.rex.mysetup.update_setup({"read_pet": False})
 28
 29    sb.rex.generate_mesh(
 30        query="(SURF>200) & (INFLUENCE=='Influence nulle ou faible')",
 31        area_error_th=0.2,
 32    )
 33
 34    ##########################################################
 35    # Artificial rainfall
 36    nrow = sb.rex.mymesh.mesh["nrow"]
 37    ncol = sb.rex.mymesh.mesh["ncol"]
 38    prcp = np.random.randint(-100, 100, size=(nrow, ncol, 24 * 30))
 39    prcp = np.where(prcp < 0, 0.0, prcp)
 40    ###########################################################
 41
 42    def model_creation(obj, prcp, qnoise):
 43
 44        obj.atmos_data_connector(input_prcp=prcp, input_dt=3600.0)
 45        obj.model()
 46
 47        obj.mysmashmodel.smash.atmos_data.pet = 0.0
 48        obj.forward_run(return_options={"q_domain": True}, invert_states=True)
 49
 50        ##########################################################
 51        # Artificial observed discharge
 52        rand = (
 53            np.random.randint(
 54                -qnoise, qnoise, obj.mysmashmodel.smash.response_data.q.shape
 55            )
 56            / 100.0
 57        )
 58        qsim = sb.rex.mysmashmodel.smash.response.q[:]
 59        obj.mysmashmodel.smash.response_data.q = qsim + qsim * rand
 60        ##########################################################
 61
 62    model_creation(sb.rex, prcp, 20.0)
 63
 64    sb.rex.mysmashmodel.mystats.fmisfit_stats()
 65    sb.rex.mysmashmodel.mystats.foutlets_stats()
 66    sb.rex.mysmashmodel.mystats.fspatial_stats()
 67    sb.rex.mysmashmodel.mystats.fquantile_stats(
 68        chunk_size=5,
 69        estimate_method="MLE",
 70        ncpu=6,
 71        fit="gumbel",
 72        compute_uncertainties=False,
 73    )
 74
 75    sb.copymodel("rex", "rex2")
 76    model_creation(sb.rex2, prcp * 1.5, 30.0)
 77
 78    sb.rex2.mysmashmodel.mystats.fmisfit_stats()
 79    sb.rex2.mysmashmodel.mystats.foutlets_stats()
 80    sb.rex2.mysmashmodel.mystats.fspatial_stats()
 81    # sb.rex2.mysmashmodel.mystats.quantile_stats = sb.rex.mystats.quantile_stats
 82    sb.rex2.mysmashmodel.mystats.fquantile_stats(
 83        chunk_size=5,
 84        estimate_method="MLE",
 85        ncpu=6,
 86        fit="gumbel",
 87        compute_uncertainties=False,
 88    )
 89
 90    sb.multimodel_statistics._get_model_list()
 91
 92    sb.multimodel_statistics.compute_multimodel_statistics()
 93
 94    sb.multimodel_statistics.compute_multimodel_statistics_misfit()
 95    sb.multimodel_statistics.compute_multimodel_statistics_outlets()
 96    sb.multimodel_statistics.compute_multimodel_statistics_outlets(obs=True)
 97    sb.multimodel_statistics.compute_multimodel_statistics_spatial()
 98    sb.multimodel_statistics.compute_multimodel_statistics_quantile()
 99
100    sb.multimodel_statistics.multimodel_quantile_stats.Quantile_1h.Q_th.data.shape
101    sb.multimodel_statistics.multimodel_quantile_stats.Quantile_1h.Q_th.data[
102        0, 60, 50, 0
103    ]
104    sb.multimodel_statistics.multimodel_quantile_stats.Quantile_1h.Q_th.data[
105        1, 60, 50, 0
106    ]
107    sb.rex.mysmashmodel.mystats.quantile_stats.Quantile_1h.Q_th[60, 50, 0]
108    sb.rex2.mysmashmodel.mystats.quantile_stats.Quantile_1h.Q_th[60, 50, 0]
109
110    from smashbox.plot import plot
111
112    fig, ax = plot.plot_image(
113        sb.multimodel_statistics.multimodel_quantile_stats.Quantile_1h.Q_th.median[
114            :, :, 0
115        ]
116    )
117    fig.show()
118    fig, ax = plot.plot_image(
119        sb.rex.mysmashmodel.mystats.quantile_stats.Quantile_1h.Q_th[:, :, 0]
120    )
121    fig.show()
122    fig, ax = plot.plot_image(
123        sb.rex2.mysmashmodel.mystats.quantile_stats.Quantile_1h.Q_th[:, :, 0]
124    )
125    fig.show()
126
127    fig, ax = plot.plot_misfit(
128        sb.multimodel_statistics.multimodel_misfit_stats.nse.mean,
129    )
130    fig.show()
131    fig, ax = plot.plot_misfit(
132        sb.rex.mysmashmodel.mystats.misfit_stats.results.nse,
133    )
134    fig.show()
135    fig, ax = plot.plot_misfit(
136        sb.rex2.mysmashmodel.mystats.misfit_stats.results.nse,
137    )
138    fig.show()