Repository navigation
Added .plot.contour() function and added docs for visualizations - #1812
dylannelson wants to merge 25 commits into
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. |
|
@cmdupuis3 I have added an option for the labels arg to be a ratio, from 0-1, as well as True (1) or False (0), which covers # 3, I think the other 2 either you have worked on or would heavily change the functions you have work on. Over chat we discussed how some of these changes can/should come after this is merged, so we can come back to that later on. |
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
| corners = np.zeros(len(n_nodes_per_face), dtype=np.int64) | ||
| for face in range(len(n_nodes_per_face)): | ||
| n = n_nodes_per_face[face] | ||
| if n < 4: | ||
| continue | ||
| x = np.empty(n) | ||
| y = np.empty(n) |
There was a problem hiding this comment.
This will likely permit some nice performance optimizations later due to the tiny numpy arrays, as per #1648. It's probably okay to push the performance optimization to a later PR, though. Feel free to mark this as resolved once you read it, I just wanted to leave this note/link for future reference when looking at 1648.
There was a problem hiding this comment.
+1, the tiny numpy arrays should be fixed.
There was a problem hiding this comment.
Would you be able to make it possible with this PR if you sat together, @dylannelson and @Sevans711 ?
My thinking is that since we solved it once, let's do it here as well rather than introducing the issue again and solving it later.
There was a problem hiding this comment.
I'm working on a fix now!
There was a problem hiding this comment.
"Fixed" on sevans/contour-maps-avoid-tiny-numpy-in-numba branch, by lifting the x and y buffer initializations to be outside the loop. Putting "fixed" in quotes because I didn't actually see any noticeable speedup when calling arr.plot.contour(). This is probably because the _fan_corners accounts for at most a few percent of the .plot.contour() calls, in all cases I tested. So, this is probably some premature optimization and not actually necessary… though it also doesn't hurt. I'd personally be fine regardless of whether this gets merged into this PR or not.
I also separately checked what happens if I keep the buffers in the loop and use prange instead of range, but that also didn't show any noticeable speedups (again, this method is only a few percent of the .plot.contour() calls).
More details: timed across tutorial examples "outCSne30-vortex" and "mpas-QU-oQU480" (ocean grid, might have nontrivial _fan_corners), as well as healpix grids (I tried example UxDataArray like: grid = ux.Grid.from_healpix(zoom=zoom); areas = grid1.compute_face_areas(as_uxarray=True); arr = (areas1* np.arange(len(areas))), with zoom values of 6, 7, 8, and 9, which reaches up to 3M faces). Also tried using both method="edges" and method="interpolated".
Even with "interpolated", .plot.contour() is taking <30 seconds on 3M faces, which to me feels "fast enough" for initially shipping the feature?
| """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?)
| "Contours are only supported for face-centered or node-centered data, " | ||
| f"but got dimension '{dim}'" | ||
| ) | ||
| return np.asarray(uxda.squeeze().values, dtype=float), dim |
There was a problem hiding this comment.
I pushed back on this squeeze() in my main review comment; leaving an inline comment to help track this topic.
There was a problem hiding this comment.
This one I'm really confused about. The original concern was "presence of "time" ... may be mishandling coords from source to the output", and when time is 1D, it can cause functions to crash. From my understanding, squeezing solves the error in this case and is the recommended solution. But the proposed change is to fix it sometimes and not other's because we want to "follow in the style of xarray when possible"? If we wanted to make user's workflows/behavior more consistent, why would having asymmetrical behaviors like that be beneficial? Feels like a lose-lose situation either way, either be different than xarray, or have similar functions error differently to the same coordinates. I feel like it would be easier to just squeeze to clean the dims and provide less resistance for the user, vs in the proposed case, which is to assume the user knows the asymmetrical relationships of plotting functions and we avoid friction for that user. Since this isn't a discrepancy with xarray that anyone has actually raised a concern about, and it doesn't seem to make the "extra dims" problem any easier for a user, could we discuss more about it with others before letting it be a blocker here?
There was a problem hiding this comment.
If we don't want to delay for a larger discussion, we should ship the feature with whatever is easiest to revert / change later without breaking any existing workflows.
If we ship with "automatically squeezes", reverting could break existing code.
If we ship with "does not automatically squeeze", then later decide to update to "automatically squeezes", no existing code would break, it would just add support for more use-cases. This one is easier to revert / change later.
I'm not sure if I understand the original concern in this case. Maybe it would be easiest to briefly talk through this out loud and report back here afterwards? Leaving a note about this in upcoming internal meeting agenda.
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
5 participants
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