Skip to content

General Improvements + bug fixes + multi scale meshing support for global meshes#87

Merged
CHLNDDEV merged 23 commits intomasterfrom
improvements
Dec 1, 2025
Merged

General Improvements + bug fixes + multi scale meshing support for global meshes#87
CHLNDDEV merged 23 commits intomasterfrom
improvements

Conversation

@krober10nd
Copy link
Copy Markdown
Collaborator

@krober10nd krober10nd commented Nov 11, 2025

Overview

This PR refactors and hardens the meshing pipeline with multiscale (global + regional) support, CRS-aware grid handling, and consistent bathymetric gradient sizing.

Key Changes

SDF

  • Ignores non-finite query points; assigns large positive distance.
  • Removed deprecated n_jobs param in cKDTree.
  • Propagates CRS / stereo metadata through domain ops.

Multiscale Meshing

  • Centralized validation (ordering, bbox, stereo, CRS).
  • Auto-detects mixed CRS grids and reprojects on the fly.
  • Honors per-nest stereo; unified blending.
  • Robust hmin handling from valid finite grid values.

Sizing Functions & Grids

  • Added Grid.reproject_to() for reliable CRS alignment.
  • compute_minimum reprojects before elementwise min.
  • Refactored bathymetric_gradient_sizing_function:
    • Computes gradients in meters with latitude-aware scaling.
    • Optional integer coarsen controls DEM downsampling.
    • Removes ad-hoc dx/2 scaling; preserves resolution.
    • Enforces finite size values; replaces invalids with median.

DEM & Orientation

  • Unified raster orientation between full and windowed reads.

Bug Fixes

  • Fixed KDTree crashes from NaN coordinates.
  • Corrected reprojection-related flips/rotations.
  • Prevented invalid or zero min_edge_length.
  • Replaced fragile CRS equality checks with explicit reprojection.

Performance / Stability

  • Graceful hmin recovery in multiscale blending.
  • Localized reprojection minimizes overhead.
  • Optional coarsening for DEM sizing efficiency.

Compatibility & Validation

  • APIs unchanged; optional params only.
  • Verified stable hmin and correct CRS alignment across single- and multiscale tests.
  • Vendored inpoly-python to avoid complications with installation on some systems.

Summary

Improves safety, CRS consistency, and multiscale robustness while eliminating silent resolution loss and preparing for advanced sizing diagnostics.

…nd multiscale_signed_distance_function in signed_distance_function.py; metadata propagation preserved.

Eliminated all set() usage on CRS objects (potentially unhashable) by dropping those validations entirely per your direction.
Centralized all CRS / stereo validation logic in _validate_multiscale_domains inside mesh_generator.py.
Added implicit global detection (EPSG:4326 + global-like bbox) and error if mixed CRS without stereo=True or not first.
Refined missing-CRS handling: now only enforced when some CRS info is present or implicit global detected; pure bbox-only workflows remain allowed.
Removed duplicate “missing CRS” message branch.
Replaced string CRS comparisons with robust pyproj.CRS.equals checks; retained get_crs_string only for error message formatting.
Edge length vs domain CRS comparison now uses CRS.equals with guarded exception handling.
…t sizing

* signed_distance_function:
    * Guard non-finite query points (NaN/Inf) in SDF and covering; assign large positive distance to invalid rows
    * Use cKDTree.query(x, k=1) without deprecated/unsupported n_jobs/workers args
* mesh_generator:
    * Stabilize _compute_forces by replacing non-finite or <=0 hbars with median valid values (prevents NaNs in updates)
* grid:
    * compute_minimum: drop same-CRS assertion; reproject each grid onto the base lattice using pyproj (flatten→transform→reshape) to avoid axis/broadcasting issues; compute min on aligned arrays; derive hmin from positive finite values only
    * Add Grid.reproject_to(target_grid) helper for CRS/lattice-aligned resampling
* edgefx:
    * bathymetric_gradient_sizing_function rewritten to:
        * Build sizing on the DEM lattice (same nx/ny/dx/dy/crs); no ad-hoc dx halving/doubling
        * Compute gradients in physical units; convert meters↔degrees for geographic CRS using mean-latitude scales
        * Enforce min/max bounds robustly; add coarsen (int≥1) for optional downsampling
        * Accept crs kwarg for backward compatibility (ignored; DEM’s CRS is authoritative)
