prepro_private.py

import numpy as np
from netCDF4 import Dataset
import urllib.request
from datetime import datetime, timedelta
import slimPre
import os

try:
    import cdsapi

    _have_cdsapi = True
except:
    _have_cdsapi = False

fmt = "%Y-%m-%d %H:%M:%S"


# --- FUNCTIONS TO FILL MASKED VALUES --- #
def fill_mask(tab, N=5):
    ma = tab.copy()
    if ma[0].mask.ndim != 2:
        return ma
    maskx, masky = np.where(ma[0, :, :].mask)
    for ipass in range(N):
        if maskx.shape[0] == 0:
            break
        # nok = number of neighbors that are not masked
        # isok = boolean stating whether neighbor is not masked
        nok = np.zeros((maskx.shape[0],), np.int32)
        v = np.zeros((ma.shape[0], maskx.shape[0]), np.float64)
        # loop through neighbors (diagonals could be added)
        for shiftx, shifty in ((-1, 0), (1, 0), (0, -1), (0, 1)):
            nx = maskx + shiftx
            ny = masky + shifty
            # check that neighbours are inside the table
            isok = (nx >= 0) & (nx < ma.shape[1]) & (ny >= 0) & (ny < ma.shape[2])
            # check that neighbours are not masked
            isok[isok] = ma[0, nx[isok], ny[isok]] != np.ma.masked
            v[:, isok] += ma[:, nx[isok], ny[isok]]
            nok += isok
        # mean of "ok" neighbors
        ma[:, maskx[nok > 0], masky[nok > 0]] = v[:, nok > 0] / nok[nok > 0]
        # we keep going with entries that did not have 'ok' neighbors
        maskx = maskx[nok == 0]
        masky = masky[nok == 0]
    return ma


def fill_mask3(tab):
    """Fill masked entry in 4D tab[time, depth, x, y] based on nearest non-masked entries"""
    val = tab.copy()
    # Horizontal first
    for d in range(val.shape[1]):
        val[:, d, :, :] = fill_mask(tab[:, d, :, :])
    # Vertical afterwards
    for i in range(val.shape[2]):
        for j in range(val.shape[3]):
            if val[0, 0, i, j] is np.ma.masked or val[0, -1, i, j] is not np.ma.masked:
                pass
            first_raw = val[0, :, i, j]
            ind_last_unmask = len(np.where(~first_raw.mask)[0])
            for idx in range(ind_last_unmask, val.shape[1]):
                val[:, idx, i, j] = val[:, ind_last_unmask - 1, i, j]
    return val


# --- USEFUL NETCDF FUNCTIONS --- #
def find_nc_variables(nc_file):
    time_var = None
    lon_var = None
    lat_var = None
    depth_var = None
    with Dataset(nc_file, "r") as f:
        for var_name, var in f.variables.items():
            if hasattr(var, "standard_name"):
                if var.standard_name == "longitude":
                    lon_var = var_name
                elif var.standard_name == "latitude":
                    lat_var = var_name
                elif var.standard_name == "time":
                    time_var = var_name
                elif var.standard_name == "depth":
                    depth_var = var_name
            elif hasattr(var, "long_name"):
                if var.long_name.lower() == "longitude":
                    lon_var = var_name
                elif var.long_name.lower() == "latitude":
                    lat_var = var_name
                elif (
                    var.long_name.lower() == "time"
                    or var.long_name.lower() == "valid time"
                ):
                    time_var = var_name
                elif var.long_name.lower() == "depth":
                    depth_var = var_name
            else:
                if var_name[:3].lower() == "lon":
                    lon_var = var_name
                elif var_name[:3].lower() == "lat":
                    lat_var = var_name
                elif var_name.lower() == "depth":
                    depth_var = var_name
                elif var_name.lower() == "time":
                    time_var = var_name
    if lon_var is None or lat_var is None:
        raise Exception("Data is not spatially dependent !!!")
    return time_var, lon_var, lat_var, depth_var


