usp-utils/post_proc: vertical cross-section plotting tool#16
Open
daniloceano wants to merge 9 commits into
Open
usp-utils/post_proc: vertical cross-section plotting tool#16daniloceano wants to merge 9 commits into
daniloceano wants to merge 9 commits into
Conversation
New post-processing tool for vertical cross-sections on the native MPAS grid, reusing the shared building blocks from mpas_viz.py (file opening, derived mesh coords, color-scale helpers, colorbar, argparse conventions). Two modes: - --levels-only: draw the zgrid interfaces + terrain along a transect (see the vertical level structure without needing a field). - -v <field>: color a 3D field (theta, qv, rho, ...) along the transect with height (m MSL) on the vertical axis and the terrain filled underneath. Transect defined by two points (--start/--end) or a constant lat/lon line (--lat/--lon). Cells are nearest-neighbour sampled (haversine), consecutive duplicates collapsed, columns placed at accumulated great-circle distance. Out-of-domain sample points are detected via the mesh resolution and skipped so an off-domain transect does not silently snap onto boundary cells. --by-index and --zmax control the vertical axis; with --zmax the color scale uses only the visible part so upper-level values do not wash out the plot. The vertical grid (zgrid) lives only in *.init.nc; a dedicated resolver points users there (static.nc has ter but not zgrid). README documents the tool. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
There was a problem hiding this comment.
Pull request overview
Adds a new post-processing plotting script to generate vertical cross-sections (transects) on the native MPAS grid (no regridding), and documents its usage alongside the existing mpas_viz.py tooling.
Changes:
- Introduces
mpas_cross_section.pyto plot either vertical level structure (--levels-only) or a colored 3D field along a transect. - Adds README documentation describing cross-section modes, transect definition options, and examples.
- Updates author attribution in the plotting scripts README to include the new tool.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 3 comments.
| File | Description |
|---|---|
usp-utils/post_proc/plotting_scripts/README.md |
Documents the new vertical cross-section tool and adds usage examples/options. |
usp-utils/post_proc/plotting_scripts/mpas_cross_section.py |
New CLI tool implementing transect sampling, grid/field validation, and plotting. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+430
to
+444
| da = ds[vname] | ||
| if 'nVertLevels' not in da.dims: | ||
| hint = "" | ||
| if vname == 'ter' or ('nCells' in da.dims and da.ndim <= 2): | ||
| hint = ("\n To see the terrain along the transect (it is " | ||
| "drawn as the filled bottom), use --levels-only.") | ||
| raise SystemExit( | ||
| f"\nERROR: '{vname}' has dims {da.dims}; a colored cross-section " | ||
| "needs a 3D field on (nCells, nVertLevels). Run without -v to " | ||
| f"list the colorable variables.{hint}") | ||
| if 'Time' in da.dims: | ||
| da = da.isel(Time=tindex) | ||
|
|
||
| field = da.values # (nCells, nVertLevels) | ||
| field2d = field[cells, :].T # (nVertLevels, ncols) |
Comment on lines
+240
to
+244
| for la, lo in zip(plat, plon): | ||
| d = _haversine_km(lat_cell, lon_cell, la, lo) | ||
| j = int(np.argmin(d)) | ||
| if max_snap_km is not None and d[j] > max_snap_km: | ||
| n_outside += 1 |
Contributor
There was a problem hiding this comment.
I suggest to be fixed
| ax.set_title(f"{title}\n{subtitle}", fontsize=11) | ||
| ax.grid(True, alpha=0.2, linewidth=0.4) | ||
|
|
||
| ds.close() |
…gh cmap mpas_cross_section.py: - Wind overlay decomposed relative to the transect orientation: in-plane arrows (along-transect + vertical, w scaled by --w-exag for a qualitative tilt) plus transect-normal symbols following the standard meteorological convention (filled dot ⊙ = towards viewer / out of page, cross ⊗ = away / into page, size ∝ magnitude). The along-transect tangent is smoothed to damp the spurious normal component from the wiggly nearest-cell path. New flags -u/-v_wind/-w, --w-exag, --wind-stride, --wind-lstride. Inputs are the cell-centered reconstructed winds (from diag/history); init.nc only carries edge-normal u. README.md: - Split into a copy-paste "Validation walkthrough" (fill-in shell variables + 5 staged checks: horizontal grid, terrain, vertical grid/initial conditions, surface forcing, model output) and per-script option references. Default colormap -> Spectral_r (red = higher values) in both mpas_cross_section.py and mpas_viz.py, for a consistent, intuitive scale across the tools. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…xtent mpas_cross_section.py: - Animate a transect: files expand into one timeline (as in mpas_viz), one step gives a still, several give an .mp4/.gif. The transect geometry is resolved once (fixed in time) and only the field/wind are re-read per frame; consistent color scale across frames. New flags --tstart/--tend, --list-times, --fps; -f accepts globs. run() refactored into build_geometry / render_field_frame. - Cartopy locator inset (top-left) drawing the transect as a red line on a coastline map, auto-framed to the mesh footprint via auto_extent(); --no-inset disables it. mpas_viz.py: - --auto-extent frames the map to the mesh footprint, or to the non-NaN region of the plotted field (e.g. sst over the ocean), instead of global. Explicit -lat_min/-lon_min/... still take priority. README.md: document animation, the locator inset and --auto-extent; make the walkthrough paths/MESH generic and reframe $BOUNDS as a zoom usable on global variable-resolution meshes too. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Files like sfc_update/diag carry the field but not latCell/lonCell, so the coordinates must come from -gf/--gridfile. compute_auto_extent now looks in the plotted file first, then in the grid file, instead of giving up (which silently fell back to a global extent, so --auto-extent appeared to do nothing). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…o plugin Root cause of the crash hit in real use: with neither the ffmpeg CLI nor the imageio-ffmpeg plugin installed, imageio.get_writer(outfile, fps=fps) for a .mp4 does not raise ImportError (the only case _stitch_pngs handled) — it silently resolves to an unrelated always-available plugin (tifffile), which then crashes deep inside with "TiffWriter.write() got an unexpected keyword argument 'fps'" instead of naming the real, missing dependency. Because the crash happened before the frame cleanup line ran, it also left ~300 temp PNGs (and an 8-byte corrupt .mp4) behind in the working tree. Fix in both mpas_viz.py and mpas_cross_section.py (they share _stitch_pngs): - Detect the video backend explicitly: ffmpeg CLI (shutil.which) first, then the imageio-ffmpeg plugin; if neither is present, exit with a clear message naming the install options (conda ffmpeg / pip imageio-ffmpeg / use .gif). - Wrap frame combination in try/finally so temp frames are always removed, even when video writing fails. - Render frames into a scratch tempfile.mkdtemp() directory instead of the current working directory, and rmtree it in a finally block, so a crash can no longer scatter debris across the repo. Also removes the leftover debris (329 orphaned frame PNGs + the corrupt t2m_anim.mp4) from the run that hit this bug, and documents the FFmpeg requirement more prominently in the README. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Real run hit this next, right after the previous fix let ffmpeg actually be invoked: matplotlib's bbox_inches='tight' crop rarely lands on an even pixel count (here 1375x841), and libx264 + yuv420p requires even width/height, so the encoder failed immediately with "width not divisible by 2" and produced no output. Add '-vf scale=trunc(iw/2)*2:trunc(ih/2)*2' to force even dimensions — the standard ffmpeg fix for this. Also drop check=True (which raised a full subprocess.CalledProcessError traceback) in favor of checking the return code and raising a short, clear SystemExit — the useful diagnostic is ffmpeg's own stderr, already streamed live to the terminal; the traceback added nothing. Cleaned up debris from the run that hit this: the empty t2m_anim.mp4(.filelist.txt) in runs/meqbr_05km/, and two orphaned /tmp/mpas_viz_frames_* directories left by hard-killed (timeout) test runs during development, not by user code. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
mpas_viz.py and mpas_cross_section.py now accept either an integer timeline index or a datetime string (e.g. '2021-11-01_00:00:00' or date-prefix '2021-11-01') for --tstart, --tend and -t. The string is matched against the xtime labels in the timeline; the first entry whose xtime starts with the supplied value is used. Integer indices continue to work unchanged. README updated with examples and option descriptions.
…hrough - All check_0N output filenames now carry $MESH as prefix so it is immediately clear which mesh/file produced each figure. - Step 5 (model output): added t2m animation with 10-m wind vectors (--stride) and converted cross-section step to an animation with wind overlay (--wind-stride / --wind-lstride). - Added inline comments with suggested --stride / --wind-stride values for three typical horizontal spacings (~50, ~15, ~5 km).
add_wind_vectors now receives pre-resolved lat/lon arrays instead of the dataset, fixing a KeyError when the data file (e.g. diag.*.nc) does not carry latitude/longitude variables. _resolve_cell_coords() is called in render_one_frame to look up latCell/lonCell (or latitude/longitude) from the data file first, then from the gridfile (-gf). A clear warning is printed and vectors are skipped (instead of crashing) when no coordinates are found.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds
mpas_cross_section.py, a post-processing tool that draws vertical cross-sections (transects) on the native MPAS grid — no regridding. It reuses the shared building blocks frommpas_viz.py(file opening, derived mesh coords, color-scale helpers, colorbar, argparse conventions);mpas_viz.pyis not modified.What it does
Two modes:
--levels-only— draw thezgridinterfaces + terrain along a transect (see the vertical level structure without needing a field).-v <field>— color a 3D field (theta,qv,rho, …) along the transect, height (m MSL) on the vertical axis, terrain filled underneath.Transect defined by two points (
--start/--end) or a constant lat/lon line (--lat/--lon). Cells are nearest-neighbour sampled (haversine), consecutive duplicates collapsed, columns placed at accumulated great-circle distance.--by-indexand--zmaxcontrol the vertical axis.Robustness
--zmax, the scale uses only values below the cap so upper-level values (e.g. stratospherictheta) don’t wash out the plot.zgridlives only in*.init.nc; a dedicated resolver points users there (*.static.nchasterbut notzgrid). Colored plots reject 2D fields and suggest--levels-only.Testing
Exercised on
runs/meqbr_05km/meqbr_05km.init.nc(nCells=76813, 55 levels): variable listing,--levels-only, coloredtheta/qv,--by-index,--zmax, and the error paths (no transect, conflicting--lat/--lon, wrong grid file, 2D field). Figures inspected visually and physically consistent.🤖 Generated with Claude Code