* geodata:
    * Normalize DEM full-read orientation (transpose + flip) so x/y/array orientation is consistent with windowed reads
Effects
* Fixes “x must be finite” KDTree error during SDF evaluation
* Correctly aligns hfun_coarse and hfun_bathy across mixed CRS without inversion
* Produces stable, positive min_edge_length for combined sizing
* No breaking API changes; adds helper and extends function signature (backwards compatible)
Copilot AI review requested due to automatic review settings November 11, 2025 23:59
Copy link
Copy Markdown

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull Request Overview

This PR introduces comprehensive support for multiscale (global + regional) meshing with CRS-aware grid handling, improved bathymetric gradient sizing, and enhanced robustness through validation and error handling. The changes centralize domain validation, add reprojection capabilities, and fix several bugs while maintaining API compatibility.

Key Changes

  • Multiscale domain validation: Added _validate_multiscale_domains() to enforce ordering, CRS compatibility, and bbox containment rules for global+regional mesh workflows
  • CRS-aware grid operations: Implemented Grid.reproject_to() and enhanced compute_minimum() to handle mixed-CRS grids through coordinate transformation
  • Improved bathymetric gradient sizing: Refactored gradient computation to use proper meter-based scaling with latitude awareness, added optional coarsen parameter for DEM downsampling, and replaced the ignored crs parameter with dem.crs

Reviewed Changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated 7 comments.

Show a summary per file
File Description
oceanmesh/signed_distance_function.py Added NaN/non-finite point handling in SDF evaluation; propagated CRS/stereo metadata through Domain, Union, Intersection, and Difference classes
oceanmesh/region.py Added validation helpers for CRS compatibility, bbox containment, and global bbox detection to support multiscale meshing
oceanmesh/mesh_generator.py Implemented comprehensive domain validation with _validate_multiscale_domains(); added hmin sanitization logic; enhanced force computation with non-finite value guards
oceanmesh/grid.py Fixed bbox typo in interpolate_to(); refactored compute_minimum() to handle CRS reprojection; added reproject_to() method for cross-CRS grid alignment
oceanmesh/geodata.py Unified raster orientation between full and windowed reads; added bbox transformation logic; improved handling of un-georeferenced rasters; retained stereo flag in Shoreline
oceanmesh/edgefx.py Refactored bathymetric gradient sizing with proper latitude-aware meter conversion; added coarsen parameter for DEM downsampling; changed crs parameter default to None (now uses dem.crs)

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

type_of_filter="lowpass",
filter_cutoffs=1000,
crs="EPSG:4326",
coarsen=1,
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The coarsen parameter default value is 1 (an integer), but the validation at line 483 checks int(coarsen) != coarsen, which will fail for float inputs like 1.0. For better user experience, consider using not isinstance(coarsen, int) or coarsen < 1 instead, or explicitly convert with coarsen = int(coarsen) if float inputs should be accepted and truncated.

Copilot uses AI. Check for mistakes.
p=3,
nnear=28,
blend_width=1000,
domain_metadata=None,
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The domain_metadata parameter is added to the function signature but is never used in the function body. If this parameter is reserved for future use, consider adding a comment explaining its purpose. Otherwise, remove it to avoid confusion.

Suggested change
domain_metadata=None,

Copilot uses AI. Check for mistakes.
triang = tri.Triangulation(points[:, 0], points[:, 1], cells)
fig, ax = plt.subplots(figsize=(10, 10))
ax.triplot(triang)
ax.triplot(triang, lw=0.1)
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The addition of lw=0.1 is a minor visualization improvement that may be subjective. While it can improve plotting for dense meshes, consider making this configurable via a keyword argument to allow users to customize based on their mesh density.

Copilot uses AI. Check for mistakes.
Comment on lines +380 to +382
ok_global, msg_global = validate_crs_compatible(gcrs, gcrs)
if not ok_global:
errors.append(msg_global)
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The function validate_crs_compatible is called with the same CRS twice (gcrs, gcrs), but it expects two different arguments: global_crs and regional_crs. This validation is meant to check if the global CRS is EPSG:4326. Since both arguments are identical, the function will always validate the CRS against itself. Consider using a dedicated validation function or pass a known EPSG:4326 CRS as the second argument for proper validation.

