Dfsu: contourf plot type results error with NaN values

[mmfontana]

Describe the bug

I am trying to plot a. dfsu file with NaN values, which are physically justified. I am using "contourf" as the plot type to clearly show the isolines. However, I am getting the following error related to the NaN values in the mesh: "ValueError: z array must not contain non-finite values within the triangulation".

To Reproduce

Please find below how to reproduce the error. Also, some try to overcome the issue. The linear interpolation on the boundaries should not happen. I have tried to use some masks (out of mikeio), but this is still not the expected behavior - see below how MIKE Zero GUI deals with it ("Shaded contour").

.dfsu file

Expected

First try

Further try

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from mikeio import Dfsu

# Options
error = True
first_try = False
further_try = False

# File name
file_name = './time_of_arrival_summer.dfsu'
# Read file
dfs = Dfsu(file_name)
ds = dfs.read()

#### Error ####
if error:
    # Create axis
    fig, ax = plt.subplots()
    # Plot with error
    ds['min_time_of_arival'][0].plot(ax=ax,plot_type='contourf',levels=[0,1,2,3,4,5,6,7,8,9,10],cmap='rainbow')
    ax.set_ylim([7168010.85,7182717.85])
    ax.set_xlim([768708.83,796419.53])

#### First try ####
if first_try:
    # Create axis
    fig, ax = plt.subplots()
    # First try
    data = ds['min_time_of_arival'][0].values.copy()
    # Replace nan values
    ds['min_time_of_arival'][0] = np.where(np.isnan(data),10**-20,data)
    ds['min_time_of_arival'][0].plot(ax=ax,plot_type='contourf',levels=[0,1,2,3,4,5,6,7,8,9,10],cmap='rainbow')
    ax.set_ylim([7168010.85,7182717.85])
    ax.set_xlim([768708.83,796419.53])
    fig.suptitle('first_try')
    plt.savefig('first_try.png')    

#### Further try ####
if further_try:
    # Further try
    # Method originally from MIKEIO
    def create_tri_only_element_table(self,data=None,geometry=None):

        """Convert quad/tri mesh to pure tri-mesh"""
        if geometry is None:
            geometry = self

        ec = geometry.element_coordinates
        if geometry.is_tri_only:
            return np.vstack(geometry.element_table), ec, data

        if data is None:
            data = []

        elem_table = [
            list(geometry.element_table[i]) for i in range(geometry.n_elements)
        ]
        tmp_elmnt_nodes = elem_table.copy()
        for el, item in enumerate(tmp_elmnt_nodes):
            if len(item) == 4:
                elem_table.pop(el)  # remove quad element

                # insert two new tri elements in table
                elem_table.insert(el, item[:3])
                tri2_nodes = [item[i] for i in [2, 3, 0]]
                elem_table.append(tri2_nodes)

                # new center coordinates for new tri-elements
                ec[el] = geometry.node_coordinates[item[:3]].mean(axis=1)
                tri2_ec = geometry.node_coordinates[tri2_nodes].mean(axis=1)
                ec = np.append(ec, tri2_ec.reshape(1, -1), axis=0)

                # use same data in two new tri elements
                data = np.append(data, data[el])

        return np.asarray(elem_table), ec, data
    
    # Create triangulation
    z = np.squeeze(ds['min_time_of_arival'][0].values.copy())
    zn_raw = ds.geometry.get_node_centered_data(z)
    zn = np.where(np.isnan(zn_raw),10**-20,zn_raw)
    elem_table, ec, z = create_tri_only_element_table(dfs,data=data,geometry=dfs.geometry)
    nc = ds.geometry.node_coordinates
    triang = tri.Triangulation(nc[:,0],nc[:,1],elem_table)
    # Create axis
    fig, axs = plt.subplots(1,3,figsize=(20,5),sharey=True)
    # Masks
    isbad = np.isnan(zn_raw)
    masks = {'relaxed':np.all(np.where(isbad[triang.triangles],True,False),axis=1),
            'moderate':np.sum(np.where(isbad[triang.triangles],True,False),axis=1)>1,
            'restrictive':np.any(np.where(isbad[triang.triangles],True,False),axis=1)}
    # Plots
    i = 0
    # Loop over masks
    for mask in masks:
        # Different masks
        triang.set_mask(masks[mask])
        time_of_arrival_ax = axs[i].tricontourf(triang,zn,levels=[0,1,2,3,4,5,6,7,8,9,10],cmap='rainbow',extend='max')
        axs[i].set_ylim([7168010.85,7182717.85])
        axs[i].set_xlim([768708.83,796419.53])
        axs[i].set_title(mask)
        i += 1
    fig.colorbar(time_of_arrival_ax,ax=axs.ravel().tolist())
    fig.suptitle('further_try')
    plt.savefig('further_try.png')

System information:

• Python version: 3.8.10
• mikeio version: 1.0b2

[jsmariegaard]

Hi @mmfontana - sorry for the delay in getting back to you. We are really busy. I have not yet had the time to look into this but thanks to you, we are aware and will try to look at it soon. Thanks a lot for submitting this issue.

[mmfontana]

Hello @jsmariegaard. Thank you for your message and help. It is fine, I will plot through the GUI. But it would be nice to have it for the next Oil Spill projects (several scenarios). Thanks.

[jouArc]

Hello,
I'm also really interested in this bug correction when you have time because I use a lot mikeio and "contourf" in my work!
Many thanks! :)

[jsmariegaard]

Hello, I'm also really interested in this bug correction when you have time because I use a lot mikeio and "contourf" in my work! Many thanks! :)

Noted. Thanks for letting us know. As soon as we have time - we will look into this.

[alexander0042]

If anyone is looking for a temporary workaround, applying a filter using the elements isel argument to keep only non-zero elements gets around this for now! Something like:

dhiFilt = ds_list[m]['Current speed'].sel(time="2018-01-01 00:50:00").values > 0
dhiFiltIDX = [i for i, x in enumerate(dhiFilt) if x]

axDHI = ds_list[m]['Current speed'].sel(time="2018-01-01 00:50:00").isel(element=dhiFiltIDX).plot.contourf()

[jsmariegaard]

Thanks @alexander0042 - that's very helpful