def merge_ncfiles(file_list, variables, out_file):
    if isinstance(variables, str):
        variables = [variables]
    attr_dict, grp_attrs = {}, {}
    time_var, lon_var, lat_var, depth_var = find_nc_variables(file_list[0])
    with Dataset(file_list[0], "r") as f:
        for attr_name in f.ncattrs():
            grp_attrs[attr_name] = f.getncattr(attr_name)
        for var_name, var in f.variables.items():
            if var_name in variables:
                ndims = len(var.get_dims())
                _3d = ndims == 4
            tmp = {}
            for attr_name in var.ncattrs():
                tmp[attr_name] = var.getncattr(attr_name)
            attr_dict[var_name] = tmp
        lat = np.array(f.variables[lat_var][:])
        lon = np.array(f.variables[lon_var][:])
        if depth_var is not None:
            depth = np.array(f.variables[depth_var])

    dim_0 = (0, depth.size, lat.size, lon.size) if _3d else (0, lat.size, lon.size)
    val, _scaled = {}, {}
    for var_name in variables:
        _scaled[var_name] = "scale_factor" in attr_dict[var_name].keys()
        np_type = np.int if _scaled[var_name] else np.float32
        val[var_name] = np.empty(dim_0, dtype=np_type)

    time, tc = np.empty(0), -1e9
    for file_name in file_list:
        with Dataset(file_name, "r") as f:
            tf = np.array(f.variables[time_var][:])
            idx = tf > tc
            for var_name, arr in val.items():
                if _scaled[var_name]:
                    np_type, sf = np.int, attr_dict[var_name]["scale_factor"]
                    offset = attr_dict[var_name]["add_offset"]
                else:
                    np_type, sf = np.float32, 1
                    offset = 0
                nc_var = f.variables[var_name]
                ncarr = ((np.array(nc_var[:]) - offset) / sf).astype(np_type)
                if np.ma.array(nc_var[0]).mask.ndim > 1 and _scaled[var_name]:
                    if _3d:
                        maskx, masky, maskz = np.where(np.ma.array(nc_var[0]).mask)
                        ncarr[:, maskx, masky, maskz] = attr_dict[var_name][
                            "_FillValue"
                        ]
                    else:
                        maskx, masky = np.where(np.ma.array(nc_var[0]).mask)
                        ncarr[:, maskx, masky] = attr_dict[var_name]["_FillValue"]
                val[var_name] = np.row_stack([arr, ncarr[idx, :]])
        time = np.append(time, tf[idx])
        tc = time[-1]

    if "valid_min" in attr_dict[time_var].keys():
        attr_dict[time_var]["valid_min"] = time.min()
        attr_dict[time_var]["valid_max"] = time.max()

    var_dict = {}
    with Dataset(out_file, "w") as f:
        for attr_name, attr_value in grp_attrs.items():
            setattr(f, attr_name, attr_value)
        f.createDimension(lon_var, lon.size)
        f.createDimension(lat_var, lat.size)
        f.createDimension(time_var, time.size)
        var_dict[time_var] = f.createVariable(time_var, "f8", (time_var,))
        var_dict[time_var][:] = time[:]
        var_dict[lon_var] = f.createVariable(lon_var, "f8", (lon_var,))
        var_dict[lon_var][:] = lon[:]
        var_dict[lat_var] = f.createVariable(lat_var, "f8", (lat_var,))
        var_dict[lat_var][:] = lat[:]
        if depth_var is not None:
            f.createDimension(depth_var, depth.size)
            var_dict[depth_var] = f.createVariable(depth_var, "f8", (depth_var,))
            var_dict[depth_var][:] = depth[:]

        dim = (
            (time_var, depth_var, lat_var, lon_var)
            if _3d
            else (time_var, lat_var, lon_var)
        )
        for var_name, arr in val.items():
            if _scaled[var_name]:
                nc_type, sf = "i2", attr_dict[var_name]["scale_factor"]
            else:
                nc_type, sf = "f4", 1
            if "_FillValue" in attr_dict[var_name].keys():
                var_dict[var_name] = f.createVariable(
                    var_name, nc_type, dim, fill_value=attr_dict[var_name]["_FillValue"]
                )
            else:
                var_dict[var_name] = f.createVariable(var_name, nc_type, dim)
            var_dict[var_name][:] = arr[:]
        for var_name in attr_dict.keys():
            for attr_name, attr_value in attr_dict[var_name].items():
                if attr_name not in ["_FillValue", "valid_range"]:
                    setattr(var_dict[var_name], attr_name, attr_value)