Suggested change
ok_global, msg_global = validate_crs_compatible(gcrs, gcrs)
if not ok_global:
errors.append(msg_global)
try:
if CRS.from_user_input(gcrs).to_epsg() != 4326:
errors.append(f"Global domain CRS must be EPSG:4326, got {gcrs_str}.")
except Exception:
errors.append(f"Global domain CRS could not be parsed: {gcrs_str}.")

Copilot uses AI. Check for mistakes.
bbox = (xmin, xmax, ymin, ymax)

# Warn if user requested output CRS different from raster CRS (no reprojection performed here)
if (src_crs is not None) and (desired_crs is not None) and (src_crs != desired_crs) and not ((src_crs is None) or src.transform == Affine.identity):
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The boolean condition has a redundant check. After verifying (src_crs is not None) earlier in the expression, the final clause not ((src_crs is None) or ...) contains a redundant src_crs is None check that can never be true. Simplify to: if (src_crs is not None) and (desired_crs is not None) and (src_crs != desired_crs) and (src.transform != Affine.identity):

Suggested change
if (src_crs is not None) and (desired_crs is not None) and (src_crs != desired_crs) and not ((src_crs is None) or src.transform == Affine.identity):
if (src_crs is not None) and (desired_crs is not None) and (src_crs != desired_crs) and (src.transform != Affine.identity):

Copilot uses AI. Check for mistakes.
Comment on lines +376 to +378
gcrs_str = get_crs_string(gcrs)
except Exception:
gcrs_str = str(gcrs)
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable gcrs_str is not used.

Suggested change
gcrs_str = get_crs_string(gcrs)
except Exception:
gcrs_str = str(gcrs)
get_crs_string(gcrs)
except Exception:
str(gcrs)

Copilot uses AI. Check for mistakes.
Comment on lines +376 to +378
gcrs_str = get_crs_string(gcrs)
except Exception:
gcrs_str = str(gcrs)
Copy link

Copilot AI Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable gcrs_str is not used.

Suggested change
gcrs_str = get_crs_string(gcrs)
except Exception:
gcrs_str = str(gcrs)
get_crs_string(gcrs)
except Exception:
str(gcrs)

Copilot uses AI. Check for mistakes.
…ation; add test

Summary of changes:

Docs (README):
Inserted new section “Global mesh generation with regional refinement” with full workflow example.
Added TOC entry for the new section.
Clarified unit usage (degrees for EPSG:4326 sizing grids) and realistic blend_width guidance (now 1.0e6 meters with explanatory comment).
Added note on automatic CRS/projection handling and support for projected regional CRSs.
Docs (code):
Expanded generate_multiscale_mesh docstring with workflow, automatic coordinate handling, example reference, stereo flag behavior, and notes.
Validation:
Made global domain CRS check explicitly require EPSG:4326 (removed self-comparison).
Strengthened multiscale domain validation messaging (CRS, bbox containment, stereo usage).
Mesh generation robustness:
Sanitized sizing grids’ hmin fallback and non-finite handling.
Hardened _compute_forces against invalid sizing values.
Added automatic recomputation of min_edge_length from positive finite grid values.
Sizing & grid utilities:
Refactored bathymetric_gradient_sizing_function to preserve DEM resolution, operate in meters, support coarsening, enforce bounds.
Added CRS-aware grid combination via compute_minimum reprojection alignment and Grid.reproject_to helper.
SDF / domain:
Added non-finite input sanitization for KDTree queries and removed deprecated cKDTree parameters.
Propagated stereo and CRS metadata through multiscale signed distance functions.
Multiscale meshing:
Simplified stereo handling—global domain first with stereo=True; regional domains in lat/lon (or projected) transparently transformed.
Added multiscale sizing blending improvements and stereo flag propagation to final union.
Tests:
Added test_global_regional_multiscale.py for end-to-end global + Australia regional refinement, with assertions on regional/transition presence and optional visualization output.
Adjusted iteration counts and thresholds for performance and CI stability.
Misc:
Updated README examples and inline comments to match new unit conventions and projection workflow.
…containment; runnable e2e test; plotting/logging fixes

README

Add “Global mesh generation with regional refinement” section with full workflow example
Insert TOC entry for the new section
Clarify units: sizing grids follow CRS (degrees for EPSG:4326); note internal conversions
Update blend_width guidance (use ~1e5–1e6 m; example uses 1.0e6) and comments
Add note that regional domains may use projected CRSs (e.g., UTM); blending across CRSs is automatic
mesh_generator.py

