def square_domain(R, metadata, method="pad", inverse=False):
    """
    Either pad or crop a field to obtain a square domain.

    Parameters
    ----------
    R: array-like
        Array of shape (m,n) or (t,m,n) containing the input fields.
    metadata: dict
        Metadata dictionary containing the x1, x2, y1, y2,
        xpixelsize, ypixelsize,
        attributes as described in the documentation of
        :py:mod:`pysteps.io.importers`.
    method: {'pad', 'crop'}, optional
        Either pad or crop.
        If pad, an equal number of zeros is added to both ends of its shortest
        side in order to produce a square domain.
        If crop, an equal number of pixels is removed
        to both ends of its longest side in order to produce a square domain.
        Note that the crop method involves an irreversible loss of data.
    inverse: bool, optional
        Perform the inverse method to recover the original domain shape.
        After a crop, the inverse is performed by padding the field with zeros.

    Returns
    -------
    R: array-like
        the reshape dataset
    metadata: dict
        the metadata with updated attributes.
    """

    R = R.copy()
    R_shape = np.array(R.shape)
    metadata = metadata.copy()

    if not inverse:
        if len(R.shape) < 2:
            raise ValueError("The number of dimension must be > 1")
        if len(R.shape) == 2:
            R = R[None, None, :]
        if len(R.shape) == 3:
            R = R[None, :]
        if len(R.shape) > 4:
            raise ValueError("The number of dimension must be <= 4")

        if R.shape[2] == R.shape[3]:
            return R.squeeze()

        orig_dim = R.shape
        orig_dim_n = orig_dim[0]
        orig_dim_t = orig_dim[1]
        orig_dim_y = orig_dim[2]
        orig_dim_x = orig_dim[3]

        if method == "pad":
            new_dim = np.max(orig_dim[2:])
            R_ = np.ones((orig_dim_n, orig_dim_t, new_dim, new_dim)) * R.min()

            if orig_dim_x < new_dim:
                idx_buffer = int((new_dim - orig_dim_x) / 2.0)
                R_[:, :, :, idx_buffer : (idx_buffer + orig_dim_x)] = R
                metadata["x1"] -= idx_buffer * metadata["xpixelsize"]
                metadata["x2"] += idx_buffer * metadata["xpixelsize"]

            elif orig_dim_y < new_dim:
                idx_buffer = int((new_dim - orig_dim_y) / 2.0)
                R_[:, :, idx_buffer : (idx_buffer + orig_dim_y), :] = R
                metadata["y1"] -= idx_buffer * metadata["ypixelsize"]
                metadata["y2"] += idx_buffer * metadata["ypixelsize"]

        elif method == "crop":
            new_dim = np.min(orig_dim[2:])
            R_ = np.zeros((orig_dim_n, orig_dim_t, new_dim, new_dim))

            if orig_dim_x > new_dim:
                idx_buffer = int((orig_dim_x - new_dim) / 2.0)
                R_ = R[:, :, :, idx_buffer : (idx_buffer + new_dim)]
                metadata["x1"] += idx_buffer * metadata["xpixelsize"]
                metadata["x2"] -= idx_buffer * metadata["xpixelsize"]

            elif orig_dim_y > new_dim:
                idx_buffer = int((orig_dim_y - new_dim) / 2.0)
                R_ = R[:, :, idx_buffer : (idx_buffer + new_dim), :]
                metadata["y1"] += idx_buffer * metadata["ypixelsize"]
                metadata["y2"] -= idx_buffer * metadata["ypixelsize"]

        else:
            raise ValueError("Unknown type")

        metadata["orig_domain"] = (orig_dim_y, orig_dim_x)
        metadata["square_method"] = method

        R_shape[-2] = R_.shape[-2]
        R_shape[-1] = R_.shape[-1]

        return R_.reshape(R_shape), metadata

    elif inverse:
        if len(R.shape) < 2:
            raise ValueError("The number of dimension must be > 2")
        if len(R.shape) == 2:
            R = R[None, None, :]
        if len(R.shape) == 3:
            R = R[None, :]
        if len(R.shape) > 4:
            raise ValueError("The number of dimension must be <= 4")

        method = metadata.pop("square_method")
        shape = metadata.pop("orig_domain")

        if R.shape[2] == shape[0] and R.shape[3] == shape[1]:
            return R.squeeze(), metadata

        R_ = np.zeros((R.shape[0], R.shape[1], shape[0], shape[1]))

        if method == "pad":
            if R.shape[2] == shape[0]:
                idx_buffer = int((R.shape[3] - shape[1]) / 2.0)
                R_ = R[:, :, :, idx_buffer : (idx_buffer + shape[1])]
                metadata["x1"] += idx_buffer * metadata["xpixelsize"]
                metadata["x2"] -= idx_buffer * metadata["xpixelsize"]

            elif R.shape[3] == shape[1]:
                idx_buffer = int((R.shape[2] - shape[0]) / 2.0)
                R_ = R[:, :, idx_buffer : (idx_buffer + shape[0]), :]
                metadata["y1"] += idx_buffer * metadata["ypixelsize"]
                metadata["y2"] -= idx_buffer * metadata["ypixelsize"]

        elif method == "crop":
            if R.shape[2] == shape[0]:
                idx_buffer = int((shape[1] - R.shape[3]) / 2.0)
                R_[:, :, :, idx_buffer : (idx_buffer + R.shape[3])] = R
                metadata["x1"] -= idx_buffer * metadata["xpixelsize"]
                metadata["x2"] += idx_buffer * metadata["xpixelsize"]

            elif R.shape[3] == shape[1]:
                idx_buffer = int((shape[0] - R.shape[2]) / 2.0)
                R_[:, :, idx_buffer : (idx_buffer + R.shape[2]), :] = R
                metadata["y1"] -= idx_buffer * metadata["ypixelsize"]
                metadata["y2"] += idx_buffer * metadata["ypixelsize"]

        R_shape[-2] = R_.shape[-2]
        R_shape[-1] = R_.shape[-1]

        return R_.reshape(R_shape), metadata
Either pad or crop a field to obtain a square domain.

Parameters
----------
R: array-like
Array of shape (m,n) or (t,m,n) containing the input fields.
metadata: dict
Metadata dictionary containing the x1, x2, y1, y2,
xpixelsize, ypixelsize,
attributes as described in the documentation of
:py:mod:`pysteps.io.importers`.
method: {'pad', 'crop'}, optional
Either pad or crop.
If pad, an equal number of zeros is added to both ends of its shortest
side in order to produce a square domain.
If crop, an equal number of pixels is removed
to both ends of its longest side in order to produce a square domain.
Note that the crop method involves an irreversible loss of data.
inverse: bool, optional
Perform the inverse method to recover the original domain shape.
After a crop, the inverse is performed by padding the field with zeros.

Returns
-------
R: array-like
the reshape dataset
metadata: dict
the metadata with updated attributes.
pysteps/utils/dimension.py