# --- FUNCTIONS TO DOWNLOAD ERA, MERCATOR AND EREEFS FORCINGS --- #
def cds_download(request, output_file):
    if _have_cdsapi:
        c = cdsapi.Client()
        c.retrieve("reanalysis-era5-single-levels", request, output_file)
    else:
        print(
            "Warning : cdsapi not found -> wind and atmospheric pressure cannot be downloaded :-("
        )


def download_wind(
    lon_min, lon_max, lat_min, lat_max, initial_time, final_time, path_output
):

    d0 = datetime.strptime(initial_time, fmt)
    d1 = datetime.strptime(final_time, fmt)

    request_dict = {
        "product_type": "reanalysis",
        "format": "netcdf",
        "variable": [
            "10m_u_component_of_wind",
            "10m_v_component_of_wind",
            "mean_sea_level_pressure",
        ],
        "time": ["%02i:00" % i for i in range(24)],
        "area": [lat_max + 1, lon_min - 1, lat_min - 1, lon_max + 1],
    }

    year, month_from = d0.year, d0.month
    if d0.month == d1.month and d0.year == d1.year:
        request_dict["day"] = [i for i in range(d0.day, d1.day + 1)]
    else:
        request_dict["day"] = [i for i in range(1, 32)]

    file_list = []
    while year <= d1.year:
        request_dict["year"] = year
        month_to = d1.month if year == d1.year else 12
        request_dict["month"] = [i for i in range(month_from, month_to + 1)]
        filename = path_output + "/wind_tmp_%i.nc" % (year - d0.year)
        file_list += [filename]
        cds_download(request_dict, filename)
        month_from = 1
        year += 1
    outfile = path_output + "/wind_%s_%s.nc" % (
        d0.strftime("%Y%m%d"),
        d1.strftime("%Y%m%d"),
    )

    if len(file_list) == 1:
        os.rename(file_list[0], outfile)
    else:
        merge_ncfiles(file_list, ["msl", "u10", "v10"], outfile)
    os.system("rm " + path_output + "/wind_tmp*")