Docstring: expand generate_multiscale_mesh description with automatic coordinate handling, stereo propagation, example reference, and notes
Validation: explicitly require global domain CRS = EPSG:4326 (remove self-compare)
Stereo containment: when global domain is stereo, transform regional bbox to stereo before bbox containment check
Robustness: retain positive hmin fallback on sizing grids; propagate stereo=True to final blend if union is stereo
Improve error messages for multiscale validation
test_global_regional_multiscale.py

End-to-end global (stereo) + Australia regional refinement test
Convert mesh points to lat/lon for regional/transition classification
Add assertions: non-empty mesh, quality thresholds, regional presence, hmin ordering, regional areas smaller than global
Fix plotting to show global mesh in lat/lon with proper extents; optional PNG save under tests/output
Enable verbose logging for easier diagnosis; add main guard for direct execution
Runtime-tuned params (blend_width ~200 km, reduced iterations)
__init__.py

Export simp_qual at package level for use in tests/examples
Overall: Documents the mixed global+regional workflow, enforces correct CRS/stereo usage, fixes bbox containment across projections, and adds a runnable test and plotting/logging improvements.
… consistency and new tests

Shoreline now prefers a Region object (bbox + CRS) while preserving tuple backward compatibility; warns on explicit CRS overrides.
Added _infer_crs_from_coordinates heuristic to auto-switch from default EPSG:4326 when a projected-looking tuple bbox is passed.
Optimized initialization to avoid unnecessary shapefile reads for geographic bboxes.
Improved _read error message with active/native CRS, bbox, and remediation guidance.
Refactored DEM.__init__ to use Region-derived CRS cleanly and avoid tuple attribute misuse.
Updated all README examples (including global/multiscale) to use Region-first pattern with notes on projected CRS.
Added tests in test_geodata.py for Region precedence, CRS inference (projected vs geographic), logging, and enhanced error guidance.
Maintains backward compatibility; no breaking changes
Copy link
Copy Markdown

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull Request Overview

Copilot reviewed 23 out of 37 changed files in this pull request and generated 8 comments.

Comments suppressed due to low confidence (1)

tests/test_global_regional_multiscale.py:213

  • Variable global_only_quality is not used.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

filter_cutoffs=1000,
crs="EPSG:4326",
coarsen=1,
crs=None,
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changing the crs parameter default from 'EPSG:4326' to None is a breaking change that may affect existing code expecting the old default behavior. Consider deprecation warnings or documentation to help users migrate.

Suggested change
crs=None,
crs="EPSG:4326",

Copilot uses AI. Check for mistakes.
Comment on lines +657 to +658
_kwargs = dict(kwargs)
_kwargs.pop("max_iter", None)
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The code removes max_iter from kwargs to avoid duplicate parameter errors, but this pattern is fragile. Consider using explicit parameter passing or restructuring to avoid this defensive programming pattern.

Copilot uses AI. Check for mistakes.
regional_quality = quality[in_regional]
# global_only_quality intentionally unused; kept for potential comparative diagnostics
# noqa below suppresses flake8 unused variable warning.
global_only_quality = quality[~in_regional] # noqa: F841
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] Variable global_only_quality is computed but intentionally unused with a noqa comment. If this is for potential future diagnostics as the comment suggests, consider removing it until actually needed to keep the code cleaner.

Copilot uses AI. Check for mistakes.
finite_mask = np.isfinite(x).all(axis=1)

# Default: very large positive distance (outside domain)
d = np.full(n, 1.0e12, dtype=float)
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] Magic number 1.0e12 used for large positive distance. Consider defining this as a named constant (e.g., LARGE_POSITIVE_DISTANCE = 1.0e12) for better maintainability and documentation of intent.

Copilot uses AI. Check for mistakes.
setup.py Outdated
Comment on lines +21 to +25
# Build system: This package uses pybind11's setuptools integration to compile
# C++ extensions (HamiltonJacobi, delaunay_class, fast_geometry). CMake is NOT
# invoked for building oceanmesh itself; on Windows it may be required only by
# external tools (e.g., vcpkg) to build CGAL dependencies. Optional inpoly
# acceleration is enabled when the [fast] extra pulls in Cython.
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Excellent inline documentation explaining the build system architecture, CMake usage scope, and optional acceleration features. This helps maintainers understand the complex build setup.

