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