def download_currents(
    lon_min,
    lon_max,
    lat_min,
    lat_max,
    initial_time,
    final_time,
    user,
    password,
    path_output,
):
    d = datetime.strptime("2018-12-25 00:00:00", fmt)
    di = datetime.strptime(initial_time, fmt)
    df = datetime.strptime(final_time, fmt)
    variables = ["zos", "uo", "vo", "so", "thetao"]
    file_name_fmt = (
        "mercator_%s_" + di.strftime("%Y%m%d") + "_" + df.strftime("%Y%m%d") + ".nc"
    )
    #    fmt_da = "mercator_%s_"+di.strftime("%Y%m%d")+"_"+ df.strftime("%Y%m%d")+"_DA.nc"
    delta = timedelta(days=1)
    if df > d:
        motu_address = "http://nrt.cmems-du.eu/motu-web/Motu"
        service_id = "GLOBAL_ANALYSIS_FORECAST_PHY_001_024-TDS"
        product_id = "global-analysis-forecast-phy-001-024"
    else:
        motu_address = "http://my.cmems-du.eu/motu-web/Motu"
        service_id = "GLOBAL_REANALYSIS_PHY_001_030-TDS"
        product_id = "global-reanalysis-phy-001-030-daily"
    request = (
        "python3 -m motuclient --motu %s --service-id %s --product-id %s --longitude-min %d --longitude-max %d --latitude-min %d --latitude-max %d --date-min '%s' --date-max '%s' --depth-min 0.0 --depth-max 6000.0 --variable %s --out-dir '%s' --out-name '%s' --user "
        + user
        + " --pwd "
        + password
    )

    downloaded_files = {}
    for var in variables:
        file_name = file_name_fmt % (var)
        downloaded_files[var] = path_output + "/" + file_name
        os.system(
            request
            % (
                motu_address,
                service_id,
                product_id,
                np.floor(lon_min - 1),
                np.ceil(lon_max + 1),
                np.floor(lat_min - 1),
                np.ceil(lat_max + 1),
                (di - delta).strftime(fmt),
                (df + delta).strftime(fmt),
                var,
                path_output,
                file_name,
            )
        )
        if var in ["uo", "vo"]:
            os.system(
                "ncwa -a depth %s %s"
                % (
                    path_output + "/" + file_name,
                    path_output + "/" + file_name[:-3] + "_DA.nc",
                )
            )


def download_tpxo(path_tpxo_h, path_tpxo_u, user, password):
    if not os.path.isfile(path_tpxo_h):
        url = "ftp://geo07.elie.ucl.ac.be//export/miro/vvallaeys/slim_data/tides/h_tpxo9.nc"
        urllib.request.urlretrieve(url, path_tpxo_h)
    if not os.path.isfile(path_tpxo_u):
        url = "ftp://geo07.elie.ucl.ac.be//export/miro/vvallaeys/slim_data/tides/u_tpxo9.nc"
        urllib.request.urlretrieve(url, path_tpxo_u)


def write_elapsed_time(seconds):
    m, s = divmod(seconds, 60)
    return "%3d min %02.3f sec" % (m, s)


# --- FUNCTION TO COMPUTE TIDES COMPONENTS --- #
def getTPXO9(lon_min, lon_max, lat_min, lat_max, full_H, full_U, path_output):
    if lon_min * lon_max < 0.0:
        print("Error : not working around Greenwich now !")
        exit(-1)

    # Info nx : j 0 10800 lon 0 180 puis -180 à 0
    # Info ny : i 0 5400 lat -90 90

    y_min = np.floor((lat_min + 90) * 30).astype("i8")
    y_max = np.ceil((lat_max + 90) * 30).astype("i8")
    if lon_min > 0.0:
        x_min = np.floor(lon_min * 30).astype("i8")
        x_max = np.ceil(lon_max * 30).astype("i8")
    else:
        x_min = np.floor((lon_min + 180) * 30).astype("i8") + 5400
        x_max = np.ceil((lon_max + 180) * 30).astype("i8") + 5400

    tmp_H = "h_tpxo9_raw.nc"
    tmp_U = "u_tpxo9_raw.nc"

    print("Step 1/8 : Cut the local area from global TPXO 9")
    os.system(
        "ncks -d nx,"
        + str(x_min)
        + ","
        + str(x_max)
        + " -d ny,"
        + str(y_min)
        + ","
        + str(y_max)
        + " "
        + full_H
        + " "
        + tmp_H
    )
    os.system(
        "ncks -d nx,"
        + str(x_min)
        + ","
        + str(x_max)
        + " -d ny,"
        + str(y_min)
        + ","
        + str(y_max)
        + " "
        + full_U
        + " "
        + tmp_U
    )

    with Dataset(tmp_H, "r") as topex:
        ha = np.ma.array(topex.variables["ha"][:])
        hp = np.ma.array(topex.variables["hp"][:])
        x = np.array(topex.variables["lon_z"][:])
        y = np.array(topex.variables["lat_z"][:])
        conc = np.array(topex.variables["con"][:])

    nx = x.shape[0]
    ny = y.shape[1]
    nc = ha.shape[0]
    nt = conc.shape[1]

    ha = np.ma.masked_where(ha <= 1e-6, ha)
    hp = np.ma.masked_where(np.abs(hp) <= 1e-6, hp)

    print("Step 2/8 : Fix masked elevation amplitude")
    ha_new = fill_mask(ha)
    print("Step 3/8 : Fix masked elevation phase (cos)")
    cos = np.cos(hp * np.pi / 180.0)
    cos_new = fill_mask(cos)
    print("Step 4/8 : Fix masked elevation phase (sin)")
    sin = np.sin(hp * np.pi / 180.0)
    sin_new = fill_mask(sin)
    hp_new = np.arctan2(sin_new, cos_new) * 180.0 / np.pi

    with Dataset(path_output + "/h_tpxo9_zone.nc", "w") as new_topex:
        new_topex.createDimension("nx", nx)
        new_topex.createDimension("ny", ny)
        new_topex.createDimension("nc", nc)
        new_topex.createDimension("nt", nt)
        longitude = new_topex.createVariable("lon_z", "f4", ("nx", "ny"))
        latitude = new_topex.createVariable("lat_z", "f4", ("nx", "ny"))
        Ha = new_topex.createVariable(
            "ha", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        Hp = new_topex.createVariable(
            "hp", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        C = new_topex.createVariable("con", "c", ("nc", "nt"))

        longitude[:] = x
        latitude[:] = y
        Ha[:] = ha_new
        Hp[:] = hp_new
        C[:] = conc

    with Dataset(tmp_U, "r") as topex:
        up = np.ma.array(topex.variables["up"][:])
        vp = np.ma.array(topex.variables["vp"][:])
        ua = np.array(topex.variables["ua"][:])
        Ua = np.array(topex.variables["Ua"][:])
        va = np.array(topex.variables["va"][:])
        Va = np.array(topex.variables["Va"][:])
        xu = np.array(topex.variables["lon_u"][:])
        yu = np.array(topex.variables["lat_u"][:])
        xv = np.array(topex.variables["lon_v"][:])
        yv = np.array(topex.variables["lat_v"][:])
        conc = np.array(topex.variables["con"][:])

    nx = xu.shape[0]
    ny = yu.shape[1]
    nc = ua.shape[0]
    nt = conc.shape[1]

    up = np.ma.masked_where(np.abs(up) <= 1e-6, up)
    vp = np.ma.masked_where(np.abs(vp) <= 1e-6, vp)

    print("Step 5/8 : Fix masked eastward velocity phase (cos)")
    cos = np.cos(up * np.pi / 180.0)
    cos_new = fill_mask(cos)
    print("Step 6/8 : Fix masked eastward velocity phase (sin)")
    sin = np.sin(up * np.pi / 180.0)
    sin_new = fill_mask(sin)
    up_new = np.arctan2(sin_new, cos_new) * 180.0 / np.pi
    print("Step 7/8 : Fix masked northward velocity phase (cos)")
    cos = np.cos(vp * np.pi / 180.0)
    cos_new = fill_mask(cos)
    print("Step 8/8 : Fix masked northward velocity phase (sin)")
    sin = np.sin(vp * np.pi / 180.0)
    sin_new = fill_mask(sin)
    vp_new = np.arctan2(sin_new, cos_new) * 180.0 / np.pi

    with Dataset(path_output + "/u_tpxo9_zone.nc", "w") as new_topex:
        new_topex.createDimension("nx", nx)
        new_topex.createDimension("ny", ny)
        new_topex.createDimension("nc", nc)
        new_topex.createDimension("nt", nt)
        lonU = new_topex.createVariable("lon_u", "f4", ("nx", "ny"))
        latU = new_topex.createVariable("lat_u", "f4", ("nx", "ny"))
        lonV = new_topex.createVariable("lon_v", "f4", ("nx", "ny"))
        latV = new_topex.createVariable("lat_v", "f4", ("nx", "ny"))
        uA = new_topex.createVariable(
            "ua", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        UA = new_topex.createVariable(
            "Ua", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        uP = new_topex.createVariable(
            "up", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        vA = new_topex.createVariable(
            "va", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        VA = new_topex.createVariable(
            "Va", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        vP = new_topex.createVariable(
            "vp", "f8", ("nc", "nx", "ny"), fill_value=-9999.0
        )
        C = new_topex.createVariable("con", "c", ("nc", "nt"))

        lonU[:] = xu
        latU[:] = yu
        lonV[:] = xv
        latV[:] = yv
        uA[:] = ua
        UA[:] = Ua
        uP[:] = up_new
        vA[:] = va
        VA[:] = Va
        vP[:] = vp_new
        C[:] = conc

    os.system("rm " + tmp_H + " " + tmp_U)


# --- FETCH FORCINGS --- #
def fetch_forcings(
    lon_min,
    lon_max,
    lat_min,
    lat_max,
    initial_time,
    final_time,
    path_output,
    forcings=[],
    path_to_tpxo=None,
    path_to_gebco=None,
    user="ehanert",
    password="merc@tor4SLIM",
    on_geo_server=False,
):
    known_forcings = [
        "tides",
        "wind",
        "mercator",
    ]
    if isinstance(forcings, str):
        forcings = [forcings]
    if len(forcings) == 0:
        forcings = known_forcings
    else:
        unknown = []
        for i in range(len(forcings)):
            forcings[i] = forcings[i].lower()
            if forcings[i] not in known_forcings:
                unknown += [forcings[i]]
        if len(unknown) == len(forcings):
            forcings = known_forcings
            print(
                "Warning : No known forcing name found in argument 'forcings', dowloading all forcings by default ..."
            )
        elif len(unknown) > 0:
            print(
                "Warning : Unknown forcing names: "
                + ", ".join(["'%s'" % i for i in unknown])
            )
            print(
                "Info    : Known forcing names are: "
                + ", ".join(["'%s'" % i for i in known_forcings])
            )

    slimPre.make_directory(path_output)

    if "wind" in forcings:
        print("Info    : Importing wind...")
        download_wind(
            lon_min, lon_max, lat_min, lat_max, initial_time, final_time, path_output
        )

    if "tides" in forcings:
        print("Info    : Importing tides...")
        if on_geo_server:
            path_tpxo_h = "/export/miro/vvallaeys/slim_data/tides/h_tpxo9.nc"
            path_tpxo_u = "/export/miro/vvallaeys/slim_data/tides/u_tpxo9.nc"
        else:
            path_tpxo_h = (
                path_output if path_to_tpxo is None else path_to_tpxo
            ) + "/h_tpxo.nc"
            path_tpxo_u = (
                path_output if path_to_tpxo is None else path_to_tpxo
            ) + "/u_tpxo.nc"

        download_tpxo(path_tpxo_h, path_tpxo_u, user, password)
        getTPXO9(
            lon_min, lon_max, lat_min, lat_max, path_tpxo_h, path_tpxo_u, path_output
        )

    if "mercator" in forcings:
        # download mercator files
        print("Info    : Importing Mercator...")
        download_currents(
            lon_min,
            lon_max,
            lat_min,
            lat_max,
            initial_time,
            final_time,
            user,
            password,
            path_output,
        )
    print("fetching forcings -> done")


# --- INTERPOLATION OF FORCINGS ON MESH --- #
def convert_time(time, time_units):
    unit_name = time_units.split(" since ")[0]
    ref_date = time_units.split(" since ")[-1].split(".")[0]
    t0 = slimPre.slim_private._parse_time(ref_date)
    dt = None
    if unit_name == "days":
        dt = 86400.0
    elif unit_name == "hours":
        dt = 3600.0
    elif unit_name == "minutes":
        dt = 60.0
    elif unit_name == "seconds":
        dt = 1
    else:
        raise Exception("Unknown units: %s !!!" % unit_name)
    return t0 + dt * time


def interpolate_on_mesh(coords, time_obj, nc_file, var_name):
    """temporal and spatial interpolation on mesh

    arguments:
        * coord
            array of coordinates
        * time_obj
            slimPre.Time at which data must be interpolated
        * nc_file
            path to netcdf file containing the data to be treated for preprocessing
        * var_name
            string of the name of the variable to be preprocessed
    """
    # check coords
    coords = np.array(coords)
    if coords.ndim != 2:
        raise ValueError("Arguments 'coords' must have 2 dimensions !!!")
    elif coords.shape[0] == 0:
        return np.empty((time_obj._time.size, 0))
    elif coords.shape[1] == 2:
        x, y = coords[:, 0], coords[:, 1]
        z = None
    elif coords.shape[1] == 3:
        x, y, z = coords[:, 0], coords[:, 1], coords[:, 2]
    else:
        raise ValueError(
            "Wrong dimension for argument 'coords': it must have 2 (2d mesh) 3 (3d mesh) columns !!!"
        )
    time_var, lon_var, lat_var, depth_var = find_nc_variables(nc_file)

    # read file
    with Dataset(nc_file, "r") as f:
        u_arr = np.ma.array(f.variables[var_name][:])
        if time_var:
            time = np.array(f.variables[time_var][:])
            time_units = f.variables[time_var].units
            time[:] = convert_time(time[:], time_units)
        else:
            u_arr = u_arr[None, :]
        depth_dependent = depth_var is not None and u_arr.ndim == 4
        if not depth_dependent:
            u = fill_mask(u_arr)
        else:
            print("interpolation of depth-dependent variable")
            depths = np.array(f.variables[depth_var][:])
            u = fill_mask3(u_arr)
        lon = np.array(f.variables[lon_var][:])
        lat = np.array(f.variables[lat_var][:])

    if depth_dependent and z is None:
        raise ValueError("z coordinates must be given for 3d interpolation !!!")
    ox, dx = lon[0], lon[1] - lon[0]
    oy, dy = lat[0], lat[1] - lat[0]

    # some sanity checks...
    if x.max() < lon.min():
        x[:] += 360.0
    if (
        lon.min() > x.min()
        or lon.max() < x.max()
        or lat.min() > y.min()
        or lat.max() < y.max()
    ):
        raise Exception("data doesn't fully cover mesh geographic extent !!!")
    if time_var is None or time_obj is None:
        if not depth_dependent:
            return slimPre.interpolate_from_structured_grid(
                x, y, ox, oy, dx, dy, u[0], fill_value=0
            )
        else:
            return slimPre.interpolate_from_structured_grid_3D(
                x, y, -z, ox, oy, dx, dy, depths, u[0]
            )
    # some sanity checks again...
    time_vector = time_obj._time
    if time[0] > time_vector[0] or time[-1] < time_vector[-1]:
        raise Exception("time range of data doesn't fully cover time_vector !!!")
    # temporal  and spatial interpolation
    ur = np.empty((time_vector.size, x.size))
    for j in range(time_vector.size):
        i = time[time <= time_vector[j]].size - 1
        xi = (
            (time_vector[j] - time[i]) / (time[i + 1] - time[i])
            if i + 1 < time.size
            else 0
        )
        U = xi * u[(i + 1) % time.size, :] + (1 - xi) * u[i]
        if not depth_dependent:
            ur[j, :] = slimPre.interpolate_from_structured_grid(
                x, y, ox, oy, dx, dy, U, fill_value=0
            )
        else:
            ur[j, :] = slimPre.interpolate_from_structured_grid_3D(
                x, y, -z, ox, oy, dx, dy, depths, U
            )
    return ur