Copilot uses AI. Check for mistakes.
Comment on lines +909 to +911
def _transform_bbox_to_src(_bbox_vals, _src_crs, _dst_crs):
if _src_crs is None or _dst_crs is None or _src_crs == _dst_crs:
return _bbox_vals
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] Nested function _transform_bbox_to_src defined inside __init__ adds complexity. Consider extracting to a module-level private function for better testability and readability.

Copilot uses AI. Check for mistakes.
assert int(len(bbox) / 2), "`dim` must be 2"


def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The noqa C901 comment indicates this function exceeds complexity limits. Consider refactoring into smaller validation functions (e.g., _validate_crs_compatibility, _validate_stereo_flags, _validate_bbox_containment) to improve maintainability and testability.

Copilot uses AI. Check for mistakes.
Comment on lines +260 to +263
print(
f"[inpoly2] Compiled kernel import failed: {_ex}. "
"Falling back to pure-Python implementation."
)
Copy link

Copilot AI Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Print statement may execute during import.

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FYI oceanmesh's GPL3 License infringes this License's terms

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, I plan on writing a quick inpoly predicate and replacing it avoiding the license issue.

# Regional sizing functions have been wrapped or transformed earlier so fh(p2)
# evaluates correctly on lat/lon even though points are maintained in stereo space.
p1 = p[bars].sum(1) / 2
x, y = to_lat_lon(p1[:, 0], p1[:, 1])
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since that you now use CRS and transformations, I would try to drop the to_lat_lon() and to_stereo() functions which are hardcoded.

@tomsail
Copy link
Copy Markdown
Contributor

tomsail commented Nov 24, 2025

Hi Keith, I just had a quick look.
I think the stereographic handling needs to be revised (see the quality of the boundary triangles in polar locations).
This mainly my fault really, because I have hardcoded the transformations from stereo to lat/lon.

Having had experience on seamsh, I think the best - to simplify things - would be to convert global coastlines systematically to use the ccrs.NorthPolarStereo() projection.

This would simplify the scaling factor:
image
source

Has done in seamsh stereographic example

By doing so, you would:

  1. get rid of the hardcoded function for lat/lon <> stereo conversions
  2. have a unique scaling factor that applies well to the same projection. thus avoiding scaling errors and ensuring consistency.
  3. be able to provide different CRS for the global coastlines

@CHLNDDEV
Copy link
Copy Markdown
Owner

Hi Keith, I just had a quick look. I think the stereographic handling needs to be revised (see the quality of the boundary triangles in polar locations). This mainly my fault really, because I have hardcoded the transformations from stereo to lat/lon.

Having had experience on seamsh, I think the best - to simplify things - would be to convert global coastlines systematically to use the ccrs.NorthPolarStereo() projection.

This would simplify the scaling factor: image source

Has done in seamsh stereographic example

By doing so, you would:

  1. get rid of the hardcoded function for lat/lon <> stereo conversions
  2. have a unique scaling factor that applies well to the same projection. thus avoiding scaling errors and ensuring consistency.
  3. be able to provide different CRS for the global coastlines

Okay sounds good, do you want to add that to this PR?

@CHLNDDEV
Copy link
Copy Markdown
Owner

CHLNDDEV commented Dec 1, 2025

@tomsail I'll tackle the improvements to global scale meshing on the subsequent PR to this one. I have some ideas.

@CHLNDDEV CHLNDDEV merged commit b9bc6a9 into master Dec 1, 2025
5 checks passed
@CHLNDDEV CHLNDDEV deleted the improvements branch December 1, 2025 14:40
@tomsail
Copy link
Copy Markdown
Contributor

tomsail commented Dec 1, 2025

No worries. Thanks Keith and kudos on this PR
Always happy to help and review. Been quite busy to do anything more..

@krober10nd
Copy link
Copy Markdown
Collaborator Author

No worries. Thanks Keith and kudos on this PR Always happy to help and review. Been quite busy to do anything more..

Thanks Thomas! Totally understand. I greatly appreciate the comments. I'll ask for some more if that's alright on the next one that I'll start probably this evening or tomorrow evening.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants