Repository navigation
Added .plot.contour() function and added docs for visualizations - #1812
Conversation
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
dylannelson
added
documentation
ASV BenchmarkingBenchmark Comparison ResultsBenchmarks that have improved:
Benchmarks that have stayed the same:
Benchmarks that have got worse:
|
dylannelson
changed the title
.plot.contour() function and added docs for visualizations
|
Graphical review: Contours are looking good overall, so these are my fine-tuning preferences:
|
HoloViews and both plotting backends handle each path separately, so a noisy field with ~77k lines on a 786k-face mesh took 181 s to render with Matplotlib and 14.5 s with Bokeh, plus 2 s to build the element. With one path per level the lines are drawn as gaps at almost no cost: 0.14 s, 0.13 s and 0.07 s. The level is still the value of each path, so hover is unchanged, and the rendered images match. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Projecting every line with project_geometry took ~0.1 ms per line, or 7.4 s for the ~77k lines of a noisy field on a 786k-face mesh. A line with a projected point well inside the projection's domain has extent on the map, so all points are projected at once and only the remaining lines, such as those along the edge of the map or off it, are projected as lines. The lines kept are the same as before for 10 projections (global, hemispheric, regional and shifted central longitudes) on four datasets, and the check takes 0.04-1.1 s instead of 2.3-4.4 s on a 197k-face mesh. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Both were done with several NumPy calls per line, which for the ~77k lines of a noisy field on a 786k-face mesh took 1.1 s in _edge_contours and 0.5 s in _label_points. The lines are now handled as one array of points: _edge_contours unwraps them together and splits only the lines that cross the antimeridian (1.8 s -> 1.0 s including joining), and _label_points finds the middles only of the lines it labels (0.52 s -> 0.05 s). Longitudes are now unwrapped by whole multiples of 360, so a node on the antimeridian lands exactly on +-180 instead of a rounding error past it. That removes the zero-length pieces and repeated end points that clipping at the antimeridian gave for such lines. The lines are otherwise the same, and the labels are at the same points; lines of exactly equal length are now labeled in the order they were found. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
@cmdupuis3
|
|
@dylannelson (1) Yeah, I was thinking some kind of range, because I think it would enable more flexibility for showing certain granular features without overwhelming the plot. One other suggestion: usually for synoptic charts, you'll see the label carves out the contour behind it, so you can center the label on the contour without having the actual line clutter it up. That might actually be a use case for neighborhood filters; you can select where the labels will go, then delete contour line segments anywhere in a radius from them. |
_join_segments walked the segments in pure Python, which took 0.89 s for the ~650k segments of a noisy field on a 786k-face mesh. The kernel keeps the same walk, with the segments at each node in CSR arrays and a per-node cursor to the first unused one, and gives the same lines in the same order (checked on random graphs with branching nodes, self-loops and repeated segments, and on the segments of four meshes). It takes 0.074 s, and _edge_contours as a whole goes from 1.0 s to 0.2 s. The lines are returned as one array of nodes with offsets, the layout _edge_contours already built from the lists. The kernel holds no Python objects, so it is compiled with nogil=True and can run in threads (8 calls on 8 threads ran 4x faster than in sequence). Also applies ruff format to a line from the projection filter commit. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
(2) Actually, I'm thinking the tiny contour filter should be presented as an optional radius. Looking at your Maria contour plots, there are some small bits that I'm fairly sure are real over some of the Caribbean islands. I'm adding performance edits on cmd/contour-maps-opt, each optimization is committed separately. I'm seeing speedups on the order of 100x or more so far. Some of them require borrowing from cartopy, so you can merge from 32b29ec for the non-cartopy stuff and I can maybe add a PR for the rest if we want to hold off on that. |
rajeeja
left a comment
rajeeja
left a comment
There was a problem hiding this comment.
Requesting changes for the two correctness cases in the inline comments.
|
@dylannelson Yeah, 2 and 3 are more QOL fixes, they can probably be another PR |
Sevans711
left a comment
Sevans711
left a comment
There was a problem hiding this comment.
This looks really great overall, and I'm excited to use this new feature!
I have a few thoughts, mostly minor, and left some inline comments as well.
-
Re: summary of before and after, most cases look good to me, the only one I would push back on is the
temp.isel(time=[0], lev=[0])case. I would personally say this should crash with DimensionError, just liketempdoes, because the time and lev dimensions do exist, they just happen to have size 1. Crashing here would match xarray's behavior for specific plotting methods (likearr.plot.methodname()), which do not squeeze() automatically (see #1724 for more details). -
Re: the known concern about projection centered on 180, I agree it is out of scope for this PR; it is in fact already being tracked at #1795 (which has an open PR, so hopefully it will be fixed soon)!
-
(Non-blocking, could go into a followup PR) Is there any way to get the contour labels to be rotated to align with the local contour line?
-
(Non-blocking style question) would it make sense to add
%%timeat the top of some of the cells? E.g. this would follow the style of https://uxarray.readthedocs.io/en/latest/user-guide/mpl.html where %%time is placed at the top of cells relevant to explicit mentions of computation speed. -
Please add a bit more explanation for this sentence: "With
lat=the samples cover the whole latitude circle, sosteps=14401places one every 0.025° of longitude, about the width of a mesh cell". How did you derive / check these values? It's okay to not explain further the "about the width of a mesh cell" part if you just checked it by eye, but it would be helpful to at least show where the number 14401 comes from. -
(Non-blocking, can go into a later PR) Contours should probably infer the projection from the other plot object's projection if they are being combined via
*, or at least provide a clearer warning/help if the projections differ, but this can be pushed to a later PR. It feels like this would either be very easy to set up or very challenging to set up, but I don't know which it is, yet. -
Contours using a projection by themselves display strangely, see example below. Is there any way to fix this? It feels like a basic enough behavior that it would also be appropriate to expect a fix as part of this PR, but if it is especially challenging to fix I could also accept it getting pushed to a later PR.
Example of strange display when just showing contours, why is the plot on the left so squished? I would have expected the plot sizes to be the same.
import uxarray as ux
import cartopy.crs
arr = ux.tutorial.open_dataset('outCSne30-vortex')['psi']
proj = cartopy.crs.Robinson()
contours = arr.plot.contour(projection=proj, show_legend=False)
data_and_contours = arr.plot(projection=proj, cmap="viridis") * contours
contours + data_and_contours
| """Contour lines for data on an unstructured grid, without regridding. | ||
|
|
||
| Two methods are provided: | ||
|
|
||
| - ``"interpolated"``: the data is treated as values at the corners of a | ||
| triangulation of the grid (the grid's own faces for node-centered data, its | ||
| dual mesh for face-centered data) and contoured with linear interpolation | ||
| inside each triangle. | ||
| - ``"edges"``: a contour is made of the grid edges that separate faces above a | ||
| level from faces at or below it. No interpolation is done. | ||
| """ |
There was a problem hiding this comment.
This file has a lot of code, it makes me think you may already have in your mind a high-level overview of what problems are being solved in this file, but it isn't fully clear yet to me at a high level like this. Could you provide a few sentences or bullet points in the file docstring here about the main problems being solved / the difficulties with building contour maps?
For example, one thing is needing to figure out labels and where they go, another is needing to deal with lines crossing the antimeridian, another is needing to join segments into continuous lines (why?), and another is needing to triangulate faces (why?).
This doesn't need to be a long explanation, but even just one sentence about each problem would make it much easier to understand what this file is doing at a high level, so that in the future developers (or LLMs) don't need to parse the whole thing to understand its goals, they can use the docstring as a guide.
There was a problem hiding this comment.
I feel like those comments would belong to the actual methods that solve those issues though? My POV, I wouldn't even mention triangulation at all, it should just say "we're making contours with two methods" and leave any complexity out.
There was a problem hiding this comment.
I agree with Chris, the majority of the explanation I think should fall within each function. I can see how this is only a partial explanation though, from a previous version that only had a few functions, and so I removed it, as the other methods are plentiful now and have their own explanations.
Meanwhile constants.py exists in the same directory and it's literally only 2 lines of code:
# threshold for number of faces before recompensing rasterization
N_FACE_THRESHOLD = 10000matplotlib.py exists in the same directory and has 0 doc string/overview at all
accessor.py is even more broad and less self explanatory, and has no high level overview:
Not sure why we need a full explanation at the top, or what basis there is for needing one now.
There was a problem hiding this comment.
Linking as relevant: #1598. I believe the team eventually wants to see docstrings at the top of files, even if they are short? I would classify these example files as "missing a docstring", not as evidence to put no docstring at all. Even if the decision here ends up being to keep all the details out of the docstring, I would personally really appreciate at least keeping something like the first line you had originally, like """Contour lines for data on an unstructured grid, without regridding.""".
More details: the main reason I asked for more explanation here is the top-level complexity of this file; the other examples considered have much less top-level complexity.
- uxarray/plot/constants.py has only 2 lines, so it's not hard to quickly understand all of its contents.
- uxarray/plot/matplotlib.py has only 4 top-level methods/classes:
_ensure_dimensions,_RasterAxAttrs,_get_points_from_axis, and_nearest_neighbor_resample, so it's not too difficult to parse. - uxarray/plot/accessor.py similarly has only 4 top-level methods/classes:
_ensure_hvplot_imported,GridPlotAccessor,UxDataArrayPlotAccessor,UxDatasetPlotAccessor. - The proposed uxarray/plot.contour.py has 16 top-level methods:
_spatial_values,_wrap_longitude,_fan_corners,_ear_clipping,_fan_triangles,_triangulation,_split_at_antimeridian,_contour_levels,_interpolated_contours,_join_segments,_edge_contours,_drop_lines_outside_projection,_label_fraction,_label_points,_apply_options,_compute_contours. Some are called from other top-level methods, and some are called from outside the file. The tree of calls goes a few layers deep._compute_contourscalls_spatial_values,_contour_levels,_edge_contours,_interpolated_contours_edge_contourscalls_wrap_longitude,_join_segments,_split_at_antimeridian_interpolated_contourscalls_triangulation,_split_at_antimeridian_triangulationcalls_wrap_longitude,_fan_triangles_fan_trianglescalls_fan_cornersand_ear_clipping
There was a problem hiding this comment.
(Leaving a note about this for upcoming internal meeting, maybe it would be easiest to chat about it out loud and report back here with a decision?)
There was a problem hiding this comment.
I believe the team eventually wants to see docstrings at the top of files, even if they are short?
Honestly I think they're just a Claude-ism, I just leave them in if they're actually helpful, but otherwise I wouldn't expect one to be there at all.
There was a problem hiding this comment.
I believe the team eventually wants to see docstrings at the top of files, even if they are short?
Honestly I think they're just a Claude-ism, I just leave them in if they're actually helpful, but otherwise I wouldn't expect one to be there at all.
#1598 describes more about the benefits of including such docstrings, and has been triaged as To-Do. If there's disagreement about whether it should be done in general, could we move that discussion over there so that it doesn't get lost?
In this particular PR, though, I'd argue that """Contour lines for data on an unstructured grid, without regridding.""" is quite helpful, and has very little cost, it's just one line.
Upon looking at this again, I'm realizing part of the reason I wanted a high-level overview is that the docstrings for the 16 helper functions don't clearly specify their expected entry points and outputs. For example, if I isolate just one function and its docstring:
def _ear_clipping(x: np.ndarray, y: np.ndarray) -> list[tuple[int, int, int]]:
"""Triangles of a simple polygon, as positions of its corners."""It's hard to tell exactly what the intended inputs and outputs mean. I understand the rough concepts, but what shapes should x and y have? Why do you only need x and y to specify a polygon and not also z? How long with the output list be?
As an alternative to a detailed top-of-file docstring, I would also be happy to resolve this point without needing to discuss at internal meetings if the helper functions' docstrings were all updated to clearly answer questions like these. This would make it easy to understand each function as a standalone object, without needing a high-level map to understand how it relates to all 15 other functions in order to review the file.
There was a problem hiding this comment.
Reportback from internal discussion: please update the docstrings of helper functions to improve clarity, especially with regards to the shape and meaning of their inputs and outputs!
Meanwhile, there is no need to put any docstring at the top of the file. (My personal preference would still be to at least include the one-liner, but I this is a non-blocking style preference here.)
There was a problem hiding this comment.
Nearly all of the new internal function have descriptions that are 2-3x longer now and state how they are used
Ear clipping now reads:
def _ear_clipping(x: np.ndarray, y: np.ndarray) -> list[tuple[int, int, int]]:
"""Triangles of a simple polygon, as positions of its corners.
Used by ``_fan_triangles`` for the rare dented faces in which no corner can see
every other corner (a C-shaped face, for example), so that no fan of triangles
covers them. It cuts off one "ear" at a time: three corners in a row whose triangle
turns the same way as the polygon and holds no other corner. Plain Python is fast
enough, as it only runs on those faces."""
Lifts the buffers out of the loop in order to avoid tiny numpy arrays in the _fan_corners numba method
Commented on that line, I can see the desire for consistency, but don't understand why we need to follow xarray in this scenario since it seems like it wouldn't address the main concern of reducing the extra dims which cause problems.
Awesome, and yeah!
Probably! That would be an interesting addition. Maybe for the next update?
Not really sure where this would be useful since none of this is that computationally heavy, or likely to be looped heavily, and is 100% dependent on whatever CI pipeline bot runs it, so not really a consistent metric. But if you want to add it, you sure can, may just be high variance noise though.
360 degrees divided by .025 degree steps = 14400 (+1 for looping back around). But yeah that's a bit convoluted and not needed. Removed the explanation.
Agreed! If people want they can set a variable to a particular projection and use that. But yeah not sure if I can get one map object to inherit the projection from another. I'll look into that.
Ha, and I thought I caught all the weird display quirks. I didn't do many contours just on their own with additional maps to the side. Especially since I tended to put them with/on other data. Great catch, thanks. I'm working on a fix, but this one is weird, hope to test it further and upload soon. |
|
Thank you for looking into all of my comments, I really appreciate it! Replying below for each of the points from my main message (partly as a note to myself also for when I re-review).
|
Conflicts resolved: - uxarray/plot/contour.py: kept both _unwrap_lines and _grid_extents, and both _paths_by_level and _label_fraction. _label_points combines the line lengths found for all lines at once with labeling only a fraction of the levels. - uxarray/plot/accessor.py: one NaN-separated path per level, keyed by the fixed "level" dimension, with labels for a fraction of the levels. - test/test_plot.py: kept both the projection filter test and the frame test. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
Local timing data for cmd/contour-maps-opt: The gain is large wherever there are many contour lines, and close to nothing with only a few dozen. Times are
In the real-data rows, each pair is Bokeh / Matplotlib. The biggest gain is in rendering. One path per level instead of one per line takes the noisy-mesh render from 47–224 s to under 0.1 s without a projection. Some more optimizations are available on cmd/contour-maps-opt-2, but I'm deferring that to another PR due to possible licensing questions with cartopy. |
Merging contour-maps-opt into contour-maps
…numpy-in-numba avoid tiny numpy in contour._fan_corners
cmdupuis3
left a comment
cmdupuis3
left a comment
There was a problem hiding this comment.
One major issue, grids given in meters produce contours as if they were in degrees:
import warnings; warnings.simplefilter("ignore")
import numpy as np, xarray as xr, uxarray as ux
from uxarray.plot.contour import _compute_contours, _grid_extents
# 5 x 4 patch of 100 m quads in projected coordinates
X0, Y0 = 1_500_000.0, 800_000.0
xs, ys = np.meshgrid(X0 + 100.0 * np.arange(6), Y0 + 100.0 * np.arange(5))
faces = [[j*6+i, j*6+i+1, (j+1)*6+i+1, (j+1)*6+i] for j in range(4) for i in range(5)]
g = ux.Grid(xr.Dataset({
"node_lon": xr.DataArray(xs.ravel(), dims=["n_node"], attrs={"standard_name": "projection_x_coordinate", "units": "m"}),
"node_lat": xr.DataArray(ys.ravel(), dims=["n_node"], attrs={"standard_name": "projection_y_coordinate", "units": "m"}),
"face_node_connectivity": xr.DataArray(np.array(faces), dims=["n_face", "n_max_face_nodes"],
attrs={"cf_role": "face_node_connectivity", "start_index": 0, "_FillValue": -1}),
}), source_grid_spec="UGRID")
dist = ux.UxDataArray(xr.DataArray(g.node_lon.values - X0, dims=["n_node"]), uxgrid=g)
print(_compute_contours(dist, [250.0])[0][1][0]) # branch: x=130 expected: x=1500250
col = ux.UxDataArray(xr.DataArray(np.tile(np.arange(5.0), 4), dims=["n_face"]), uxgrid=g)
print([l[0, 0] for _, l in _compute_contours(col, [0.5, 1.5, 2.5, 3.5], "edges")])
# branch: [-20, 80, -180, -80]
# expected: [1500100, 1500200, 1500300, 1500400]
print(_grid_extents(g)) # branch: (-180, 800000, 180, 800400)
Curious how other functions handle this (because this didn't come up in any of the test files because it doesn't appear as if any are in meters), I had claude run through this on other functions and this may just be a broader problem we haven't been accounting for in function design. From claude:
My takeaway: |
|
@dylannelson Yeah, since this is a bigger issue than I thought, it should probably be a separate PR. I would guess this has more to do with units checking and conversion to the unit sphere for |
|
Closes #366 Closes #1243
Overview
dataarray.plot.contour()All plots are holoviews interactive plots by default. If matplotlib output is desired,
backend="matplotlib"gives a static image. Note that you can still get a static version of any holoviews plot by clicking thesavebutton in the top right of the plot. I think this arg is likely to have the highest potential for bugs, and we could even likely do without it, but I wanted to include because I'd prefer if we can improve it with time, and explore what it may be useful for (mostly flexibility in outputs). But if it's better to proceed without it let me know.The following new pages/edits to existing pages have been added to the docs:
plot.contour()documentation on readthedocs.plot.contour().plot.contour()plot.contour()Notebooks on readthedocsExpected Usage
Old plot.contour
Here's a kind of summary with before/after, as well as info on how this works on different data types, and scenarios where it doesn't work and why:
uxdsis aUxDataset;tempis aUxDataArraywith dims(time, lev, n_face)on the same grid..plot.contour()called on.plotmainuxdsUxDatasetAttributeErroruxds["var"].plot.contour()uxds.uxgridGridAttributeErroruxds["bottomDepth"]UxDataArray(n_face)ValueError: DataArray must be 2dhv.Overlay)UxDataArray(n_node)ValueError: DataArray must be 2dhv.Overlay)UxDataArray(n_edge)ValueError: DataArray must be 2dDataCenteringErrortemp.isel(time=[0], lev=[0])UxDataArray(time=1, lev=1, n_face)ValueError: DataArray must be 2dhv.Overlay)temp.isel(lev=0)UxDataArray(time=4, n_face)DimensionErrormain, whose plot was not a map. Select one time:.isel(time=0)tempUxDataArray(time, lev, n_face)ValueError: DataArray must be 2dDimensionErrortemp.mean("n_face"),temp.weighted_mean()UxDataArray(time, lev)QuadContourSet)temp.isel(time=0).mean("n_face")UxDataArray(lev)ValueError: DataArray must be 2dzonal_mean()/azimuthal_mean()resultsxr.DataArrayQuadContourSet)uxds["bottomDepth"].to_xarray()xr.DataArray(n_face)ValueError: DataArray must be 2dUxDataArrayinstead.Only one call that returned something on
mainraises with this PR: data with a grid dimension and a second dimension longer than 1, for example(time=4, n_face). Onmainxarray drew it with the face number as an axis, which is not a map. Everything else that plotted onmainplots the same way.Known Concerns, not planned for this PR:
PR Checklist
General
Testing & Benchmarking
Documentation and Examples
docs/api.rst; internal (private) function names start with an underscore (_)AI Disclosure
AI Usage:
Claude Code