polyline_mapslice fails for mapformat1

[veenstrajelmer]

When applying polyline_mapslice to a uds with mapformat1, it raises KeyError: 'Mesh2D_s1'. This is probably only relevant for 2D models since it finds the wl/bl variables there. Since mapformat1 does not have the gridname as varname prefix, this fails.

MWE:

import dfm_tools as dfmt
import matplotlib.pyplot as plt
plt.close('all')
import numpy as np

file_nc = r"p:\1230882-emodnet_hrsm\GTSMv5.0\runs\_testcases_mapformat_ncformat\GM43_testmap_mf1_ncf3\output\gtsm_model_0004_map.nc"
uds = dfmt.open_partitioned_dataset(file_nc)

uds_sel = uds.ugrid.sel(x=slice(None,45), y=slice(20,None)).isel(time=-1)
fig,ax = plt.subplots()
uds_sel.s1.ugrid.plot(ax=ax)
line_array = np.array([[33.93558625, 25.31932949],
       [37.7394547 , 25.99706362]])
mapslice = dfmt.polyline_mapslice(uds_sel, line_array)

Can be solved by editing get_xzcoords_onintersection, but then somewhat neater:

        varn_wl = f'{gridname}_s1' #mapformat4
        varn_bl = f'{gridname}_flowelem_bl' #mapformat4
        if (varn_wl not in uds.data_vars) & (varn_bl not in uds.data_vars): #if "nNetElem" in uds.dims:
            varn_wl = 's1' #mapformat1
            varn_bl = 'FlowElem_bl' #mapformat1
        data_frommap_wl3_sel = uds_sel[varn_wl].to_numpy()
        data_frommap_bl_sel = uds_sel[varn_bl].to_numpy()

It also happens for 3D models since it fails to find the layer/interface dimensions, they are not set in the topology attributes:

import matplotlib.pyplot as plt
plt.close("all")
import numpy as np
import dfm_tools as dfmt

map_netcdf = r'p:\archivedprojects\1210301-storten-diepe-delen\03_modelling\02_simulations_FM\2018\02_simulations\PvH_a19refDFMsmallREVps_wlnorm_hw_sigr_coarse2_007f02\DFM_OUTPUT_g34i_2d\g34i_2d_0000_map.nc'
uds = dfmt.open_partitioned_dataset(map_netcdf)
uds_sel = uds.isel(time=-1)
fig,ax = plt.subplots()
uds_sel.s1.ugrid.plot()
line_array = np.array([[ 74737.27889613, 379753.94930325],
       [ 71240.55764847, 375812.95795849]])
xr_crs_ugrid = dfmt.polyline_mapslice(uds_sel, line_array)
fig,ax = plt.subplots()
xr_crs_ugrid.s1.ugrid.plot()

Therefore, a 3D dataset is incorrectly interpreted as being 2D. Fix in get_vertical_dimensions() by adding:

    elif "laydim" in uds.dims and "wdim" in uds.dims:
        # workaround for mapformat1
        layer_dimension = "laydim"
        interface_dimension = "wdim"
        return layer_dimension, interface_dimension

This in turn raises KeyError: 'layers present, but unknown layertype, expected one of variables: mesh2d_flowelem_zw, mesh2d_layer_sigma, mesh2d_layer_z' since the layertype is not derived in a generic way. Renaming these as a workaround raises KeyError: 'depth' since that key is not present in the formula_terms. A link to FlowElem_bl should be added here, but keep in mind the sign is opposite as expected so add the correct standard_name from reconstruct_zw_zcc_fromsigma. Seems best to cluster all mapformat1 preprocessing in a separate preprocessing function, to start with:

if set(["laydim","wdim"]).issubset(ds.dims):
    # set layer/interface dimension attrs in mesh_topology variable
    varn_mesh_topology = list(ds.filter_by_attrs(cf_role="mesh_topology").data_vars)[0]
    vertical_attrs = {"layer_dimension":"laydim",
                      "interface_dimension":"wdim",
                      }
    ds[varn_mesh_topology] = ds[varn_mesh_topology].assign_attrs(vertical_attrs)
if set(["LayCoord_cc","LayCoord_w"]).issubset(ds.data_vars):
    # rename layer variables so they can be recognized as layer/interface formula_terms variables
    ds = ds.rename_vars(LayCoord_cc="LayCoord_cc_layer", LayCoord_w="LayCoord_w_interface")

Keep in mind that the renamed varnames should also be updated in formula_terms attrs.

Todo:

• Generalize finding of sigma layer variables #764
• implement further fix
• add local unittest for 2D and 3D
• update whatsnew