Skip to content

Commit 6f7e3b3

Browse files
asinghvi17simone-silvestrijuliasloan25claude
authored
General improvements (#60)
Co-authored-by: Simone Silvestri <silvestri.simone0@gmail.com> Co-authored-by: Anshul Singhvi <anshulsinghvi@gmail.com> Co-authored-by: Julia Sloan <51397186+juliasloan25@users.noreply.github.com> Co-authored-by: Claude Opus 4.6 <noreply@anthropic.com>
1 parent 54ea967 commit 6f7e3b3

20 files changed

+677
-183
lines changed

CLAUDE.md

Lines changed: 57 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
# CLAUDE.md
22

3-
Guidance for Claude Code when working with this repository.
3+
This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
44

55
## Project Overview
66

7-
ConservativeRegridding.jl performs area-weighted conservative regridding between polygon grids. "Conservative" means the area-weighted mean is preserved. The package computes intersection areas between source/destination grid cell pairs to create averaging weights.
7+
ConservativeRegridding.jl performs area-weighted conservative regridding between polygon grids on planes or spheres. "Conservative" means the area-weighted mean is preserved during regridding. The package computes intersection areas between source/destination grid cell pairs, stores them as a sparse matrix, and uses that matrix for fast forward and backward regridding.
88

99
## Common Commands
1010

11-
**Important**: Always use `julia --project=docs` when running scripts to access all dependencies.
11+
**Important**: Always use `julia --project=docs` when running scripts and investigating — the docs environment has diagnostic packages not available in test.
1212

1313
```bash
1414
# Run all tests
@@ -22,89 +22,92 @@ julia --project=test -e 'include("test/usecases/simple.jl")'
2222
julia --project=test -e 'include("test/trees/grids.jl")'
2323
julia --project=test -e 'include("test/trees/quadtree_cursors.jl")'
2424

25-
# Start Julia REPL with the package loaded
25+
# Run examples/scripts or load the package interactively
26+
julia --project=docs examples/speedy_to_speedy.jl
2627
julia --project=docs -e 'using ConservativeRegridding'
2728
```
2829

29-
The project uses separate Project.toml files in `test/`, `docs/`, and `examples/`.
30+
The repo uses Julia's workspace feature (`[workspace]` in Project.toml) with separate environments in `test/`, `docs/`, and `examples/`.
3031

3132
## Architecture
3233

33-
### Core Types and Functions
34+
### Three-Layer Design
3435

35-
**`Regridder`** (`src/regridder/regridder.jl`): Main type storing intersection areas as a sparse matrix, grid cell areas, and work arrays.
36+
```
37+
Regridder (sparse matrix + area vectors)
38+
↓ constructed by
39+
Intersection Areas (dual DFS candidate search + parallel intersection computation)
40+
↓ operates on
41+
Trees Module (grid representations + quadtree cursors for spatial indexing)
42+
```
3643

37-
- `transpose(regridder)` returns a regridder for the reverse direction (shares underlying data)
38-
- `normalize!` scales the intersection matrix by its maximum value
44+
### Regridder (`src/regridder/`)
3945

40-
**`Regridder(dst, src; normalize=true, intersection_operator=..., threaded=True())`**: Constructor that:
41-
1. Determines manifold (Planar or Spherical) from input grids
42-
2. Converts grids to spatial trees via `Trees.treeify`
43-
3. Performs multithreaded dual depth-first search to find intersecting polygon pairs
44-
4. Computes intersection areas via `DefaultIntersectionOperator` (user-overridable)
46+
**`Regridder`** (`regridder.jl`): Stores a sparse intersection matrix, source/destination area vectors, and temporary work arrays.
47+
- Constructor: `Regridder(dst, src)` auto-detects manifold, treeifies grids, runs dual DFS, computes intersections
48+
- `transpose(regridder)` returns reverse-direction regridder sharing underlying data (no copy)
49+
- `normalize!` scales intersection matrix by its maximum value
4550

46-
**`regrid!(dst_field, regridder, src_field)`** (`src/regridder/regrid.jl`): Performs regridding: `dst = (A * src) / dst_areas`
51+
**`regrid!`** (`regrid.jl`): `dst = (A * src) / dst_areas`. Handles dense and non-contiguous arrays, copies to temporary arrays when needed for BLAS performance.
4752

48-
**`DefaultIntersectionOperator(manifold)`**: Dispatches to the appropriate intersection algorithm:
53+
**`intersection_areas`** (`intersection_areas.jl`): Two-phase approach:
54+
1. Dual DFS through spatial trees to find candidate cell pairs (extent-based pruning)
55+
2. Parallel intersection area computation on candidates, partitioned into `nthreads * 4` chunks via ChunkSplitters
56+
57+
**Intersection operators** dispatch on manifold:
4958
- Planar: `FosterHormannClipping`
5059
- Spherical: `ConvexConvexSutherlandHodgman`
5160

52-
### Trees Submodule (`src/trees/`)
61+
### Trees Module (`src/trees/Trees.jl`)
5362

54-
Provides quadtree-based spatial indexing for matrix-shaped grids.
63+
**Grid types** (`grids.jl`), all `<: AbstractCurvilinearGrid`, parameterized by manifold `M`:
64+
- `ExplicitPolygonGrid{M}`: Wraps a matrix of pre-computed polygons
65+
- `CellBasedGrid{M}`: Builds polygons on-the-fly from (n+1)×(m+1) corner point matrix
66+
- `RegularGrid{M}`: Regular lon/lat grids from 1D coordinate vectors
5567

56-
**`treeify(manifold, grid)`** (`src/trees/interfaces.jl`): Converts grid representations into `SpatialTreeInterface`-compliant trees:
57-
- Matrices of polygons -> `ExplicitPolygonGrid` + `TopDownQuadtreeCursor`
58-
- Matrices of points -> `CellBasedGrid` + `TopDownQuadtreeCursor`
59-
- Iterables of polygons -> `FlatNoTree`
60-
- Existing spatial trees -> pass-through
68+
**Required interface** for `AbstractCurvilinearGrid`: `getcell(grid, i, j)`, `ncells(grid, dim)`, `cell_range_extent(grid, irange, jrange)`
6169

62-
**Grid types** (`src/trees/grids.jl`), all parameterized by manifold `M`:
63-
- `ExplicitPolygonGrid{M}`: Wraps a matrix of pre-computed polygons
64-
- `CellBasedGrid{M}`: Builds polygons on-the-fly from corner points
65-
- `RegularGrid{M}`: For regular lon/lat grids defined by 1D vectors
70+
**Quadtree cursors** (`quadtree_cursors.jl`): Implement `SpatialTreeInterface` on top of grids:
71+
- `TopDownQuadtreeCursor`: Recursive subdivision by index ranges (primary cursor used)
72+
- `QuadtreeCursor`: Bottom-up traversal with level-based indexing
6673

67-
**`AbstractCurvilinearGrid`** interface (`src/trees/interfaces.jl`):
68-
- `getcell(grid, i, j)` -> polygon at grid position
69-
- `ncells(grid, dim)` -> number of cells in dimension
70-
- `cell_range_extent(grid, irange, jrange)` -> bounding extent for cell range
74+
**Specialized cursors** (`specialized_quadtree_cursors.jl`):
75+
- `IndexOffsetQuadtreeCursor`: For multi-grid scenarios (e.g., cubed spheres), applies an index offset to grid-local indices
76+
- `ReorderedTopDownQuadtreeCursor`: For custom element orderings (space-filling curves in ClimaCore)
7177

72-
**Cursor types** (`src/trees/quadtree_cursors.jl`):
73-
- `QuadtreeCursor`: Cursor with explicit index ranges
74-
- `TopDownQuadtreeCursor`: Top-level cursor wrapping a grid
75-
Sit on top of `AbstractCurvilinearGrid`s to provide quadtree descent.
78+
**Wrappers** (`wrappers.jl`):
79+
- `KnownFullSphereExtentWrapper`: Returns full-sphere extent to skip expensive extent computation
80+
- `CubedSphereToplevelTree`: Vector of per-face cursors with global indexing
7681

77-
**Wrappers** (`src/trees/wrappers.jl`):
78-
- `KnownFullSphereExtentWrapper`: Avoids expensive extent computation for full-sphere trees
82+
**`treeify(manifold, grid)`** (`interfaces.jl`): Dispatches input to the right tree type:
83+
- Matrices of polygons → `TopDownQuadtreeCursor(ExplicitPolygonGrid(...))`
84+
- Matrices of points → `TopDownQuadtreeCursor(CellBasedGrid(...))`
85+
- Iterables of polygons → `FlatNoTree`
86+
- Tuple of vectors → `TopDownQuadtreeCursor(RegularGrid(...))`
87+
- Existing spatial trees → pass-through
7988

8089
### Manifold Support
8190

8291
Grids operate on manifolds from GeometryOps:
8392
- `Planar()`: Cartesian 2D, uses `Extents.Extent` for bounding boxes
8493
- `Spherical()`: Unit sphere (lon/lat), uses `SphericalCap` for bounding regions
8594

86-
The manifold affects extent computation, intersection algorithms, and area calculations.
95+
The manifold affects extent computation, intersection algorithms, and area calculations. If source and destination have different manifolds, it promotes to Spherical.
8796

8897
### Multithreading
8998

90-
`multithreaded_dual_depth_first_search` (`src/utils/MultithreadedDualDepthFirstSearch.jl`) provides parallel tree traversal, spawning tasks when nodes satisfy an area criterion.
91-
92-
### Package Extensions
93-
94-
- **ConservativeRegriddingOceananigansExt**: Oceananigans.jl grid integration
95-
- **ConservativeRegriddingClimaCoreExt**: ClimaCore.jl grid integration (cubed sphere)
96-
- **ConservativeRegriddingInterfacesExt**: Interfaces.jl contracts
99+
`multithreaded_dual_query` (`src/utils/MultithreadedDualDepthFirstSearch.jl`): Parallel dual-tree traversal. Spawns tasks when both nodes are leaves or when nodes satisfy an area criterion (avoids spawning excessive tasks for small regions). Intersection area computation is separately parallelized via ChunkSplitters partitioning.
97100

98-
### Key Dependencies
101+
### Package Extensions (`ext/`)
99102

100-
- **GeometryOps**: Polygon intersection, area calculations, `SpatialTreeInterface`
101-
- **GeoInterface**: Geometry type wrappers
102-
- **SparseArrays**: Sparse matrix storage for intersection weights
103-
- **StableTasks/ChunkSplitters**: Multithreaded computation
103+
All follow the same pattern: implement `Trees.treeify()` for domain-specific grid types.
104104

105-
### Mathematical Model
105+
- **OceananigansExt**: LatitudeLongitudeGrid, RectilinearGrid, TripolarGrid → vertex matrices wrapped in KnownFullSphereExtentWrapper
106+
- **ClimaCoreExt**: Cubed sphere topologies → per-face cursors with index remapping
107+
- **HealpixExt**, **RingGridsExt**: HEALPix and SpeedyWeather grids
108+
- **SpeedyWeatherExt**: Requires *both* RingGrids and SpeedyWeather loaded (dual weak dependency)
109+
- **InterfacesExt**: Interfaces.jl contracts
106110

107-
Given intersection matrix A, source field s, destination field d, and area vectors a_d, a_s:
108-
- Forward: `d = (A * s) / a_d`
111+
### API Surface
109112

110-
Grid cell areas are derived from row/column sums of A when grids fully overlap.
113+
The package uses `@public` from SciMLPublic for API visibility: `Regridder`, `regrid`, `regrid!`, `areas` are public. Grid types and tree types are exported.

Project.toml

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -32,15 +32,16 @@ SpeedyWeather = "9e226e20-d153-4fed-8a5b-493def4f21a9"
3232

3333
[extensions]
3434
ConservativeRegriddingClimaCoreExt = "ClimaCore"
35+
ConservativeRegriddingHealpixExt = "Healpix"
3536
ConservativeRegriddingInterfacesExt = "Interfaces"
3637
ConservativeRegriddingOceananigansExt = "Oceananigans"
3738
ConservativeRegriddingRingGridsExt = "RingGrids"
38-
ConservativeRegriddingHealpixExt = "Healpix"
3939
ConservativeRegriddingSpeedyWeatherExt = ["RingGrids", "SpeedyWeather"]
4040

4141
[compat]
4242
ChunkSplitters = "3"
43-
ClimaCore = "0.14"
43+
ClimaCore = "0.14.49"
44+
Oceananigans = "0.104"
4445
DocStringExtensions = "0.8, 0.9"
4546
Extents = "0.1"
4647
GeoInterface = "1"

docs/Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
[deps]
2+
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
23
Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
34
ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
45
ClimaCore = "d414da3d-4745-48bb-8d80-42e94e092884"
@@ -35,5 +36,4 @@ YAXArrayBase = "90b8fcef-0c2d-428d-9c56-5f86629e9d14"
3536
[sources]
3637
ConservativeRegridding = {path = ".."}
3738
GeometryBasics = {rev = "as/fix-gi-convert-staticarray", url = "https://github.com/JuliaGeometry/GeometryBasics.jl"}
38-
Oceananigans = {rev = "FPivot", url = "https://github.com/briochemc/Oceananigans.jl"}
3939
SphericalSpatialTrees = {rev = "main", url = "https://github.com/meggart/SphericalSpatialTrees.jl"}

examples/isea20_sst.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ function Trees.ncells(tree::SST.ISEACircleTree)
2222
return prod(SST.gridsize(tree))
2323
end
2424

25+
# vec(values)[i] <=> Trees.getcell(tree, i)
2526
function Trees.getcell(tree::SST.ISEACircleTree, i::Int)
2627
return GI.Polygon(SA[GI.LinearRing(SST.index_to_polygon_unitsphere(i, tree))])
2728
end

examples/lowlevel/climacore_exp.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -260,7 +260,7 @@ quadtrees = map(lin2carts, cart2lins, all_coords) do lin2cart, cart2lin, coords
260260
ReorderedTopDownQuadtreeCursor(Trees.CellBasedGrid(GO.Spherical(; radius = mesh.domain.radius), coords), Reorderer2D(cart2lin, lin2cart))
261261
end
262262
quadtrees = map(1:6, all_coords) do face_idx, coords
263-
Trees.FaceAwareQuadtreeCursor(Trees.CellBasedGrid(GO.Spherical(; radius = mesh.domain.radius), coords), face_idx)
263+
Trees.IndexOffsetQuadtreeCursor(Trees.CellBasedGrid(GO.Spherical(; radius = mesh.domain.radius), coords), (face_idx - 1) * ne^2)
264264
end
265265

266266

examples/oceananigans_new_api.jl

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -44,10 +44,13 @@ src_polys = collect(Trees.getcell(Trees.treeify(src_field))) |> x-> GO.transfor
4444
dst_polys = collect(Trees.getcell(Trees.treeify(dst_field))) |> x-> GO.transform(GO.UnitSpherical.GeographicFromUnitSphere(), x) .|> GI.convert(LibGEOS) |> vec
4545

4646
f, a, p = poly(src_polys; color = vec(interior(src_field)), axis = (; type = GlobeAxis), colorrange = extrema(interior(src_field)), highclip = :red)
47-
f, a, p = poly(dst_polys; color = vec(interior(dst_field)), axis = (; type = GlobeAxis), colorrange = extrema(interior(src_field)), highclip = :red)
47+
p2 = poly!(a, src_polys; zlevel = 100_000, color = :transparent, strokewidth = .5, strokecolor = :black, transparency = true)
4848

49-
lines!(a, GeoMakie.coastlines(); zlevel = 100_000, color = :orange)
49+
f, a, p = poly(dst_polys; color = vec(interior(dst_field)), axis = (; type = GlobeAxis), colorrange = extrema(interior(src_field)), highclip = :red)
50+
p2 = poly!(a, dst_polys; zlevel = 100_000, color = :transparent, strokewidth = .5, strokecolor = :black, transparency = true)
5051

52+
lp = lines!(a, GeoMakie.coastlines(); zlevel = 100_000, color = :orange)
53+
lp.color = :gray
5154
# ## Metrics
5255
# First, set up the analytical destination field.
5356
analytical_dst_field = CenterField(dst_grid)

0 commit comments

Comments
 (0)