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