Repository navigation
Added .plot.contour() function and added docs for visualizations - #1812
dylannelson wants to merge 15 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:
|
|
@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. |
|
(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
| The contour lines are computed on the unstructured grid itself, without | ||
| regridding the data, and are returned as HoloViews elements that can be | ||
| overlaid on other plots with ``*``. |
There was a problem hiding this comment.
I don't think it needs to be individual plotting methods' docstrings' job to mention how holoviews elements can be combined, at least not this high up in the docstring. It is also already discoverable by the examples in the docstring! Also, more importantly, using * is not the only want to combine holoviews elements; ror example, you can use + to create a layout with multiple plots.
I'd recommend to either remove the mention of "can be overlaid on other plots with *", or, if you prefer, add a mention about the ability to use + as well.
There was a problem hiding this comment.
I don't really agree with this entirely, I think that you will likely want to always put contours on top of other data (*), and not really ever add it to a separate plot on it's own (+) so I don't know why the added (+) information would be needed. I think it's of the few plot types that is almost exclusively used on top of other data, hence the distinction. But if I'm in the minority I can delete it
There was a problem hiding this comment.
You do see pure contour plots occasionally, it's the standard with pressure, streamfunction, geopotential height, etc.
There was a problem hiding this comment.
I still imagine many would still at least have at least coastlines, which require "*", but I get it, it's gone now
| projection : ccrs.Projection, optional | ||
| The map projection to use. Pass the same projection as the plot the | ||
| contours are overlaid on. |
There was a problem hiding this comment.
Minor nit, the contours might not be overlaid on a plot at all!
| def _fan_corners( | ||
| face_node_connectivity: np.ndarray, | ||
| n_nodes_per_face: np.ndarray, | ||
| lon: np.ndarray, | ||
| lat: np.ndarray, | ||
| ) -> np.ndarray: | ||
| """For each face, the position of a corner that can see every other corner, so that | ||
| a fan of triangles from it covers the face exactly. 0 is tried first; -1 means that | ||
| no corner can. Faces that touch or contain a pole, or have no area, keep 0: a pole is | ||
| one point drawn at an arbitrary longitude, so their shape in longitude and latitude | ||
| says nothing about dents.""" |
There was a problem hiding this comment.
I think you don't need to add support for this? Uxarray only supports convex grids right now, see #1615. If this could be relevant on a convex grid, please clarify. If this is only to handle non-convex grids, please clarify that in the docstring instead. Since it already exists and is an internal method, I won't push for removing the function right now, since it could easily be removed or edited later if needed; just updating comment/docstring would be sufficient for clarity here.
| 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.
| "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?
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