smashbox.tools.geo_toolbox

   1import numpy as np
   2from scipy.ndimage import zoom
   3import rasterio
   4
   5
   6# Convert coordinates to mesh matrix row-col (from smash)
   7def xy_to_rowcol(x, y, xmin, ymax, xres, yres):
   8    row = int((ymax - y) / yres)
   9    col = int((x - xmin) / xres)
  10    return row, col
  11
  12
  13# Convert mesh matrix row-col to geographic coordinates (from smash)
  14def rowcol_to_xy(row, col, xmin, ymax, xres, yres):
  15    x = int(col * xres + xmin)
  16    y = int(ymax - row * yres)
  17    return x, y
  18
  19
  20def get_bbox_from_smash_mesh(mesh):
  21    """
  22    Description
  23    -----------
  24    Compute the bbox from a Smash mesh dictionary
  25    Parameters
  26    ----------
  27    mesh: dict
  28        dict of smash mesh
  29    return
  30    ------
  31    dict()
  32        the bounding box of the smash mesh
  33    """
  34
  35    if "xres" in mesh and "yres" in mesh:
  36        dx = mesh["xres"]
  37        dy = mesh["yres"]
  38    else:
  39        dx = np.mean(mesh["dx"])
  40        dy = np.mean(mesh["dy"])
  41
  42    if "ncol" in mesh:
  43        ncol = mesh["ncol"]
  44        nrow = mesh["nrow"]
  45    else:
  46        nrow = mesh["active_cell"].shape[0]
  47        ncol = mesh["active_cell"].shape[1]
  48
  49    left = mesh["xmin"]
  50    right = mesh["xmin"] + ncol * dx
  51    bottom = mesh["ymax"] - nrow * dy
  52    top = mesh["ymax"]
  53    bbox = {"left": left, "bottom": bottom, "right": right, "top": top}
  54
  55    return bbox
  56
  57
  58def get_bbox_from_ascii_data(ascii_data):
  59    """
  60
  61    Description
  62    -----------
  63
  64    Compute the boundingbox from ascii-data stored in a dictionnary (must contain the header)
  65
  66    Parameters
  67    ----------
  68
  69    ascii_data: dict()
  70        Un dictionnaire contenant un format ascii: {ascii_data,xllcorner,yllcorner,ncol,nrows,cellsize}, la clé ascii_data contient un sous dictionnaire contenant des clés associés a des sources de données {clé:numpy_array}
  71
  72    return
  73    ------
  74
  75    dict()
  76        Un dictionnaire contenant la boundingbox (coordonnées): {left, bottom, right, top}
  77
  78    """
  79    ascii_data_l = {}
  80    for key, value in ascii_data.items():
  81        ascii_data_l.update({key.lower(): value})
  82
  83    bbox = {
  84        "left": ascii_data_l["xllcorner"],
  85        "bottom": ascii_data_l["yllcorner"],
  86        "right": ascii_data_l["xllcorner"]
  87        + ascii_data_l["ncols"] * ascii_data_l["cellsize"],
  88        "top": ascii_data_l["yllcorner"]
  89        + ascii_data_l["nrows"] * ascii_data_l["cellsize"],
  90    }
  91
  92    return bbox
  93
  94
  95def check_bbox_consistency(bbox_model_active_cell, bbox_param):
  96
  97    if bbox_model_active_cell["left"] < bbox_param["left"]:
  98        print(
  99            f"Warning: Model domain is larger than the domain of the parameter. {bbox_model_active_cell['left']}<{bbox_param['left']} (bbox model left < bbox param left). Expect lacuna (-99.) in model parameters"
 100        )
 101
 102    if bbox_model_active_cell["right"] > bbox_param["right"]:
 103        print(
 104            f"Warning: Model domain is larger than the domain of the parameter. {bbox_model_active_cell['right']}>{bbox_param['right']} (bbox model left < bbox param left). Expect lacuna (-99.) in model parameters"
 105        )
 106
 107    if bbox_model_active_cell["bottom"] < bbox_param["bottom"]:
 108        print(
 109            f"Warning: Model domain is larger than the domain of the parameter. {bbox_model_active_cell['bottom']}<{bbox_param['bottom']} (bbox model left < bbox param left). Expect lacuna (-99.) in model parameters"
 110        )
 111
 112    if bbox_model_active_cell["top"] < bbox_param["top"]:
 113        print(
 114            f"Warning: Model domain is larger than the domain of the parameter. {bbox_model_active_cell['top']}>{bbox_param['top']} (bbox model left < bbox param left). Expect lacuna (-99.) in model parameters"
 115        )
 116
 117
 118def intersection_bbox(bbox1, bbox2):
 119    """
 120
 121    Description
 122    -----------
 123
 124    Function which compute the bounding boxes intersection of 2 input bbox. It return the working bbox
 125
 126    Parameters
 127    ----------
 128
 129    bbox1: dict()
 130        containing the first bbox informations
 131    bbox2 : dict()
 132        containing the second bbox informations
 133
 134    returns
 135    ----------
 136
 137    dict()
 138        containing the bbox union
 139
 140    Examples
 141    ----------
 142
 143    dataset=gdal_raster_open(filename)
 144    possible_bbox=intersection_bbox(bbox,bbox_dataset)
 145
 146    """
 147    left = max(bbox1["left"], bbox2["left"])
 148    bottom = max(bbox1["bottom"], bbox2["bottom"])
 149    right = min(bbox1["right"], bbox2["right"])
 150    top = min(bbox1["top"], bbox2["top"])
 151    if (left < right) and (bottom < top):
 152        bbox_intersection = {
 153            "left": left,
 154            "bottom": bottom,
 155            "right": right,
 156            "top": top,
 157        }
 158        return bbox_intersection
 159    else:
 160        print("Impossible bounding boxes intersection")
 161        return {"left": 0, "bottom": 0, "right": 0, "top": 0}
 162
 163
 164def get_window_from_bbox(mesh, bbox):
 165    """
 166    Function to get the mesh window from a defined bbox
 167
 168    Parameters
 169    ----------
 170    dataset: gdal object
 171    bbox : dict containing the bbox
 172    ----------
 173    returns
 174    dic containing the computed windows
 175
 176    Examples
 177    ----------
 178    dataset=gdal_raster_open(filename)
 179    bbox_dataset=get_bbox(dataset)
 180    window=get_window_from_bbox(dataset,bbox_dataset)
 181
 182    """
 183
 184    if "xres" in mesh and "yres" in mesh:
 185        dx = mesh["xres"]
 186        dy = mesh["yres"]
 187    else:
 188        dx = np.mean(mesh["dx"])
 189        dy = np.mean(mesh["dy"])
 190
 191    col_off = (bbox["left"] - mesh["xmin"]) / dx
 192    row_off = (mesh["ymax"] - bbox["top"]) / dy
 193    ncols = (bbox["right"] - bbox["left"]) / dx
 194    nrows = (bbox["top"] - bbox["bottom"]) / dy
 195
 196    if (col_off < 0) or (row_off < 0):
 197        raise Exception(
 198            "The requested bounding box exceeds the limits of the raster domain."
 199        )
 200
 201    window = {
 202        "row_off": int(row_off),
 203        "col_off": int(col_off),
 204        "nrows": int(nrows),
 205        "ncols": int(ncols),
 206    }
 207
 208    return window
 209
 210
 211def get_cropped_window_from_bbox_intersection(
 212    bbox_intersection, bbox_origin, dx, dy
 213):
 214    """
 215
 216    Description
 217    -----------
 218
 219    Function to compute the domain to crop between a bbox intersection (included into bbox_origin) and the origin bbox. This function return a window such as the domain with bbox_intersection can be cropped using the retruned window according the bbox_origin
 220
 221    Parameters
 222    ----------
 223
 224    bbox_intersection: dict
 225        A bbox that intersect bbox_origin
 226
 227    bbox_origin: dict
 228        a bbox from which we want to extract data
 229
 230    dx: float
 231        size of the grid in the x direction
 232
 233    dy: float
 234        size of the grid in the y direction
 235
 236    Return
 237    ------
 238
 239    dict()
 240        a window dictionnary containing information to crop a matrix: {row_off, col_off, nrows, ncols}
 241
 242
 243    """
 244    if (
 245        (bbox_intersection["left"] < bbox_origin["left"])
 246        or (bbox_intersection["bottom"] < bbox_origin["bottom"])
 247        or (bbox_intersection["right"] > bbox_origin["right"])
 248        or (bbox_intersection["top"] > bbox_origin["top"])
 249    ):
 250        print(
 251            "The domain of bbox_intersection is not included in the domain of bbox_out"
 252        )
 253        window = {"row_off": 0, "col_off": 0, "nrows": 0, "ncols": 0}
 254        return window
 255
 256    col_off = (bbox_intersection["left"] - bbox_origin["left"]) / dx
 257    row_off = (bbox_origin["top"] - bbox_intersection["top"]) / dy
 258
 259    ncols = (bbox_intersection["right"] - bbox_intersection["left"]) / dx
 260    nrows = (bbox_intersection["top"] - bbox_intersection["bottom"]) / dy
 261
 262    window = {
 263        "row_off": int(row_off),
 264        "col_off": int(col_off),
 265        "nrows": int(nrows),
 266        "ncols": int(ncols),
 267    }
 268
 269    return window
 270
 271
 272def write_array_to_geotiff(filename, array, xmin, ymax, xres, yres):
 273
 274    metadata = {
 275        "driver": "GTiff",
 276        "dtype": "float64",
 277        "nodata": None,
 278        "width": array.shape[1],
 279        "height": array.shape[0],
 280        "count": 1,
 281        "crs": None,
 282        "transform": rasterio.Affine(xres, 0.0, xmin, 0.0, -yres, ymax),
 283    }
 284
 285    rasterio_write_tiff(filename=filename, matrix=array, metadata=metadata)
 286
 287
 288def rasterio_write_tiff(filename="test.tif", matrix=np.zeros(10), metadata={}):
 289
 290    with rasterio.Env():
 291        with rasterio.open(filename, "w", compress="lzw", **metadata) as dst:
 292            dst.write(matrix, 1)
 293
 294
 295def crop_array(
 296    array,
 297    bbox_in,
 298    res_in,
 299    bbox_out,
 300    res_out,
 301    order=0,
 302    cval=-99.0,
 303    grid_mode=True,
 304):
 305    """
 306
 307    Description
 308    --------------
 309
 310    Crop a part of a numpy array with an input bbox and resolution to a new bbox with a new resolution
 311
 312    Parameters
 313    ----------
 314
 315    array : numpy.array()
 316        Input gridded numpy array shape=(n,m)
 317
 318    bbox_in : dict()
 319        bounding box of the input array. Dictionnary {"left":,"top":,"right":,"bottom":}
 320
 321    res_in : dict()
 322        resolution of the input array in x and y direction. Disctionnary {"dx":, "dy":}
 323
 324    bbox_out : dict()
 325        bounding box of the output array. Dictionnary {"left":,"top":,"right":,"bottom":}
 326
 327    res_out : dict()
 328        resolution of the output array in x and y direction. Disctionnary {"dx":, "dy":}
 329
 330    order : int()
 331        order of the resampling cubic interpolation
 332
 333    cval : float() | np.nan
 334        fill value for the extended boundaries
 335
 336    grid_mode : bool()
 337        True | False. if True coordinate start from the edge of the cell. If False coordinate starts from the center of the cell.
 338
 339    Return
 340    ------
 341
 342    numpy.array()
 343        Cropped and resampled array according bbox_out and res_out
 344
 345    """
 346
 347    # intersection bbox
 348    bbox_intersection = intersection_bbox(bbox_in, bbox_out)
 349
 350    # ---------------------- fait une découpe grossière du domaine -------
 351
 352    # 1- crop array on bbox_intersection+-res_in at res in: shrink the domain, speed up futur resampling on large domain ?
 353    if res_in["dx"] >= res_out["dx"] and res_in["dy"] >= res_out["dy"]:
 354        res_shrinked = {"dx": res_in["dx"], "dy": res_in["dy"]}
 355    elif res_in["dx"] < res_out["dx"] and res_in["dy"] < res_out["dy"]:
 356        res_shrinked = {"dx": res_out["dy"], "dy": res_out["dy"]}
 357    elif res_in["dx"] < res_out["dx"] and res_in["dy"] >= res_out["dy"]:
 358        res_shrinked = {"dx": res_out["dx"], "dy": res_in["dy"]}
 359    elif res_in["dx"] >= res_out["dx"] and res_in["dy"] < res_out["dy"]:
 360        res_shrinked = {"dx": res_in["dx"], "dy": res_out["dy"]}
 361
 362    bbox_in_shrinked = {
 363        "left": bbox_in["left"]
 364        + int(
 365            max(bbox_intersection["left"] - bbox_in["left"], 0)
 366            / res_shrinked["dx"]
 367        )
 368        * res_shrinked["dx"],
 369        "right": bbox_in["right"]
 370        - int(
 371            max(bbox_in["right"] - bbox_intersection["right"], 0)
 372            / res_shrinked["dx"]
 373        )
 374        * res_shrinked["dx"],
 375        "bottom": bbox_in["bottom"]
 376        + int(
 377            max(bbox_intersection["bottom"] - bbox_in["bottom"], 0)
 378            / res_shrinked["dy"]
 379        )
 380        * res_shrinked["dy"],
 381        "top": bbox_in["top"]
 382        - int(
 383            max(bbox_in["top"] - bbox_intersection["top"], 0)
 384            / res_shrinked["dy"]
 385        )
 386        * res_shrinked["dy"],
 387    }
 388    bbox_intersection_shrinked = intersection_bbox(bbox_in, bbox_in_shrinked)
 389
 390    windows_wrap = get_window_from_bbox(
 391        mesh={
 392            "xmin": bbox_in["left"],
 393            "ymax": bbox_in["top"],
 394            "xres": res_shrinked["dx"],
 395            "yres": res_shrinked["dy"],
 396        },
 397        bbox=bbox_intersection_shrinked,
 398    )
 399
 400    # Erase input array and bbox_in
 401    array = array[
 402        windows_wrap["row_off"] : windows_wrap["row_off"]
 403        + windows_wrap["nrows"],
 404        windows_wrap["col_off"] : windows_wrap["col_off"]
 405        + windows_wrap["ncols"],
 406    ]
 407
 408    bbox_in = bbox_intersection_shrinked
 409
 410    # ---------------------------------------------------------------------------------------
 411
 412    # 2- resample the array to res_out
 413    resampled_array = resample_array(
 414        array,
 415        res_in=res_in,
 416        res_out=res_out,
 417        order=order,
 418        cval=cval,
 419        grid_mode=grid_mode,
 420    )
 421
 422    # 3- crop the array on the intersection of bbox_in and bbox_out
 423    # window of bbox_intersection in the domain of bbox_in (bbox_prcp, matrix to read)
 424    window_intersection = get_cropped_window_from_bbox_intersection(
 425        bbox_intersection, bbox_in, res_out["dx"], res_out["dy"]
 426    )
 427
 428    # reading the part of the matrix (array_in)
 429    cropped_array = resampled_array[
 430        window_intersection["row_off"] : window_intersection["row_off"]
 431        + window_intersection["nrows"],
 432        window_intersection["col_off"] : window_intersection["col_off"]
 433        + window_intersection["ncols"],
 434    ]
 435
 436    # allocate out array: shape of bbox_out
 437    array_out = (
 438        np.zeros(
 439            shape=(
 440                int((bbox_out["top"] - bbox_out["bottom"]) / res_out["dx"]),
 441                int((bbox_out["right"] - bbox_out["left"]) / res_out["dy"]),
 442            )
 443        )
 444        - 99.0
 445    )
 446
 447    # window of bbox_intersection in the domain of bbox_smash
 448    window_intersection_out = get_cropped_window_from_bbox_intersection(
 449        bbox_intersection, bbox_out, res_out["dx"], res_out["dy"]
 450    )
 451
 452    # prcp domain smaller than smash domain:
 453    if window_intersection_out["ncols"] > cropped_array.shape[1]:
 454        window_intersection_out["ncols"] = cropped_array.shape[1]
 455    if window_intersection_out["nrows"] > cropped_array.shape[0]:
 456        window_intersection_out["nrows"] = cropped_array.shape[0]
 457
 458    # copy crop_array (input matrix cropped) in array_out
 459    array_out[
 460        window_intersection_out["row_off"] : window_intersection_out["row_off"]
 461        + window_intersection_out["nrows"],
 462        window_intersection_out["col_off"] : window_intersection_out["col_off"]
 463        + window_intersection_out["ncols"],
 464    ] = cropped_array
 465
 466    return array_out
 467
 468
 469def resample_array(
 470    array,
 471    res_in={"dx": 1, "dy": 1},
 472    res_out={"dx": 1, "dy": 1},
 473    order=0,
 474    cval=-99.0,
 475    grid_mode=True,
 476):
 477    """
 478
 479    Parameters
 480    ----------
 481
 482    array: numpy.array()
 483        Input gridded numpy array shape=(n,m)
 484
 485    res_in: dict()
 486        resolution of the input array in x and y direction. Disctionnary {"dx":, "dy":}
 487
 488    res_out: dict()
 489        resolution of the output array in x and y direction. Disctionnary {"dx":, "dy":}
 490
 491    order: int()
 492        order of the resampling cubic interpolation
 493
 494    cval: float() | np.nan
 495        fill value for the extended boundaries
 496
 497    grid_mode: bool()
 498        True | False. if True coordinate start from the edge of the cell. If False coordinate starts from the center of the cell.
 499
 500    Return
 501    ------
 502
 503    numpy.array()
 504        Cropped and resampled array according bbox_out and res_out
 505
 506    """
 507
 508    ratio_x = res_in["dx"] / res_out["dx"]
 509    ratio_y = res_in["dy"] / res_out["dy"]
 510    resampled_array = zoom(
 511        array,
 512        (ratio_y, ratio_x),
 513        order=order,
 514        mode="grid-constant",
 515        cval=-99.0,
 516        grid_mode=True,
 517    )
 518
 519    return resampled_array
 520
 521
 522# def read_geotiff_to_ascii(path=""):
 523
 524#     if not os.path.exists(path):
 525#         raise ValueError(f"{path} does not exist.")
 526
 527#     with rasterio.open(path) as ds:
 528#         transform = ds.get_transform()
 529#         xmin = transform[0]
 530#         ymax = transform[3]
 531#         xres = transform[1]
 532#         yres = -transform[5]
 533
 534#         ascii_data = {
 535#             "ncols": ds.width,
 536#             "nrows": ds.height,
 537#             "cellsize": xres,
 538#             "xllcorner": xmin,
 539#             "yllcorner": ymax - ds.height * yres,
 540#             "nodata_value": ds.nodata,
 541#             "data": ds.read(indexes=1),
 542#         }
 543
 544#     return ascii_data
 545
 546
 547# def write_smashparam_to_geotiff(model, path):
 548
 549#     if not os.path.exists(path):
 550#         os.mkdir(path)
 551
 552#     list_param = list(model.rr_parameters.keys)
 553
 554#     for param in list_param:
 555
 556#         array = model.rr_parameters.values[:, :, list_param.index(param)]
 557
 558#         write_array_to_geotiff(
 559#             os.path.join(path, param + ".tif"),
 560#             array,
 561#             model.mesh.xmin,
 562#             model.mesh.ymax,
 563#             model.mesh.xres,
 564#             model.mesh.yres,
 565#         )
 566
 567
 568# # format standart smashop ?
 569# # faire une fonction qui extrait des params depuis smash vers un dict_ascii_format ?
 570# def load_param_from_tiffformat(
 571#     model=None, path_to_parameters="", bounding_condition=False
 572# ):
 573
 574#     # check metadata first
 575
 576#     list_param = model.rr_parameters.keys
 577
 578#     mesh_out = object_handler.read_object_as_dict(model.mesh)
 579
 580#     bbox_out = get_bbox_from_smash_mesh(mesh_out)
 581#     res_out = {"dx": np.mean(mesh_out["dx"]), "dy": np.mean(mesh_out["dy"])}
 582
 583#     for param in list_param:
 584
 585#         if os.path.exists(os.path.join(path_to_parameters, param + ".tif")):
 586
 587#             ascii_data = read_geotiff_to_ascii(
 588#                 os.path.join(path_to_parameters, param + ".tif")
 589#             )
 590
 591#             bbox_in = get_bbox_from_ascii_data(ascii_data)
 592
 593#             check_bbox_consistency(bbox_out, bbox_in)
 594
 595#             res_in = {"dx": ascii_data["cellsize"], "dy": ascii_data["cellsize"]}
 596
 597#             print(f"</> In read_param_from_asciiformat, reading param {param}")
 598
 599#             cropped_param = crop_array(
 600#                 ascii_data["data"],
 601#                 bbox_in,
 602#                 res_in,
 603#                 bbox_out,
 604#                 res_out,
 605#                 resampling_method="scipy-zoom",
 606#                 order=0,
 607#                 cval=-99.0,
 608#                 grid_mode=True,
 609#             )
 610
 611#             pos = np.where(list_param == param)[0][0]
 612#             model.rr_parameters.values[:, :, pos] = cropped_param
 613
 614#             print("</> Removing values lower than 0")
 615#             model.rr_parameters.values[:, :, pos] = np.where(
 616#                 model.rr_parameters.values[:, :, pos] < 0,
 617#                 0,
 618#                 model.rr_parameters.values[:, :, pos],
 619#             )
 620
 621#         else:
 622
 623#             print(
 624#                 f"</> Error: in load_param_from_tiffformat, missing parameter {param} in {path_to_parameters}"
 625#             )
 626
 627
 628# # format standart smashop ?
 629# # faire une fonction qui extrait des params depuis smash vers un dict_ascii_format ?
 630# def load_param_from_asciiformat(model, dict_ascii_parameters, bounding_condition=False):
 631
 632#     list_param = model.rr_parameters.keys
 633
 634#     mesh_out = object_handler.read_object_as_dict(model.mesh)
 635
 636#     bbox_in = get_bbox_from_ascii_data(dict_ascii_parameters)
 637#     res_in = {
 638#         "dx": dict_ascii_parameters["cellsize"],
 639#         "dy": dict_ascii_parameters["cellsize"],
 640#     }
 641
 642#     bbox_out = get_bbox_from_smash_mesh(mesh_out)
 643#     res_out = {"dx": np.mean(mesh_out["dx"]), "dy": np.mean(mesh_out["dy"])}
 644
 645#     check_bbox_consistency(bbox_out, bbox_in)
 646
 647#     for param in list_param:
 648
 649#         if param in dict_ascii_parameters["ascii_data"].keys():
 650
 651#             print(f"</> In read_param_from_asciiformat, reading param {param}")
 652
 653#             cropped_param = crop_array(
 654#                 param,
 655#                 bbox_in,
 656#                 res_in,
 657#                 bbox_out,
 658#                 res_out,
 659#                 resampling_method="scipy-zoom",
 660#                 order=0,
 661#                 cval=-99.0,
 662#                 grid_mode=True,
 663#             )
 664
 665#             pos = np.where(list_param == param)[0][0]
 666#             model.rr_parameters.values[:, :, pos] = cropped_param
 667#         else:
 668#             print(
 669#                 f"</> In read_param_from_asciiformat, skipping missing parameter {param}"
 670#             )
 671
 672
 673# def transfert_params_to_model(
 674#     from_mesh=dict(), with_param=(dict | object), to_model=None
 675# ):
 676
 677#     if isinstance(with_param, dict):
 678#         pass
 679#     else:
 680#         with_param = object_handler.read_object_as_dict(with_param)
 681
 682#     structure = to_model.setup.structure
 683#     structure_parameters = smash._constant.STRUCTURE_RR_PARAMETERS.copy()
 684#     # copy() obligatoire car la liste contenu dans la copie du dictionnaire n'est pas une copy mais un pointeur ...
 685#     parameters_list = structure_parameters[structure].copy()
 686
 687#     mesh_out = object_handler.read_object_as_dict(to_model.mesh)
 688
 689#     bbox_in = get_bbox_from_smash_mesh(from_mesh)
 690
 691#     if "xres" in from_mesh and "yres" in from_mesh:
 692#         dx = from_mesh["xres"]
 693#         dy = from_mesh["yres"]
 694#     else:
 695#         dx = np.mean(from_mesh["dx"])
 696#         dy = np.mean(from_mesh["dy"])
 697
 698#     res_in = {
 699#         "dx": dx,
 700#         "dy": dy,
 701#     }
 702
 703#     bbox_out = get_bbox_from_smash_mesh(mesh_out)
 704
 705#     if "xres" in mesh_out and "yres" in mesh_out:
 706#         dx = mesh_out["xres"]
 707#         dy = mesh_out["yres"]
 708#     else:
 709#         dx = np.mean(mesh_out["dx"])
 710#         dy = np.mean(mesh_out["dy"])
 711
 712#     res_out = {
 713#         "dx": dx,
 714#         "dy": dy,
 715#     }
 716
 717#     for i in range(len(parameters_list)):
 718
 719#         if "values" in with_param:
 720#             param = with_param["values"][:, :, i]
 721#             pos = i
 722#         else:
 723#             param = with_param[parameters_list[i]]
 724#             pos = np.where(to_model.rr_parameters.keys == parameters_list[i])[0][0]
 725
 726#         cropped_param = crop_array(
 727#             param,
 728#             bbox_in,
 729#             res_in,
 730#             bbox_out,
 731#             res_out,
 732#             resampling_method="scipy-zoom",
 733#             order=0,
 734#             cval=-99.0,
 735#             grid_mode=True,
 736#         )
 737
 738#         to_model.rr_parameters.values[:, :, pos] = cropped_param
 739
 740
 741# def resample_array(
 742#     array,
 743#     res_in={"dx": 1, "dy": 1},
 744#     res_out={"dx": 1, "dy": 1},
 745#     # resampling_method="scipy-zoom",
 746#     order=3,
 747#     cval=-99.0,
 748#     grid_mode=True,
 749# ):
 750# """
 751
 752# Description
 753# -----------
 754
 755# Resample an array with input resolution to the output resolution. It use zoom function form scipy package.
 756
 757# Parameters
 758# ----------
 759
 760# array: numpy.array()
 761#     Input gridded numpy array shape=(n,m)
 762
 763# res_in: dict()
 764#     resolution of the input array in x and y direction. Disctionnary {"dx":, "dy":}
 765
 766# res_out: dict()
 767#     resolution of the output array in x and y direction. Disctionnary {"dx":, "dy":}
 768
 769# resampling_method: str()
 770#     scipy-zoom | numpy ; scipy-zoom use the zoom function from scipy (see function resample_array_scipy_zoom for other input arg). numpy only resample using basic function np.repeat ans basic calculation (not suitable for every case).
 771
 772# Return
 773# ------
 774
 775# numpy.array()
 776#     Resampled array according res_out
 777
 778# """
 779# if resampling_method == "numpy":
 780#     resampled_array = resample_array_numpy(array, res_in=res_in, res_out=res_out)
 781
 782#     if resampled_array is None:
 783#         print("</>Fallback to scipy-zoom resampling (scipy.ndimage.zoom method)")
 784#         resampling_method = "scipy-zoom"
 785
 786# if resampling_method == "scipy-zoom":
 787#     resampled_array = resample_array_scipy_zoom(
 788#         array, res_in, res_out, order=0, cval=-99.0, grid_mode=True
 789#     )
 790
 791# return resampled_array
 792
 793
 794# def resample_array_numpy(array, res_in={"dx": 1, "dy": 1}, res_out={"dx": 1, "dy": 1}):
 795#     """
 796
 797#     Description
 798#     -----------
 799
 800#     resample an array with input resolution to the output resolution. It simpy use numpy functions.
 801
 802#     Parameters
 803#     ----------
 804
 805#     array: numpy.array()
 806#         Input gridded numpy array shape=(n,m)
 807
 808#     res_in: dict()
 809#         resolution of the input array in x and y direction. Disctionnary {"dx":, "dy":}
 810
 811#     res_out: dict()
 812#         resolution of the output array in x and y direction. Disctionnary {"dx":, "dy":}
 813
 814#     Return
 815#     ------
 816
 817#     numpy.array()
 818#         Resampled array according res_out
 819
 820#     """
 821#     # 1Resampling
 822#     if (res_in["dx"] >= res_out["dx"]) and (res_in["dy"] >= res_out["dy"]):
 823#         ratio_x = int(res_in["dx"] / res_out["dx"])
 824#         ratio_y = int(res_in["dy"] / res_out["dy"])
 825#         # resampling input matrix to dx,dy smash si résolution demandée est inférieur
 826#         resampled_array = np.repeat(np.repeat(array, ratio_x, axis=0), ratio_y, axis=1)
 827
 828#     elif (res_in["dx"] < res_out["dx"]) and (res_in["dy"] < res_out["dy"]):
 829
 830#         ratio_x = int(res_out["dx"] / res_in["dx"])
 831#         ratio_y = int(res_out["dy"] / res_in["dy"])
 832#         # resampling input matrix si résolution demandée est supérieur:
 833#         xsize = int(np.ceil(array.shape[0] / ratio_x))
 834#         ysize = int(np.ceil(array.shape[1] / ratio_y))
 835#         resampled_array = np.full(shape=(xsize, ysize), fill_value=0)
 836
 837#         xr = 0
 838#         for x in range(0, xsize):
 839#             yr = 0
 840#             for y in range(0, ysize):
 841#                 buffer = array[xr : xr + ratio_x, yr : yr + ratio_y]
 842#                 average = np.mean(buffer, axis=(0, 1))
 843#                 resampled_array[x, y, :] = average
 844#                 yr = yr + ratio_y
 845#             xr = xr + ratio_x
 846#     else:
 847#         print(
 848#             "</> average resampling method impossible, different resampling resolution in x and y axis."
 849#         )
 850#         return None
 851
 852#     return resampled_array
 853
 854
 855# def import_parameters(model: Model, path_to_parameters: FilePath):
 856#     """
 857#     Description
 858#     -----------
 859#     Read a geotif, resample if necessarry then clip it on the bouning box of the smash mesh
 860
 861#     Parameters
 862#     ----------
 863#     model: object
 864#         SMASH model object
 865#     path_to_parameters: str
 866#         Path to the directory which contain the geotiff files (parameters)
 867
 868#     return
 869#     ------
 870#     np.ndarray
 871#         The data clipped on the SMASH bounding box
 872#     """
 873#     list_param = model.rr_parameters.keys
 874
 875#     for param in list_param:
 876#         if os.path.exists(os.path.join(path_to_parameters, param + ".tif")):
 877#             cropped_param = _rasterio_read_param2(
 878#                 path=os.path.join(path_to_parameters, param + ".tif"), mesh=model.mesh
 879#             )
 880
 881#             pos = np.argwhere(list_param == param).item()
 882#             model.rr_parameters.values[:, :, pos] = cropped_param
 883
 884#         else:
 885#             raise ValueError(f"Missing parameter {param} in {path_to_parameters}")
 886
 887
 888# def _rasterio_read_param(path: FilePath, mesh: MeshDT):
 889#     """
 890#     Description
 891#     -----------
 892#     Read a geotif, resample if necessarry then clip it on the bouning box of the smash mesh
 893
 894#     Parameters
 895#     ----------
 896#     path: str
 897#         Path to a geotiff file.
 898#     mesh: object
 899#         object of the smash mesh
 900
 901#     return
 902#     ------
 903#     np.ndarray
 904#         The data clipped on the SMASH bounding box
 905#     """
 906#     bounds = get_bbox_from_smash_mesh(mesh)
 907#     xres = mesh.xres
 908#     yres = mesh.yres
 909
 910#     # Open the larger raster
 911#     with rasterio.open(path) as dataset:
 912#         x_scale_factor = dataset.res[0] / xres
 913#         y_scale_factor = dataset.res[1] / yres
 914
 915#         # resampling first to avoid spatial shifting of the parameters
 916#         data = dataset.read(
 917#             out_shape=(
 918#                 dataset.count,
 919#                 int(dataset.height * y_scale_factor),
 920#                 int(dataset.width * x_scale_factor),
 921#             ),
 922#             resampling=Resampling.nearest,
 923#         )
 924
 925#     bbox_dataset = {
 926#         "left": dataset.bounds.left,
 927#         "bottom": dataset.bounds.bottom,
 928#         "right": dataset.bounds.right,
 929#         "top": dataset.bounds.top,
 930#     }
 931
 932#     bbox_intersection = intersection_bbox(bbox_dataset, bounds)
 933
 934#     window_intersection = get_cropped_window_from_bbox_intersection(
 935#         bbox_intersection, bbox_dataset, xres, yres
 936#     )
 937
 938#     cropped_array = data[
 939#         0,
 940#         window_intersection["row_off"] : window_intersection["row_off"]
 941#         + window_intersection["nrows"],
 942#         window_intersection["col_off"] : window_intersection["col_off"]
 943#         + window_intersection["ncols"],
 944#     ]
 945
 946#     # allocate out array: shape of bbox_out
 947#     array_out = np.zeros(
 948#         shape=(
 949#             int((bounds["top"] - bounds["bottom"]) / xres),
 950#             int((bounds["right"] - bounds["left"]) / yres),
 951#         )
 952#     )
 953
 954#     # window of bbox_intersection in the domain of bbox_smash
 955#     window_intersection_out = get_cropped_window_from_bbox_intersection(
 956#         bbox_intersection, bounds, xres, yres
 957#     )
 958
 959#     # copy crop_array (input matrix cropped) in array_out
 960#     array_out[
 961#         window_intersection_out["row_off"] : window_intersection_out["row_off"]
 962#         + window_intersection_out["nrows"],
 963#         window_intersection_out["col_off"] : window_intersection_out["col_off"]
 964#         + window_intersection_out["ncols"],
 965#     ] = cropped_array
 966
 967#     return array_out
 968
 969
 970# def _rasterio_read_param2(path: FilePath, mesh: MeshDT):
 971#     """
 972#     Description
 973#     -----------
 974#     Read a geotif, resample if necessarry then clip it on the bouning box of the smash mesh
 975
 976#     Parameters
 977#     ----------
 978#     path: str
 979#         Path to a geotiff file.
 980#     mesh: object
 981#         object of the smash mesh
 982
 983#     return
 984#     ------
 985#     np.ndarray
 986#         The data clipped on the SMASH bounding box
 987#     """
 988#     input_bbox = get_bbox_from_smash_mesh(mesh)
 989
 990#     xres = mesh.xres
 991#     yres = mesh.yres
 992
 993#     output_crs = rasterio.CRS.from_epsg(mesh.epsg)
 994
 995#     # Open the larger raster
 996#     with rasterio.open(path) as dataset:
 997
 998#         x_scale_factor = dataset.res[0] / xres
 999#         y_scale_factor = dataset.res[1] / yres
1000
1001#         transform = dataset.transform
1002#         height = dataset.height
1003#         width = dataset.width
1004#         crs = dataset.crs
1005
1006#         # resampling first to avoid spatial shifting of the parameters
1007#         data = dataset.read(
1008#             out_shape=(
1009#                 dataset.count,
1010#                 int(dataset.height * y_scale_factor),
1011#                 int(dataset.width * x_scale_factor),
1012#             ),
1013#             resampling=Resampling.nearest,
1014#         )
1015
1016#     # Use a memory dataset
1017#     with rasterio.io.MemoryFile() as memfile:
1018
1019#         with memfile.open(
1020#             driver="GTiff",
1021#             height=height,
1022#             width=width,
1023#             count=1,
1024#             dtype=data.dtype,
1025#             transform=transform,
1026#             crs=crs,
1027#         ) as dataset:
1028#             dataset.write(data[0, :, :], 1)
1029
1030#             new_width = int((input_bbox["right"] - input_bbox["left"]) / yres)
1031#             new_height = int((input_bbox["top"] - input_bbox["bottom"]) / xres)
1032#             new_transform = rasterio.transform.from_bounds(
1033#                 west=input_bbox["left"],
1034#                 south=input_bbox["bottom"],
1035#                 east=input_bbox["right"],
1036#                 north=input_bbox["top"],
1037#                 width=new_width,
1038#                 height=new_height,
1039#             )
1040
1041#             # Target array
1042#             new_array = np.empty((new_height, new_width), dtype=np.float32)
1043
1044#             # reproject dataset
1045#             rasterio.warp.reproject(
1046#                 source=rasterio.band(dataset, 1),
1047#                 destination=new_array,
1048#                 src_transform=transform,
1049#                 src_crs=crs,
1050#                 dst_transform=new_transform,
1051#                 dst_crs=output_crs,
1052#                 resampling=Resampling.nearest,
1053#             )
1054
1055#     return new_array