Deltares / xugrid

Xarray and unstructured grids
https://deltares.github.io/xugrid/
MIT License
53 stars 8 forks source link

plotting of fine grained meshes is sub-optimal because of edge drawing #203

Open veenstrajelmer opened 5 months ago

veenstrajelmer commented 5 months ago

When plotting a ugrid with large resolution, sometimes narrow branches with fine cells are drawn too wide because the edges are drawn. I have created a minimal reproducible example with one of the datasets included in xugrid:

import matplotlib.pyplot as plt
plt.close("all")
import xugrid as xu

uds = xu.data.adh_san_diego()
uds_sel = uds.ugrid.sel(x=slice(470000, 490000), y=slice(3605000, None))

vmin, vmax = 0, 10

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10,7), sharex=True, sharey=True)

uds_sel.elevation.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax1) #edgecolors="face", linewidths=None
uds_sel.elevation.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax2, edgecolors="face", linewidths=None) #defaults
uds_sel.elevation.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax3, edgecolors="face", linewidths=0)
uds_sel.elevation.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax4, edgecolors="none", linewidths=0)

for ax in [ax1, ax2, ax3, ax4]:
    ax.scatter(478100, 3616500, s=300, edgecolors="r", facecolors='none')

fig.tight_layout()

This is not the best example. However, it is clear that the empty area between cells is wider in some cases than in others. The lower-left example shows the most whitespace but also has white edges all over the grid: image

A better example would be with:

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

file_nc = r"p:\11209585-dtp\voor_Caroline\01_scripts\post\python_differences\225deg-Car004-cropped.nc"
uds = dfmt.open_partitioned_dataset(file_nc)
uds_sel = uds.isel(time=-1).ugrid.sel(x=slice(3, None), y=slice(None, 53))

vmin, vmax = 0, 10

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10,7), sharex=True, sharey=True)

uds_sel.mesh2d_s1.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax1) #edgecolors="face", linewidths=None
uds_sel.mesh2d_s1.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax2, edgecolors="face", linewidths=None) #defaults
uds_sel.mesh2d_s1.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax3, edgecolors="face", linewidths=0)
uds_sel.mesh2d_s1.ugrid.plot(vmin=vmin, vmax=vmax, ax=ax4, edgecolors="none", linewidths=0)

fig.tight_layout()

image

Huite commented 5 months ago

Since you mention it's a better example, maybe you could post the pictures for the file nonetheless? (The p:\11209585-dtp\voor_Caroline\01_scripts\post\python_differences\225deg-Car004-cropped.nc I mean)

Having the MWE is helpful, but those pictures probably clarify it. I atleast have trouble seeing what the problem exactly is in these pictures.

Also, what's the question/issue exactly? Do you suggest a different default? Or do none of the option above look satisfactory?

veenstrajelmer commented 5 months ago

@Huite the mwe and figure are now also added to the issue description for the other example. Unfortunately none of the alternatives result in a perfect plot. Furthermore, I think this is actually a matplotlib issue, but since most matplotlib applications do not draw these tiny polygons I guess we are the only ones that are affected by it.

Huite commented 5 months ago

Ah okay I think I understand: the rivers look really fat.

If I blow up the dpi:

fig, ax = plt.subplots(figsize=(10, 7), dpi=600)
uds_sel.mesh2d_s1.ugrid.plot(vmin=vmin, vmax=vmax, edgecolors="none")

The rivers are drawn accurately: image

But you also get these edge artifacts: image

If you do edgecolors="face", it smears like crazy:

image

Huite commented 5 months ago

So in essence: you either have too small polygons, or too thick lines.

You want the lines to fill up the holes between the polygons, but not so much that it smears everything. The easiest way to achieve this seems to fiddle with the linewidths parameter.

It's a bit hard to find what linewidths means exactly, but apparently it's points: https://discourse.matplotlib.org/t/what-are-the-units-for-linewidths-parameter/22982/3

Where 1 point is 0.353 mm or 1/72 of an inch: https://en.wikipedia.org/wiki/Point_(typography)

fig, ax = plt.subplots(figsize=(7, 7), dpi=300)
uds_sel.mesh2d_s1.ugrid.plot(vmin=vmin, vmax=vmax, linewidths=0.1)

At 300 dpi, 0.1 is not quite enough; if you zoom in: image

But 0.3 seems to do the trick: image

0.3 seems to work fairly well across the board here, also at the default 100 dpi: image

Total extent is about 200 km in the y direction, with 7 inches times 72 is 514 points. A single point is then approximately 400 m on the map. So we'd expect the rivers to be buffered 200 m at both sides, increasing the total width by ~400 m.

This analysis doesn't help that much though, because I wouldn't know how to compute the widths of the gaps in the North sea between the polygons...

So probably best to fiddle with the linewidths: get them a small as possible, without the gaps appearing.

veenstrajelmer commented 5 months ago

Thanks for these efforts. However, it is a bit cumbersome to fiddle with linewidths, since they will differ for each model output file and zoom level right? For instance in the adh_san_diego plot this still gives white lines with linewidth=0.3. My preference would be that matplotlib does not draw lines/edges when using edgecolors="none" or linewidth=0. Furthermore, the faces are perfectly adjacent, so they should be drawn without white spaces in between I would say.

Huite commented 5 months ago

I agree, but it's a somewhat workable solution for now.

I tried switching to other backends, but the white lines are drawn consistently.

It's a matplotlib issue, so you could try raising an issue there, or ask on the matplotlib discourse forum first: https://discourse.matplotlib.org/

It might be partly a stylistic choice. I think I recall seeing it on other patches as well, so I'm not entirely sure it's just the PolyCollection, it could be all patches. I don't have my hopes up that they'll be eager to change this though, I reckon it has been like this forever.

As you mentioned, it's partly caused by having so many small elements, which is apparently triggering it. This also means it's a little slow at drawing the plot: I'm pretty sure it's looping over every polygon while rasterizing to pixels. It's been on my list to add some datashader trimesh support: https://datashader.org/user_guide/Trimesh.html In that case, datashader is doing the rasterization (via numba). It doesn't draw the edges, and then just class matplotlib to do an imshow (if I recall correctly).

Huite commented 5 months ago

From matplotlib: here's a good explanation https://github.com/matplotlib/matplotlib/issues/9574#issuecomment-339532262

This is a fairly difficult to solve problem in general. See e.g. https://computergraphics.stackexchange.com/questions/1824/what-is-illustrators-vector-rasterization-process.

It's not, strictly speaking, because the edges of the polygons are not matching. It's because the renderer only sees one patch at a time. When it antialiases the pixels on the edges of the polygon, the only thing it can do is assign partial occupancy by using partial transparency. When the renderer sees the second polygon, it will do the same. So an edge pixel which was supposed to have, say, left 50% red and right 50% green gets rendered as (white background) (50% transparency red) (50% transparency green), and thus the white "leaks through". This is known as a conflation artefact (I think due to the "conflation" between coverage and alpha).

AFAIK the only way to solve this is to emit your polygons as mesh objects, so that the renderer actually "knows" that the two edges are supposed to be the same (and even then you need a renderer that actually knows how to use that information). See also #7841.