From b2245f3b254c0e042866bd60edce7b61c90bee64 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Sun, 21 Jun 2026 18:48:10 -0400 Subject: [PATCH 01/20] Bump actions/checkout from 6 to 7 (#423) Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 8d1cc5c5ea..806ba61725 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -33,7 +33,7 @@ jobs: version: '1' arch: arm64 steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - uses: julia-actions/cache@v3 - uses: julia-actions/setup-julia@v3 with: @@ -58,7 +58,7 @@ jobs: statuses: write actions: write steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - uses: julia-actions/cache@v3 - uses: julia-actions/setup-julia@v3 with: @@ -79,7 +79,7 @@ jobs: statuses: write actions: write steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - uses: julia-actions/cache@v3 - uses: julia-actions/setup-julia@v3 with: From dbc5d3794211dac6d999af8fb6edfa25efaf0679 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 30 Jun 2026 14:35:46 -0400 Subject: [PATCH 02/20] Fix spherical Sutherland-Hodgman dropping thin grazing slivers (#425) Co-authored-by: Claude Opus 4.8 (1M context) --- src/methods/clipping/sutherland_hodgman.jl | 14 +++---------- test/methods/clipping/sutherland_hodgman.jl | 22 +++++++++++++++++++++ 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/src/methods/clipping/sutherland_hodgman.jl b/src/methods/clipping/sutherland_hodgman.jl index f93444ee67..a4649074c0 100644 --- a/src/methods/clipping/sutherland_hodgman.jl +++ b/src/methods/clipping/sutherland_hodgman.jl @@ -206,19 +206,11 @@ function _sh_spherical_intersection( intersection_dir = intersection_dir / dir_norm - # Two candidate intersection points (antipodal) + # The two great circles meet at an antipodal pair; keep the crossing on the + # subject arc p1→p2, i.e. the one in the same hemisphere as the midpoint p1+p2. cand1 = UnitSpherical.UnitSphericalPoint{T}(intersection_dir) cand2 = UnitSpherical.UnitSphericalPoint{T}(-intersection_dir) - - # Return the candidate that lies on the subject arc - if UnitSpherical.point_on_spherical_arc(cand1, p1, p2) - return cand1 - elseif UnitSpherical.point_on_spherical_arc(cand2, p1, p2) - return cand2 - end - - # Fallback: neither candidate is on the arc (shouldn't happen if orient detected a crossing) - return UnitSpherical.UnitSphericalPoint{T}(p1) + return cand1 ⋅ p1 + cand1 ⋅ p2 ≥ 0 ? cand1 : cand2 end # Clip polygon against a single edge using Sutherland-Hodgman rules (spherical version) diff --git a/test/methods/clipping/sutherland_hodgman.jl b/test/methods/clipping/sutherland_hodgman.jl index cb038502ec..0f5f13cfec 100644 --- a/test/methods/clipping/sutherland_hodgman.jl +++ b/test/methods/clipping/sutherland_hodgman.jl @@ -357,6 +357,28 @@ import GeoInterface as GI ) @test spherical_area(result) == 0.0 end + + @testset "Grazing sliver overlap is symmetric" begin + # Two real grid cells (OctaHEALPix × HEALPix) overlapping in a thin sliver + # where one corner grazes the other's edge: area must be symmetric. + octa = spherical_polygon([ + (-75.59999999999997, -22.93262592760755), + (-76.15384615384615, -19.86735476489602), + (-79.2, -22.93262592760755), + (-78.75, -25.944479772370016), + ]) + healpix = spherical_polygon([ + (-77.34375, -24.624318352164078), + (-78.04687500000003, -25.282603043311997), + (-77.34375, -25.94447977237001), + (-76.64062500000003, -25.282603043311997), + ]) + alg = GO.ConvexConvexSutherlandHodgman(GO.Spherical()) + a_oh = spherical_area(GO.intersection(alg, octa, healpix)) + a_ho = spherical_area(GO.intersection(alg, healpix, octa)) + @test a_oh ≈ a_ho rtol = 1e-6 + @test 0 < a_oh < 0.1 * spherical_area(healpix) # sliver, not the whole cell + end end end From 57813b972e25ebb1fcb340aa1ba92238dff13478 Mon Sep 17 00:00:00 2001 From: Will Bodeau Date: Tue, 30 Jun 2026 12:03:36 -0700 Subject: [PATCH 03/20] Tutorial revisions (#415) Co-authored-by: Anshul Singhvi Co-authored-by: Alex Gardner --- docs/Project.toml | 1 + docs/make.jl | 1 + docs/src/tutorials/creating_geometry.md | 295 +++++----------------- docs/src/tutorials/geodesic_paths.md | 17 -- docs/src/tutorials/geospatial_geometry.md | 278 ++++++++++++++++++++ docs/src/tutorials/spatial_joins.md | 29 ++- 6 files changed, 361 insertions(+), 260 deletions(-) delete mode 100644 docs/src/tutorials/geodesic_paths.md create mode 100644 docs/src/tutorials/geospatial_geometry.md diff --git a/docs/Project.toml b/docs/Project.toml index b18a68afcc..d51bc58119 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -40,6 +40,7 @@ MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Proj = "c94c279d-25a6-4763-9509-64d165bea63e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" diff --git a/docs/make.jl b/docs/make.jl index 571228e34c..2e63c523d5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -134,6 +134,7 @@ makedocs(; "API Reference" => "api.md", "Tutorials" => [ "Creating Geometry" => "tutorials/creating_geometry.md", + "Geospatial Geometry" => "tutorials/geospatial_geometry.md", "Spatial Joins" => "tutorials/spatial_joins.md", ], "Explanations" => [ diff --git a/docs/src/tutorials/creating_geometry.md b/docs/src/tutorials/creating_geometry.md index f248587239..9da04e2f15 100644 --- a/docs/src/tutorials/creating_geometry.md +++ b/docs/src/tutorials/creating_geometry.md @@ -1,25 +1,27 @@ # Creating Geometry -In this tutorial, we're going to: -1. [Create and plot geometries](@ref creating-geometry) -2. [Plot geometries on a map using `GeoMakie` and coordinate reference system (`CRS`)](@ref plot-geometry) -3. [Create geospatial geometries with embedded coordinate reference system information](@ref geom-crs) -4. [Assign attributes to geospatial geometries](@ref attributes) -5. [Save geospatial geometries to common geospatial file formats](@ref save-geometry) +In this tutorial, we're going to build some basic 2D geometry. +This [follows the Simple Features hierarchy for geospatial geometry](https://juliageo.org/GeoInterface.jl/stable/background/sf/): +1. [Create `Point`s and `MultiPoint`s](#create-points-and-multipoints) +2. [Connect `Point`s into `LineString`s](#connecting-points-into-lines) +3. [Build `LinearRing`s, `Polygon`s, and `MultiPolygon`s](#building-polygons-and-multipolygons) +Install the packages used in this tutorial: -First, we load some required packages. +````julia +using Pkg +Pkg.add(["GeoInterface", "GeometryOps", + "CoordinateTransformations", + "CairoMakie", "GeoMakie"]) +```` ````@example creating_geometry # Geospatial packages from Julia import GeoInterface as GI import GeometryOps as GO -import GeoFormatTypes as GFT -using GeoJSON # to load some data -# Packages for coordinate transformation and projection +# Coordinate transformation and projection import CoordinateTransformations -import Proj # Plotting using CairoMakie using GeoMakie @@ -27,7 +29,7 @@ using DisplayAs # hide Makie.set_theme!(Makie.MAKIE_DEFAULT_THEME) # hide ```` -## [Creating and plotting geometries](@id creating-geometry) +## Create `Point`s and `MultiPoint`s Let's start by making a single `Point`. @@ -44,9 +46,11 @@ fig, ax, plt = plot(point) Let's create a set of points, and have a bit more fun with plotting. ````@example creating_geometry -x = [-5, 0, 5, 0]; -y = [0, -5, 0, 5]; -points = GI.Point.(zip(x,y)); +xs = [-5, 0, 5, 0] +ys = [0, -5, 0, 5] + +points = GI.Point.(xs, ys) + plot!(ax, points; marker = '✈', markersize = 30) fig ```` @@ -54,21 +58,27 @@ fig `Point`s can be combined into a single `MultiPoint` geometry. ````@example creating_geometry -x = [-5, -5, 5, 5]; -y = [-5, 5, 5, -5]; -multipoint = GI.MultiPoint(GI.Point.(zip(x, y))); +xs = [-5, -5, 5, 5] +ys = [-5, 5, 5, -5] + +# zip: Create (x, y) coordinates (tuples) +# GI.Point: Turn each coordinate pair into special Point geometries +# GI.MultiPoint: Wrap all Points into a single MultiPoint geometry object +multipoint = GI.MultiPoint(GI.Point.(xs, ys)); # TODO: GeoInterfaceMakie.jl can't plot multipoints due to breaking changes # in Makie.jl. We should fix that. plot!(ax, multipoint.geom; marker = '☁', markersize = 30) fig ```` +## Connecting points into lines + Let's create a `LineString` connecting two points. ````@example creating_geometry -p1 = GI.Point.(-5, 0); -p2 = GI.Point.(5, 0); -line = GI.LineString([p1,p2]) +p1 = GI.Point.(-5, 0) +p2 = GI.Point.(5, 0) +line = GI.LineString([p1, p2]) plot!(ax, line; color = :red) fig ```` @@ -77,16 +87,20 @@ Now, let's create a line connecting multiple points (i.e. a `LineString`). This time we get a bit more fancy with point creation. ````@example creating_geometry -r = 2; -k = 10; -ϴ = 0:0.01:2pi; -x = r .* (k + 1) .* cos.(ϴ) .- r .* cos.((k + 1) .* ϴ); -y = r .* (k + 1) .* sin.(ϴ) .- r .* sin.((k + 1) .* ϴ); -lines = GI.LineString(GI.Point.(zip(x,y))); +r = 2 +k = 10 +ϴs = 0:0.01:2pi +xs = r .* (k + 1) .* cos.(ϴs) .- r .* cos.((k + 1) .* ϴs) +ys = r .* (k + 1) .* sin.(ϴs) .- r .* sin.((k + 1) .* ϴs) + +lines = GI.LineString(GI.Point.(xs, ys)) + plot!(ax, lines; linewidth = 5) fig ```` +## Building polygons and multipolygons + We can also create a single `LinearRing` trait, the building block of a polygon. A `LinearRing` is simply a `LineString` with the same beginning and endpoint, i.e., an arbitrary closed shape composed of point pairs. @@ -99,16 +113,17 @@ ring1 = GI.LinearRing(GI.getpoint(lines)); Now, let's make the `LinearRing` into a `Polygon`. ````@example creating_geometry +# Polygon fills the interior of a LinearRing, turning it into a solid shape polygon1 = GI.Polygon([ring1]); ```` Now, we can use GeometryOps and CoordinateTransformations to shift `polygon1` up, to avoid plotting over our earlier results. This is done through the [GeometryOps.transform](@ref) function. ````@example creating_geometry -xoffset = 0.; -yoffset = 50.; -f = CoordinateTransformations.Translation(xoffset, yoffset); -polygon1 = GO.transform(f, polygon1); +xoffset = 0. +yoffset = 50. +f = CoordinateTransformations.Translation(xoffset, yoffset) +polygon1 = GO.transform(f, polygon1) plot!(polygon1) fig ```` @@ -125,219 +140,37 @@ polygon2 = GI.Polygon([ring1, hole]) Shift `polygon2` to the right, to avoid plotting over our earlier results. ````@example creating_geometry -xoffset = 50.; -yoffset = 0.; -f = CoordinateTransformations.Translation(xoffset, yoffset); -polygon2 = GO.transform(f, polygon2); +xoffset = 50. +yoffset = 0. +f = CoordinateTransformations.Translation(xoffset, yoffset) +polygon2 = GO.transform(f, polygon2) plot!(polygon2) fig ```` -`Polygon`s can also be grouped together as a `MultiPolygon`. +Similar to `Point`s with `MultiPoint`s, `Polygon`s can also be grouped together as a `MultiPolygon`. ````@example creating_geometry -r = 5; -x = cos.(reverse(ϴ)) .* r .+ xoffset; -y = sin.(reverse(ϴ)) .* r .+ yoffset; -ring2 = GI.LinearRing(GI.Point.(zip(x,y))); -polygon3 = GI.Polygon([ring2]); +# Create a simple circle with a radius of 5 +r = 5 +xs = cos.(reverse(ϴs)) .* r .+ xoffset +ys = sin.(reverse(ϴs)) .* r .+ yoffset +ring2 = GI.LinearRing(GI.Point.(xs, ys)) +polygon3 = GI.Polygon([ring2]) + +# Group polygon2 (our shape with the square hole) and polygon3 (our circle) together into a MultiPolygon multipolygon = GI.MultiPolygon([polygon2, polygon3]) ```` Shift `multipolygon` up, to avoid plotting over our earlier results. ````@example creating_geometry -xoffset = 0.; -yoffset = 50.; -f = CoordinateTransformations.Translation(xoffset, yoffset); -multipolygon = GO.transform(f, multipolygon); +xoffset = 0. +yoffset = 50. +f = CoordinateTransformations.Translation(xoffset, yoffset) +multipolygon = GO.transform(f, multipolygon) plot!(multipolygon) fig ```` -Great, now we can make `Points`, `MultiPoints`, `Lines`, `LineStrings`, `Polygons` (with holes), and `MultiPolygons` and modify them using [`CoordinateTransformations`] and [`GeometryOps`]. - -## [Plot geometries on a map using `GeoMakie` and coordinate reference system (`CRS`)](@id plot-geometry) - -In geospatial sciences we often have data in one [Coordinate Reference System (CRS)](https://en.wikipedia.org/wiki/Spatial_reference_system) (`source`) and would like to display it in different (`destination`) `CRS`. `GeoMakie` allows us to do this by automatically projecting from `source` to `destination` CRS. - -Here, our `source` CRS is common geographic (i.e. coordinates of latitude and longitude), [WGS84](https://epsg.io/4326). - -````@example creating_geometry -source_crs1 = GFT.EPSG(4326) -```` - -Now let's pick a `destination` CRS for displaying our map. Here we'll pick [natearth2](https://proj.org/en/9.4/operations/projections/natearth2.html). - -````@example creating_geometry -destination_crs = "+proj=natearth2" -```` - -Let's add land area for context. First, download and open the [Natural Earth](https://www.naturalearthdata.com) global land polygons at 110 m resolution.`GeoMakie` ships with this particular dataset, so we will access it from there. - -````@example creating_geometry -land_path = GeoMakie.assetpath("ne_110m_land.geojson") -```` - -!!! note - Natural Earth has lots of other datasets, and there is a Julia package that provides an interface to it called [NaturalEarth.jl](https://github.com/JuliaGeo/NaturalEarth.jl). - -Read the land `MultiPolygon`s as a `GeoJSON.FeatureCollection`. - -````@example creating_geometry -land_geo = GeoJSON.read(land_path) -```` - -We then need to create a figure with a `GeoAxis` that can handle the projection between `source` and `destination` CRS. For GeoMakie, `source` is the CRS of the input and `dest` is the CRS you want to visualize in. - -````@example creating_geometry -fig = Figure(size=(1000, 500)); -ga = GeoAxis( - fig[1, 1]; - source = source_crs1, - dest = destination_crs, - xticklabelsvisible = false, - yticklabelsvisible = false, -); -nothing #hide -```` - -Plot `land` for context. - -````@example creating_geometry -poly!(ga, land_geo, color=:black) -fig -```` - -Now let's plot a `Polygon` like before, but this time with a CRS that differs from our `source` data - -````@example creating_geometry -plot!(multipolygon; color = :green) -fig -```` - -But what if we want to plot geometries with a different `source` CRS on the same figure? - -To show how to do this let's create a geometry with coordinates in UTM (Universal Transverse Mercator) zone 10N [EPSG:32610](https://epsg.io/32610). - -````@example creating_geometry -source_crs2 = GFT.EPSG(32610) -```` - -Create a polygon (we're working in meters now, not latitude and longitude) -````@example creating_geometry -r = 1000000; -ϴ = 0:0.01:2pi; -x = r .* cos.(ϴ).^3 .+ 500000; -y = r .* sin.(ϴ) .^ 3 .+5000000; -DisplayAs.setcontext(y, :compact => true, :displaysize => (10, 40),) # hide -```` - -Now create a `LinearRing` from `Points` - -````@example creating_geometry -ring3 = GI.LinearRing(Point.(zip(x, y))) -```` - -Now create a `Polygon` from the `LineRing` - -````@example creating_geometry -polygon3 = GI.Polygon([ring3]) -```` - -Now plot on the existing GeoAxis. - -!!! note - The keyword argument `source` is used to specify the source `CRS` of that particular plot, when plotting on an existing `GeoAxis`. - -````@example creating_geometry -plot!(ga,polygon3; color=:red, source = source_crs2) -fig -```` - -## [Create geospatial geometries with embedded coordinate reference system information](@id geom-crs) - -Great, we can make geometries and plot them on a map... now let's export the data to common geospatial data formats. To do this we now need to create geometries with embedded `CRS` information, making it a geospatial geometry. All that's needed is to include `; crs = crs` as a keyword argument when constructing the geometry. - -Let's do this for a new `Polygon` -````@example creating_geometry -r = 3; -k = 7; -ϴ = 0:0.01:2pi; -x = r .* (k + 1) .* cos.(ϴ) .- r .* cos.((k + 1) .* ϴ); -y = r .* (k + 1) .* sin.(ϴ) .- r .* sin.((k + 1) .* ϴ); -ring4 = GI.LinearRing(Point.(zip(x, y))) -```` - -But this time when we create the `Polygon` we need to specify the `CRS` at the time of creation, making it a geospatial polygon - -````@example creating_geometry -geopoly1 = GI.Polygon([ring4], crs = source_crs1) -```` - -!!! note - It is good practice to only include CRS information with the highest-level geometry. Not doing so can bloat the memory footprint of the geometry. CRS information _can_ be included at the individual `Point` level but is discouraged. - -And let's create second `Polygon` by shifting the first using CoordinateTransformations - -````@example creating_geometry -xoffset = 20.; -yoffset = -25.; -f = CoordinateTransformations.Translation(xoffset, yoffset); -geopoly2 = GO.transform(f, geopoly1); -```` - -## [Creating a table with attributes and geometry](@id attributes) - -Typically, you'll also want to include attributes with your geometries. Attributes are simply data that are attributed to each geometry. The easiest way to do this is to create a table with a `:geometry` column. Let's do this using [`DataFrames`](https://github.com/JuliaData/DataFrames.jl). - -````@example creating_geometry -using DataFrames -df = DataFrame(geometry=[geopoly1, geopoly2]) -```` - -Now let's add a couple of attributes to the geometries. We do this using [DataFrames' `!` mutation syntax](https://dataframes.juliadata.org/stable/man/getting_started/#The-DataFrame-Type) that allows you to add a new column to an existing data frame. - -````@example creating_geometry -df[!,:id] = ["a", "b"] -df[!, :name] = ["polygon 1", "polygon 2"] -df -```` - -## [Saving your geospatial data](@id save-geometry) - -There are Julia packages for most commonly used geographic data formats. Below, we show how to export that data to each of these. - -We begin with [GeoJSON](https://github.com/JuliaGeo/GeoJSON.jl), which is a [JSON](https://en.wikipedia.org/wiki/JSON) format for geospatial feature collections. It's human-readable and widely supported by most web-based and desktop geospatial libraries. - -````@example creating_geometry -import GeoJSON -fn = "shapes.json" -GeoJSON.write(fn, df) -```` - -Now, let's save as a [`Shapefile`](https://github.com/JuliaGeo/Shapefile.jl). Shapefiles are actually a set of files (usually 4) that hold geometry information, a CRS, and additional attribute information as a separate table. When you give `Shapefile.write` a file name, it will write 4 files of the same name but with different extensions. - -````@example creating_geometry -import Shapefile -fn = "shapes.shp" -Shapefile.write(fn, df) -```` - -Now, let's save as a [`GeoParquet`](https://github.com/JuliaGeo/GeoParquet.jl). GeoParquet is a geospatial extension to the [Parquet](https://parquet.apache.org/) format, which is a high-performance data store. It's great for storing large amounts of data in a single file. - -````@example creating_geometry -import GeoParquet -fn = "shapes.parquet" -GeoParquet.write(fn, df, (:geometry,)) -```` - -Finally, if there's no Julia-native package that can write data to your desired format (e.g. `.gpkg`, `.gml`, etc), you can use [`GeoDataFrames`](https://github.com/evetion/GeoDataFrames.jl). This package uses the [GDAL](https://gdal.org/) library under the hood which supports writing to nearly all geospatial formats. - -````@example creating_geometry -import GeoDataFrames -fn = "shapes.gpkg" -GeoDataFrames.write(fn, df) -```` - -And there we go, you can now create mapped geometries from scratch, manipulate them, plot them on a map, and save them in multiple geospatial data formats. +Great, now we can make `Point`s, `MultiPoint`s, `Line`s, `LineString`s, `Polygon`s (with holes), and `MultiPolygon`s and modify them using [`CoordinateTransformations`] and [`GeometryOps`]. diff --git a/docs/src/tutorials/geodesic_paths.md b/docs/src/tutorials/geodesic_paths.md deleted file mode 100644 index f6f4bf6531..0000000000 --- a/docs/src/tutorials/geodesic_paths.md +++ /dev/null @@ -1,17 +0,0 @@ -# Geodesic paths - -Geodesic paths are paths computed on an ellipsoid, as opposed to a plane. - -```@example geodesic -import GeometryOps as GO, GeoInterface as GI -using CairoMakie, GeoMakie - - -IAH = (-95.358421, 29.749907) -AMS = (4.897070, 52.377956) - - -fig, ga, _cp = lines(GeoMakie.coastlines(); axis = (; type = GeoAxis)) -lines!(ga, GO.segmentize(GO.GeodesicSegments(; max_distance = 100_000), GI.LineString([IAH, AMS])); color = Makie.wong_colors()[2]) -fig -``` \ No newline at end of file diff --git a/docs/src/tutorials/geospatial_geometry.md b/docs/src/tutorials/geospatial_geometry.md new file mode 100644 index 0000000000..d62ea5f846 --- /dev/null +++ b/docs/src/tutorials/geospatial_geometry.md @@ -0,0 +1,278 @@ +# Geospatial Geometry + +In this tutorial, we're going to: +1. [Plot geometries on a map using `GeoMakie` and coordinate reference system (`CRS`)](@ref plot-geometry) +2. [Create geospatial geometries with embedded coordinate reference system information](@ref geom-crs) +3. [Assign attributes to geospatial geometries](@ref attributes) +4. [Save geospatial geometries to common geospatial file formats](@ref save-geometry) +5. [Introduce Geodesic Paths](@ref geodesic-paths) + +Install the packages used in this tutorial: + +````julia +using Pkg +Pkg.add(["GeoInterface", "GeometryOps", "GeoFormatTypes", + "GeoJSON", "GeoParquet", "GeoDataFrames", + "CoordinateTransformations", "Proj", "DataFrames", + "CairoMakie", "GeoMakie", "Shapefile"]) +```` + +````@example geospatial_geometry +# Geospatial packages +import GeoInterface as GI +import GeometryOps as GO +import GeoFormatTypes as GFT +# Coordinate transformation and projection +import CoordinateTransformations +import Proj +# Plotting +using CairoMakie +using GeoMakie +using DisplayAs # hide +Makie.set_theme!(Makie.MAKIE_DEFAULT_THEME) # hide +# Data loading +using GeoDataFrames, DataFrames +import GeoParquet, GeoJSON, Shapefile # activate native IO packages +```` + + +## [Plot geometries on a map using `GeoMakie` and coordinate reference system (`CRS`)](@id plot-geometry) + +In geospatial sciences we often have data in one [Coordinate Reference System (CRS)](https://en.wikipedia.org/wiki/Spatial_reference_system) (`source`) and would like to display it in different (`destination`) `CRS`. `GeoMakie` allows us to do this by automatically projecting from `source` to `destination` CRS. + +Here, our `source` CRS is common geographic (i.e. coordinates of latitude and longitude), [WGS84](https://epsg.io/4326). + +````@example geospatial_geometry +source_crs1 = GFT.EPSG(4326) +```` + +Now let's pick a `destination` CRS for displaying our map. Here we'll pick [natearth2](https://proj.org/en/9.4/operations/projections/natearth2.html). + +````@example geospatial_geometry +destination_crs = "+proj=natearth2" +```` + +Let's add land area for context. First, download and open the [Natural Earth](https://www.naturalearthdata.com) global land polygons at 110 m resolution. `GeoMakie` ships with this particular dataset, so we will access it from there. +Note that this will be a path on your local machine, so you could easily point it to any other `.geojson` file you have. + +````@example geospatial_geometry +land_geojson_path = GeoMakie.assetpath("ne_110m_land.geojson") +```` + +!!! note + Natural Earth has lots of other datasets, and there is a Julia package that provides an interface to it called [NaturalEarth.jl](https://github.com/JuliaGeo/NaturalEarth.jl). + +Read this dataset in as a `GeoJSON.FeatureCollection`. + +````@example geospatial_geometry +land_features = GeoDataFrames.read(land_geojson_path) +```` + +We then need to create a figure with a `GeoAxis` that can handle the projection between `source` and `destination` CRS. For [`GeoMakie`](https://geo.makie.org/stable/), `source` is the CRS of the input and `dest` is the CRS you want to visualize in. + +````@example geospatial_geometry +fig = Figure(size=(1000, 500)); +ga = GeoAxis( + fig[1, 1]; + source = source_crs1, + dest = destination_crs, + xticklabelsvisible = false, + yticklabelsvisible = false, +); +nothing #hide +```` + +Plot `land` for context. + +````@example geospatial_geometry +poly!(ga, land_features; color=:black) +fig +```` + +Now let's plot a `Polygon` like before, but this time with a CRS that differs from our `source` data. + +````@example geospatial_geometry +function ring(radius) + ϴ = 0:0.01:2π + points = @. GI.Point(50 + radius * cos(ϴ), 50 + radius * sin(ϴ)) + return GI.LinearRing(points) +end + +function spiro(a = 22, b = 2, k = 11) + ϴ = 0:0.01:2π + points = @. GI.Point( + 50 + a*cos(ϴ) - b*cos(k*ϴ), + 50 + a*sin(ϴ) - b*sin(k*ϴ) + ) + return GI.LinearRing(points) +end + + +multipolygon = GI.MultiPolygon([ + GI.Polygon([spiro(), ring(8)]), + GI.Polygon([ring(4)]) +]) + +plot!(multipolygon; color = :green) +fig +```` + +But what if we want to plot geometries with a different `source` CRS on the same figure? + +To show how to do this let's create a geometry with coordinates in UTM (Universal Transverse Mercator) zone 10N [EPSG:32610](https://epsg.io/32610). + +````@example geospatial_geometry +source_crs2 = GFT.EPSG(32610) +```` + +Create a polygon (we're working in meters now, not latitude and longitude) +````@example geospatial_geometry +rs = 1000000; +ϴs = 0:0.01:2pi; +xs = rs .* cos.(ϴs).^3 .+ 500000; +ys = rs .* sin.(ϴs) .^ 3 .+ 5000000; +DisplayAs.setcontext(y, :compact => true, :displaysize => (10, 40),) # hide +```` + +Now create a `LinearRing` from `Points` + +````@example geospatial_geometry +ring3 = GI.LinearRing(GI.Point.(xs, ys)) +```` + +Now create a `Polygon` from the `LineRing` + +````@example geospatial_geometry +polygon3 = GI.Polygon([ring3]) +```` + +Now plot on the existing GeoAxis. + +!!! note + The keyword argument `source` is used to specify the source `CRS` of that particular plot, when plotting on an existing `GeoAxis`. + +````@example geospatial_geometry +plot!(ga,polygon3; color=:red, source = source_crs2) +fig +```` + +## [Create geospatial geometries with embedded coordinate reference system information](@id geom-crs) + +Great, we can make geometries and plot them on a map... now let's export the data to common geospatial data formats. To do this we now need to create geometries with embedded `CRS` information, making it a geospatial geometry. All that's needed is to include `; crs = crs` as a keyword argument when constructing the geometry. + +Let's do this for a new `Polygon` +````@example geospatial_geometry +rs = 3; +ks = 7; +ϴs = 0:0.01:2pi; +xs = rs .* (ks + 1) .* cos.(ϴs) .- rs .* cos.((ks + 1) .* ϴs); +ys = rs .* (ks + 1) .* sin.(ϴs) .- rs .* sin.((ks + 1) .* ϴs); +ring4 = GI.LinearRing(tuple.(xs, ys)) +```` + +But this time when we create the `Polygon` we need to specify the `CRS` at the time of creation, associating that metadata with the geometry. + +````@example geospatial_geometry +geopoly1 = GI.Polygon([ring4], crs = source_crs1) +```` + +!!! note + It is good practice to only include CRS information with the highest-level geometry. Not doing so can sometimes bloat the memory footprint of the geometry. CRS information _can_ be included at the individual `Point` level but is discouraged. + +And let's create a second `Polygon` by shifting the first using CoordinateTransformations + +````@example geospatial_geometry +xoffset = 20. +yoffset = -25. +f = CoordinateTransformations.Translation(xoffset, yoffset) +geopoly2 = GO.transform(f, geopoly1) +```` + +## [Creating a table with attributes and geometry](@id attributes) + +Typically, you'll also want to include attributes with your geometries. Attributes are simply data that are attributed to each geometry. The easiest way to do this is to create a table with a `:geometry` column. Let's do this using [`DataFrames`](https://github.com/JuliaData/DataFrames.jl). + +````@example geospatial_geometry +df = DataFrame(geometry = [geopoly1, geopoly2]) +```` + +Now let's add a couple of attributes to the geometries by adding new columns to our existing data frame. + +````@example geospatial_geometry +df.id = ["a", "b"] +df.name = ["polygon 1", "polygon 2"] +df +```` + +## [Saving your geospatial data](@id save-geometry) + +There are Julia packages for most commonly used geographic data formats. Below, we show how to export that data to each of these. + +::: tabs + +== GeoDataFrames + +In general, [`GeoDataFrames`](https://github.com/evetion/GeoDataFrames.jl) is recomended as the default way to write data to your desired format. This package uses the [GDAL](https://gdal.org/) library under the hood which supports writing to nearly all geospatial formats. + +Writing to `gpkg`: + +````@example geospatial_geometry +GeoDataFrames.write("shapes.gpkg", df) +```` + +Writing to `GeoJSON`: + +````@example geospatial_geometry +GeoDataFrames.write("file.geojson", df) +```` + +View the [`GeoDataFrames`](https://github.com/evetion/GeoDataFrames.jl) documentation for all recognized file extensions. + + +== GeoJSON + +Next, let's save as a [GeoJSON](https://github.com/JuliaGeo/GeoJSON.jl), which is a [JSON](https://en.wikipedia.org/wiki/JSON) format for geospatial feature collections. It's human-readable and widely supported by most web-based and desktop geospatial libraries. + +````@example geospatial_geometry +GeoJSON.write("shapes.json", df) +```` + +== Shapefile + +Now, let's save as a [`Shapefile`](https://github.com/JuliaGeo/Shapefile.jl). Shapefiles are actually a set of files (usually 4) that hold geometry information, a CRS, and additional attribute information as a separate table. When you give `Shapefile.write` a file name, it will write 4 files of the same name but with different extensions. + +````@example geospatial_geometry +Shapefile.write("shapes.shp", df) +```` + +== GeoParquet + +Now, let's save as a [`GeoParquet`](https://github.com/JuliaGeo/GeoParquet.jl). GeoParquet is a geospatial extension to the [Parquet](https://parquet.apache.org/) format, which is a high-performance data store. It's great for storing large amounts of data in a single file. + +````@example geospatial_geometry +GeoParquet.write("shapes.parquet", df) +```` + +::: + +## [Geodesic paths](@id geodesic-paths) + +Geodesic paths are paths computed on an ellipsoid, as opposed to a plane. The geodesic is the shortest path between two points measured along the Earth's curved surface. Because the surface is curved, that shortest path appears as a curve (not a straight line) when drawn on a flat map. + +Here, we use the `segmentize` function to add vertices along the geodesic between two points, so the line follows Earth's curved surface instead of a straight line. + +````@example geospatial_geometry +# Two points in (longitude, latitude) order: Houston (IAH) and Amsterdam (AMS). +# Geodesic methods assume lon/lat input. +IAH = (-95.358421, 29.749907) +AMS = (4.897070, 52.377956) + +# Draw coastlines for geographic context +fig, ga, _cp = lines(GeoMakie.coastlines(); axis = (; type = GeoAxis)) + +# Create our line along the Earth, accounting for curvature +lines!(ga, GO.segmentize(GO.Geodesic(), GI.LineString([IAH, AMS]); max_distance = 100_000)) +fig +```` + +And there we go, you can now create mapped geometries from scratch, manipulate them, plot them on a map, and save them in multiple geospatial data formats, as well as create geodesic paths. diff --git a/docs/src/tutorials/spatial_joins.md b/docs/src/tutorials/spatial_joins.md index 2bb63f3f86..143be81b39 100644 --- a/docs/src/tutorials/spatial_joins.md +++ b/docs/src/tutorials/spatial_joins.md @@ -4,11 +4,26 @@ Spatial joins are [table joins](https://www.geeksforgeeks.org/sql-join-set-1-inn Spatial joins can be done between any geometry types (from geometrycollections to points), just as geometrical predicates can be evaluated on any geometries. -In this tutorial, we will show how to perform a spatial join on first a toy dataset and then two Natural Earth datasets, to show how this can be used in the real world. +In this tutorial, we will demonstrate how to perform spatial joins, first using a toy dataset and then two Natural Earth datasets, to show how this can be applied in the real world + +Install the packages used in this tutorial: + +````julia +using Pkg +Pkg.add(["FlexiJoins", "DataFrames", "CairoMakie", "GeoInterfaceMakie", "GADM", "GeoInterface", "GeometryOps"]) +```` + +````@example spatialjoins +using FlexiJoins, +using DataFrames +using CairoMakie, GeoInterfaceMakie +using GADM # GADM gives us country and sublevel geometry +import GeoInterface as GI +import GeometryOps as GO +```` In order to perform the spatial join, we use **[FlexiJoins.jl](https://github.com/JuliaAPlavin/FlexiJoins.jl)** to perform the join, specifically using its `by_pred` joining method. This allows the user to specify a predicate in the following manner, for any kind of table join operation: ```julia -using FlexiJoins innerjoin((table1, table1), by_pred(:table1_column, predicate_function, :table2_column) # & add other conditions here ) @@ -43,11 +58,6 @@ The spatial join is performed using the `contains` predicate from GeometryOps, w First, we generate our data. We create two triangle polygons which, together, span the rectangle (0, 0, 1, 1), and a set of points which are randomly distributed within this rectangle. ```@example spatialjoins -import GeoInterface as GI, GeometryOps as GO -using FlexiJoins, DataFrames - -using CairoMakie, GeoInterfaceMakie - pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])]) pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])]) poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue]) @@ -92,11 +102,6 @@ Here, you can see that the colors were assigned appropriately to the scattered p Suppose I have a list of polygons representing administrative regions (or mining sites, or what have you), and I have a list of polygons for each country. I want to find the country each region is in. ```julia real -import GeoInterface as GI, GeometryOps as GO -using FlexiJoins, DataFrames, GADM # GADM gives us country and sublevel geometry - -using CairoMakie, GeoInterfaceMakie - country_df = GADM.get.(["JPN", "USA", "IND", "DEU", "FRA"]) |> DataFrame country_df.geometry = GI.GeometryCollection.(GO.tuples.(country_df.geom)) From 718544db632b74d0f31579cfdf985952ae7bc8f2 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 30 Jun 2026 15:04:00 -0400 Subject: [PATCH 04/20] Bump patch version to v0.1.41 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 35c043dfa2..104eab3d34 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeometryOps" uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" authors = ["Anshul Singhvi ", "Rafael Schouten ", "Skylar Gering ", "and contributors"] -version = "0.1.40" +version = "0.1.41" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" From cb4cb5b55f6b46146c53f03341699114a5d9f9d3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 4 Jul 2026 07:21:19 -0700 Subject: [PATCH 05/20] Fix off-by-one indexing bug in DouglasPeucker simplify (#387) Co-authored-by: Claude Opus 4.5 Co-authored-by: Claude Opus 4.8 Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> --- src/transformations/simplify.jl | 2 +- test/transformations/simplify.jl | 11 +++++++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/src/transformations/simplify.jl b/src/transformations/simplify.jl index acd9840d36..2de86b4f2e 100644 --- a/src/transformations/simplify.jl +++ b/src/transformations/simplify.jl @@ -301,7 +301,7 @@ function _simplify(alg::DouglasPeucker, points::Vector, preserve_endpoint) i = 2 # already have first and last point added start_idx, end_idx = 1, npoints max_idx, max_dist = _find_max_squared_dist(points, start_idx, end_idx) - while i ≤ min(MIN_POINTS + 1, max_points) || (i < max_points && max_dist > max_tol) + while i < min(MIN_POINTS + 1, max_points) || (i < max_points && max_dist > max_tol) # Add next point to results i += 1 results[i] = max_idx diff --git a/test/transformations/simplify.jl b/test/transformations/simplify.jl index e20e9684ce..8a6bacaa6a 100644 --- a/test/transformations/simplify.jl +++ b/test/transformations/simplify.jl @@ -21,6 +21,17 @@ datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data")) end @testset "DouglasPeucker" begin + # Test for issue #386: BoundsError when simplifying small geometries with low number/ratio + @testset "small geometry simplification (issue #386)" begin + line = GI.LineString([(0.0, 0.0), (1.0, 0.5), (2.0, 0.0), (3.0, 1.0)]) + @test_nowarn GO.simplify(line; ratio=0.1) + @test_nowarn GO.simplify(line; tol=0.1) + @test_nowarn GO.simplify(line; number=3) + # Verify the output is valid + result = GO.simplify(line; number=3) + @test GI.npoint(result) == 3 + end + poly_coords = JLD2.jldopen(joinpath(datadir, "complex_polygons.jld2"))["verts"][1:4] for c in poly_coords npoints = length(c[1]) From 95c975315dec90a9bafa15cbb5d3e75a0f5b00ec Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 5 Jul 2026 13:15:31 -0700 Subject: [PATCH 06/20] Add DimensionalData extension for `polygonize` (#429) Co-authored-by: Claude Opus 4.8 Co-authored-by: rafaqz --- Project.toml | 3 ++ .../GeometryOpsDimensionalDataExt.jl | 22 ++++++++ src/methods/polygonize.jl | 26 +++++---- test/methods/polygonize.jl | 54 +++++++++++++------ test/runtests.jl | 1 + 5 files changed, 79 insertions(+), 27 deletions(-) create mode 100644 ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl diff --git a/Project.toml b/Project.toml index 104eab3d34..fde30baf29 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [weakdeps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37" LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" @@ -34,6 +35,7 @@ GeometryOpsCore = {path = "GeometryOpsCore"} [extensions] GeometryOpsDataFramesExt = "DataFrames" +GeometryOpsDimensionalDataExt = "DimensionalData" GeometryOpsFlexiJoinsExt = "FlexiJoins" GeometryOpsLibGEOSExt = "LibGEOS" GeometryOpsMakieExt = "Makie" @@ -47,6 +49,7 @@ CoordinateTransformations = "0.5, 0.6" DataAPI = "1" DataFrames = "1" DelaunayTriangulation = "1.0.4" +DimensionalData = "0.27, 0.28, 0.29, 0.30" ExactPredicates = "2.2.8" Extents = "0.1.5" FlexiJoins = "0.1.30" diff --git a/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl b/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl new file mode 100644 index 0000000000..1affd95acb --- /dev/null +++ b/ext/GeometryOpsDimensionalDataExt/GeometryOpsDimensionalDataExt.jl @@ -0,0 +1,22 @@ +module GeometryOpsDimensionalDataExt + +import DimensionalData as DD +import GeometryOps as GO +import GeoInterface as GI + +# Polygonize an `AbstractDimArray` in the coordinate space of its `X`/`Y` lookups +# (via their interval bounds) rather than the raw integer axes. +function GO.polygonize(A::DD.AbstractDimArray; dims=(DD.X(), DD.Y()), crs=GI.crs(A), kw...) + lookups = DD.lookup(A, dims) + bounds_vecs = if all(DD.isintervals, lookups) + map(DD.intervalbounds, lookups) + else + @warn "`polygonize` is not strictly correct for `Points` sampling, as polygons cover space by definition. Treating as `Intervals`, but this may not be appropriate." + map(lookups) do l + DD.intervalbounds(DD.set(l, DD.Intervals())) + end + end + return GO.polygonize(bounds_vecs..., A; crs, kw...) +end + +end diff --git a/src/methods/polygonize.jl b/src/methods/polygonize.jl index c725a9b6c6..e2c53b7536 100644 --- a/src/methods/polygonize.jl +++ b/src/methods/polygonize.jl @@ -377,21 +377,25 @@ end function _pixel_edges(f, xs::AbstractVector{T}, ys::AbstractVector{T}, A) where T<:Tuple edges = Dict{T,Tuple{T,T}}() - # First we collect all the edges around target pixels - fi, fj = map(first, axes(A)) - li, lj = map(last, axes(A)) - for j in axes(A, 2) - y1, y2 = ys[j] - for i in axes(A, 1) + # First we collect all the edges around target pixels. + # `i`/`j` index rows/columns; the boundary guards below must match their axis + # (matters when the two axes start at different offsets). + firstrow, firstcol = map(first, axes(A)) + lastrow, lastcol = map(last, axes(A)) + # `xs`/`ys` are 1-based but `A` may have offset axes: index `A` natively, the + # bounds positionally (`length(xs) == size(A, 1)` holds). + for (yj, j) in enumerate(axes(A, 2)) + y1, y2 = ys[yj] + for (xi, i) in enumerate(axes(A, 1)) if f(A[i, j]) # This is a pixel inside a polygon # xs and ys hold pixel bounds - x1, x2 = xs[i] + x1, x2 = xs[xi] # We check the Von Neumann neighborhood to # decide what edges are needed, if any. - (j == fi || !f(A[i, j-1])) && update_edge!(edges, (x1, y1), (x2, y1)) # S - (i == fj || !f(A[i-1, j])) && update_edge!(edges, (x1, y2), (x1, y1)) # W - (j == lj || !f(A[i, j+1])) && update_edge!(edges, (x2, y2), (x1, y2)) # N - (i == li || !f(A[i+1, j])) && update_edge!(edges, (x2, y1), (x2, y2)) # E + (j == firstcol || !f(A[i, j-1])) && update_edge!(edges, (x1, y1), (x2, y1)) # S + (i == firstrow || !f(A[i-1, j])) && update_edge!(edges, (x1, y2), (x1, y1)) # W + (j == lastcol || !f(A[i, j+1])) && update_edge!(edges, (x2, y2), (x1, y2)) # N + (i == lastrow || !f(A[i+1, j])) && update_edge!(edges, (x2, y1), (x2, y2)) # E end end end diff --git a/test/methods/polygonize.jl b/test/methods/polygonize.jl index b820c5c777..e6a6dfb957 100644 --- a/test/methods/polygonize.jl +++ b/test/methods/polygonize.jl @@ -1,9 +1,13 @@ using GeometryOps, GeoInterface, Test using GeometryOpsTestHelpers +using Random -import GeometryOps as GO +import GeometryOps as GO import GeoInterface as GI +# Seed so the coordinate-dependent assertions below (see #430) stay deterministic. +Random.seed!(123) + @testset "Polygonize with xs and ys, without offsetarrays" begin @test !(@isdefined OffsetArrays) # to make sure this isn't loaded somewhere else data = rand(1:4, 100, 100) .== 1 @@ -17,7 +21,8 @@ import GeoInterface as GI end # ideally we'd have a better test to make sure this returns what we think it does data_mp_range200 = polygonize(2:2:200, 2:2:200, data) - @test length(GI.coordinates(data_mp_range200)) == length(GI.coordinates(data_mp)) + # Broken: rescaling coordinates changes the decomposition (#430). + @test_broken length(GI.coordinates(data_mp_range200)) == length(GI.coordinates(data_mp)) # this is an example that could throw floating point error range_floats = -1.333333333333343:0.041666666666666664:0.374999999999986 @@ -33,7 +38,7 @@ import OffsetArrays, DimensionalData, Rasters for i in (100, 300), j in (100, 300) @testset "bool arrays without a function return MultiPolygon" begin A = rand(Bool, i, j) - @test_nowarn multipolygon = polygonize(A); + multipolygon = @test_nowarn polygonize(A); @test multipolygon isa GeoInterface.MultiPolygon @test GeoInterface.ngeom(multipolygon) > 0 end @@ -68,26 +73,43 @@ end @testset "Polygonize with exotic arrays" begin @testset "OffsetArrays" begin - data = rand(1:4, 100, 100) .== 1 + # Fixed pattern (blocks + a hole) so the result is deterministic. + data = fill(false, 12, 14) + data[2:6, 2:5] .= true + data[8:11, 7:12] .= true + data[3:5, 9:12] .= true + data[4, 10:11] .= false evil = OffsetArrays.Origin(-100, -100)(data) + evil_mp = @test_nowarn polygonize(evil) # regression: used to throw a `BoundsError` + # An offset array uses its axis values as coordinates. + @test GO.equals(evil_mp, polygonize(axes(evil, 1), axes(evil, 2), parent(evil))) + # Broken: not a translation of `polygonize(data)` (#430). data_mp = polygonize(data) - evil_mp = @test_nowarn polygonize(evil) - evil_in_data_space_mp = GO.transform(evil_mp) do point - point .- evil.offsets # undo the offset from the OffsetArray - end - @test GO.equals(data_mp, evil_in_data_space_mp) + shifted = GO.transform(p -> p .- evil.offsets, evil_mp) + @test_broken GO.equals(data_mp, shifted) + + # Regression: edge pixels + different per-axis offsets used to throw a `BoundsError`. + edgy = OffsetArrays.Origin(-100, -50)(fill(true, 5, 7)) + edgy_mp = @test_nowarn polygonize(edgy) + @test GO.equals(edgy_mp, polygonize(axes(edgy, 1), axes(edgy, 2), parent(edgy))) end + # Offset, non-square lookups: these only match `polygonize(xs, ys, data)` if the + # extension uses the lookup values rather than the integer axes. @testset "DimensionalData" begin - data = rand(1:4, 100, 100) .== 1 - evil = DimensionalData.DimArray(data, (DimensionalData.X(1:100), DimensionalData.Y(1:100))) - data_mp = polygonize(data) - evil_mp = @test_nowarn polygonize(evil) + data = rand(1:4, 100, 50) .== 1 + evil = DimensionalData.DimArray(data, (DimensionalData.X(51:150), DimensionalData.Y(151:200))) + data_mp = polygonize(51:150, 151:200, data) + evil_mp = @test_warn "Points" polygonize(evil) # `Points` sampling warns @test GO.equals(data_mp, evil_mp) end @testset "Rasters" begin - data = rand(1:4, 100, 100) .== 1 - evil = Rasters.Raster(data; dims = (DimensionalData.X(1:100), DimensionalData.Y(1:100)), crs = Rasters.GeoFormatTypes.EPSG(4326)) - data_mp = polygonize(data) + data = rand(1:4, 100, 50) .== 1 + # `Intervals` sampling (a raster's cells cover area): no warning. + evil = DimensionalData.set( + Rasters.Raster(data; dims = (DimensionalData.X(51:150), DimensionalData.Y(151:200)), crs = Rasters.GeoFormatTypes.EPSG(4326)), + DimensionalData.X => DimensionalData.Intervals(), DimensionalData.Y => DimensionalData.Intervals(), + ) + data_mp = polygonize(51:150, 151:200, data) evil_mp = @test_nowarn polygonize(evil) @test GO.equals(data_mp, evil_mp) @test GI.crs(evil_mp) == GI.crs(evil) diff --git a/test/runtests.jl b/test/runtests.jl index 935a4bd519..86c92d3002 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,6 +34,7 @@ end @safetestset "Distance" begin include("methods/distance.jl") end @safetestset "Equals" begin include("methods/equals.jl") end @safetestset "Minimum Bounding Circle" begin include("methods/minimum_bounding_circle.jl") end +@safetestset "Polygonize" begin include("methods/polygonize.jl") end # Clipping @safetestset "Coverage" begin include("methods/clipping/coverage.jl") end @safetestset "Cut" begin include("methods/clipping/cut.jl") end From f57277b08d6c82f44691ba6446688a9c24609660 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 10:10:02 -0700 Subject: [PATCH 07/20] Add `FlexibleRTrees` and speed up `NaturalIndex` traversal (#431) Co-authored-by: Claude Fable 5 --- src/GeometryOps.jl | 2 + src/utils/FlexibleRTrees/FlexibleRTrees.jl | 44 ++++++++++ src/utils/FlexibleRTrees/bulk_loading.jl | 70 ++++++++++++++++ src/utils/FlexibleRTrees/hilbert.jl | 95 ++++++++++++++++++++++ src/utils/FlexibleRTrees/interface.jl | 87 ++++++++++++++++++++ src/utils/FlexibleRTrees/types.jl | 81 ++++++++++++++++++ src/utils/NaturalIndexing.jl | 28 +++++-- test/runtests.jl | 1 + test/utils/FlexibleRTrees.jl | 86 ++++++++++++++++++++ 9 files changed, 485 insertions(+), 9 deletions(-) create mode 100644 src/utils/FlexibleRTrees/FlexibleRTrees.jl create mode 100644 src/utils/FlexibleRTrees/bulk_loading.jl create mode 100644 src/utils/FlexibleRTrees/hilbert.jl create mode 100644 src/utils/FlexibleRTrees/interface.jl create mode 100644 src/utils/FlexibleRTrees/types.jl create mode 100644 test/utils/FlexibleRTrees.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index ac2349b86f..58d380eb34 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -55,6 +55,8 @@ include("utils/utils.jl") include("utils/NaturalIndexing.jl") using .NaturalIndexing +include("utils/FlexibleRTrees/FlexibleRTrees.jl") +using .FlexibleRTrees # Load utility modules in using .NaturalIndexing, .SpatialTreeInterface, .LoopStateMachine diff --git a/src/utils/FlexibleRTrees/FlexibleRTrees.jl b/src/utils/FlexibleRTrees/FlexibleRTrees.jl new file mode 100644 index 0000000000..4ba94f558b --- /dev/null +++ b/src/utils/FlexibleRTrees/FlexibleRTrees.jl @@ -0,0 +1,44 @@ +# # FlexibleRTrees + +#= +A packed (bulk-loaded, static) R-tree over `Extents.Extent`s of any +dimensionality, with a pluggable bulk-load algorithm — sort-tile-recursive +([`STR`](@ref)), Hilbert-packed ([`HPR`](@ref)), or none +([`Unsorted`](@ref)) — behind one tree type. + +Storage is flat: `RTree{A, E}` holds a vector of per-level extent vectors +plus a leaf permutation, and is a concrete type at any size or depth. A +bulk-load algorithm chooses only the *leaf order*, via [`loadorder`](@ref); +packing always unions consecutive runs of `nodecapacity` extents, bottom-up. +Upper levels therefore group runs of the leaf order rather than re-tiling +each level: Hilbert order is spatially local at every scale so `HPR` packs +tightly, while `STR`'s upper levels are slightly looser than a re-tiled +pointer tree's. + +The tree implements SpatialTreeInterface, so `depth_first_search` / +`dual_depth_first_search` (and everything built on them) work unchanged. +Leaf queries yield indices into the *original* input collection. + +Parts of the construction logic are adapted from +[SortTileRecursiveTree.jl](https://github.com/maxfreu/SortTileRecursiveTree.jl) (MIT). + +```julia +tree = RTree(HPR(), extents) # or STR(), Unsorted() +hits = query(tree, Extents.Extent(X = (0, 1), Y = (0, 1))) +``` +=# + +module FlexibleRTrees + +import GeoInterface as GI +import Extents +using StaticArrays: MVector + +export RTree, BulkLoadAlgorithm, STR, HPR, Unsorted, query + +include("types.jl") # `BulkLoadAlgorithm`s and the `RTree` type +include("bulk_loading.jl") # `loadorder` methods and bottom-up packing +include("hilbert.jl") # Hilbert keys for `HPR`'s sort +include("interface.jl") # SpatialTreeInterface implementation and `query` + +end # module FlexibleRTrees diff --git a/src/utils/FlexibleRTrees/bulk_loading.jl b/src/utils/FlexibleRTrees/bulk_loading.jl new file mode 100644 index 0000000000..262762be07 --- /dev/null +++ b/src/utils/FlexibleRTrees/bulk_loading.jl @@ -0,0 +1,70 @@ +# # Bulk loading + +# ## Leaf ordering + +""" + loadorder(algorithm, extents::Vector{<:Extents.Extent}, nodecapacity)::Vector{Int} + +The permutation in which `algorithm` packs `extents` into leaves. Implement +this for a new `BulkLoadAlgorithm` subtype to plug in another ordering. +""" +loadorder(::Unsorted, extents, nodecapacity) = collect(1:length(extents)) +loadorder(::HPR, extents, nodecapacity) = sortperm(_hilbert_keys(extents)) +function loadorder(::STR, extents::Vector{E}, nodecapacity) where E + centers = [_center(e) for e in extents] + perm = collect(1:length(extents)) + _str_tile!(perm, centers, 1, length(extents), 1, _ndims(E), nodecapacity) + return perm +end + +#= +One recursion level of N-dimensional sort-tile-recursive: sort the range by +the current dimension's center, cut it into `S ≈ P^(1/remaining)` slabs of +whole leaf pages, and tile each slab along the remaining dimensions. After +the last dimension the consecutive `nodecapacity`-runs are the leaf tiles. +=# +function _str_tile!(perm, centers, lo, hi, dim, ndims, nodecapacity) + len = hi - lo + 1 + len <= nodecapacity && return # a single leaf: internal order doesn't matter + sort!(view(perm, lo:hi); by = i -> @inbounds(centers[i][dim])) + dim == ndims && return + P = cld(len, nodecapacity) # leaf pages in this range + S = ceil(Int, P^(1 / (ndims - dim + 1))) # slabs along this dimension + slab = cld(P, S) * nodecapacity # items per slab (whole pages) + i = lo + while i <= hi + _str_tile!(perm, centers, i, min(i + slab - 1, hi), dim + 1, ndims, nodecapacity) + i += slab + end + return +end + +# ## Packing + +# Union consecutive runs of `nodecapacity` extents, bottom-up, until a level +# fits in one (implicit) root node. Returns levels coarsest-first. +function _pack_levels(leaves::Vector{E}, nodecapacity::Int) where E + levels = [leaves] + current = leaves + while length(current) > nodecapacity + nparents = cld(length(current), nodecapacity) + parents = Vector{E}(undef, nparents) + for p in 1:nparents + lo = (p - 1) * nodecapacity + 1 + hi = min(p * nodecapacity, length(current)) + acc = current[lo] + for j in (lo + 1):hi + acc = Extents.union(acc, @inbounds current[j]) + end + parents[p] = acc + end + push!(levels, parents) + current = parents + end + return reverse!(levels) +end + +# ## Extent helpers + +_ndims(::Type{Extents.Extent{K, V}}) where {K, V} = length(K) +_center(ext::Extents.Extent) = map(b -> (b[1] + b[2]) / 2, values(ext)) diff --git a/src/utils/FlexibleRTrees/hilbert.jl b/src/utils/FlexibleRTrees/hilbert.jl new file mode 100644 index 0000000000..ac4af301b4 --- /dev/null +++ b/src/utils/FlexibleRTrees/hilbert.jl @@ -0,0 +1,95 @@ +# # Hilbert curve encoding +# +#= +The HPR ("Hilbert-packed R-tree") bulk loader sorts items by the position of +their center on a Hilbert space-filling curve. The Hilbert curve visits every +cell of a `2^bits × … × 2^bits` grid exactly once, and consecutive cells along +the curve are always spatially adjacent — so items that are consecutive in +sorted order are close in space, in every dimension, at every scale. That +fractal locality is what makes simple "pack consecutive runs" tree +construction produce tight nodes on all levels. + +This is Skilling's transpose-based algorithm ("Programming the Hilbert +curve", AIP Conf. Proc. 707, 2004), which works in any number of dimensions — +unlike the lookup-table encoders (e.g. JTS's 2D-only `HilbertCode`). The +"transpose" is the Hilbert index of the point stored bit-interleaved across +the N input coordinates; we finish by de-interleaving it into a single +integer sort key. +=# + +""" + hilbert_key(coords::NTuple{N, UInt32}, bits::Int) -> UInt64 + +The Hilbert-curve index of a point on the `N`-dimensional `2^bits` grid, as a +sortable integer. Requires `N * bits <= 64`; each coordinate must be +`< 2^bits`. +""" +function hilbert_key(coords::NTuple{N, UInt32}, bits::Int) where N + X = MVector{N, UInt32}(coords) + M = one(UInt32) << (bits - 1) + # Inverse undo + Q = M + while Q > one(UInt32) + P = Q - one(UInt32) + for i in 1:N + if !iszero(X[i] & Q) + X[1] ⊻= P # invert + else # exchange + t = (X[1] ⊻ X[i]) & P + X[1] ⊻= t + X[i] ⊻= t + end + end + Q >>= 1 + end + # Gray encode + for i in 2:N + X[i] ⊻= X[i - 1] + end + t = zero(UInt32) + Q = M + while Q > one(UInt32) + !iszero(X[N] & Q) && (t ⊻= Q - one(UInt32)) + Q >>= 1 + end + for i in 1:N + X[i] ⊻= t + end + # Interleave the transpose into one key, most significant bits first. + key = zero(UInt64) + for b in (bits - 1):-1:0 + for i in 1:N + key = (key << 1) | UInt64((X[i] >> b) & one(UInt32)) + end + end + return key +end + +# Bits per dimension for the quantization grid: 16 where possible (the JTS +# HPRtree uses 12), fewer in high dimensions so the key fits in 64 bits. +hilbert_bits(N::Int) = min(16, 63 ÷ N) + +#= +Quantize the centers of a set of extents onto the Hilbert grid spanned by +their total extent, and return each center's Hilbert key. Degenerate +dimensions (zero span) collapse to grid coordinate 0. +=# +function _hilbert_keys(extents::Vector{E}) where E <: Extents.Extent + N = _ndims(E) + bits = hilbert_bits(N) + total = reduce(Extents.union, extents) + los = map(first, values(total)) + spans = map(b -> Float64(b[2] - b[1]), values(total)) + scale = Float64((1 << bits) - 1) + keys = Vector{UInt64}(undef, length(extents)) + for (j, e) in enumerate(extents) + c = _center(e) + q = ntuple(Val(N)) do i + span = spans[i] + frac = iszero(span) ? 0.0 : (Float64(c[i]) - Float64(los[i])) / span + UInt32(clamp(round(Int, frac * scale), 0, Int(scale))) + end + keys[j] = hilbert_key(q, bits) + end + return keys +end diff --git a/src/utils/FlexibleRTrees/interface.jl b/src/utils/FlexibleRTrees/interface.jl new file mode 100644 index 0000000000..3e4c140f88 --- /dev/null +++ b/src/utils/FlexibleRTrees/interface.jl @@ -0,0 +1,87 @@ +# # SpatialTreeInterface + +using ..SpatialTreeInterface +import ..SpatialTreeInterface: isspatialtree, isleaf, nchild, getchild, + child_indices_extents, depth_first_search + +""" + RTreeNode{A, E} + +A cursor into one node of an [`RTree`](@ref): the tree, the node's level +(0-based; the children of a level-`l` node live in `levels[l + 1]`), its +position within that level, and its extent. All SpatialTreeInterface +methods traverse the tree through these cursors. The children of one node +occupy one contiguous run of the next level's extent vector, which each +per-child method resolves once per node and then indexes into. At the leaf +level, `child_indices_extents` maps leaf slots through `tree.indices`, so +queries return indices into the original collection despite the packed +reordering. +""" +struct RTreeNode{A <: BulkLoadAlgorithm, E <: Extents.Extent} + tree::RTree{A, E} + level::Int # 0-based; children of a level-l node live in levels[l + 1] + index::Int # position within its level + extent::E +end + +Extents.extent(node::RTreeNode) = node.extent + +isspatialtree(::Type{<:RTree}) = true +isspatialtree(::Type{<:RTreeNode}) = true + +@inline _child_extents(node::RTreeNode) = node.tree.levels[node.level + 1] + +@inline function _child_range(node::RTreeNode, child_extents) + start_idx = (node.index - 1) * node.tree.nodecapacity + 1 + stop_idx = min(start_idx + node.tree.nodecapacity - 1, length(child_extents)) + return start_idx:stop_idx +end + +isleaf(node::RTreeNode) = node.level == length(node.tree.levels) - 1 + +nchild(node::RTreeNode) = length(_child_range(node, _child_extents(node))) + +function getchild(node::RTreeNode, i::Int) + child_index = (node.index - 1) * node.tree.nodecapacity + i + return RTreeNode(node.tree, node.level + 1, child_index, _child_extents(node)[child_index]) +end + +function getchild(node::RTreeNode) + extents = _child_extents(node) + tree, childlevel = node.tree, node.level + 1 + range = _child_range(node, extents) + return (RTreeNode(tree, childlevel, ci, @inbounds extents[ci]) for ci in range) +end + +function child_indices_extents(node::RTreeNode) + extents = _child_extents(node) + indices = node.tree.indices + range = _child_range(node, extents) + return ((@inbounds(indices[i]), @inbounds(extents[i])) for i in range) +end + +# The tree itself acts as the (implicit) root node. +_rootnode(tree::RTree) = RTreeNode(tree, 0, 1, tree.extent) + +isleaf(tree::RTree) = length(tree.levels) == 1 +nchild(tree::RTree) = length(tree.levels[1]) +getchild(tree::RTree) = getchild(_rootnode(tree)) +getchild(tree::RTree, i) = getchild(_rootnode(tree), i) +child_indices_extents(tree::RTree) = child_indices_extents(_rootnode(tree)) + +# ## Queries + +""" + query(tree::RTree, extent_or_geom) + +Indices (into the collection the tree was built from) of every leaf whose +extent intersects the given extent — or the extent of the given geometry — +in ascending order. +""" +query(tree::RTree, ext::Extents.Extent) = + sort!(depth_first_search(Base.Fix1(Extents.intersects, ext), tree)) +function query(tree::RTree, geom) + ext = GI.extent(geom) + isnothing(ext) && throw(ArgumentError("no extent found on $(typeof(geom))")) + return query(tree, ext) +end diff --git a/src/utils/FlexibleRTrees/types.jl b/src/utils/FlexibleRTrees/types.jl new file mode 100644 index 0000000000..dfb2151059 --- /dev/null +++ b/src/utils/FlexibleRTrees/types.jl @@ -0,0 +1,81 @@ +# # Types + +# ## Bulk-load algorithms + +""" + BulkLoadAlgorithm + +Supertype for the algorithms that decide the *leaf order* of an [`RTree`](@ref). +Packing is always "union consecutive runs of `nodecapacity`, bottom-up"; the +algorithm only chooses the order, via a [`loadorder`](@ref) method. +""" +abstract type BulkLoadAlgorithm end + +""" + STR() + +Sort-tile-recursive ordering (Leutenegger et al., 1997), generalized to any +dimensionality: sort by center along the first dimension, cut into slabs, +recurse within each slab on the remaining dimensions. +""" +struct STR <: BulkLoadAlgorithm end + +""" + HPR() + +Hilbert-packed ordering, as in JTS's `HPRtree`: sort by the Hilbert-curve +index of each extent's center. Hilbert order is spatially local at every +scale, which suits this tree's consecutive-run packing particularly well. +""" +struct HPR <: BulkLoadAlgorithm end + +""" + Unsorted() + +Keep the input order (no sort). Equivalent to natural indexing — good when +the input is already spatially coherent (e.g. the edges of a ring), and the +baseline the sorting algorithms have to beat. +""" +struct Unsorted <: BulkLoadAlgorithm end + +# ## The tree + +""" + RTree(algorithm::BulkLoadAlgorithm, data; nodecapacity = 16) + +A packed R-tree over the extents of `data` (anything `GI.extent` accepts — +geometries, or `Extents.Extent`s themselves), of any dimensionality, bulk +loaded in the order chosen by `algorithm`. + +The tree is flat and fully concrete: `levels[1]` is the coarsest level and +`levels[end]` holds the leaf extents in packed order, with `indices` mapping +each leaf slot back to its position in `data`. Queries through +SpatialTreeInterface therefore return indices into the original collection. +""" +struct RTree{A <: BulkLoadAlgorithm, E <: Extents.Extent} + algorithm::A + nodecapacity::Int + extent::E + levels::Vector{Vector{E}} # levels[1] = coarsest, levels[end] = leaf extents (packed order) + indices::Vector{Int} # leaf slot -> index into the original collection +end + +function RTree(algorithm::A, data; nodecapacity::Int = 16) where A <: BulkLoadAlgorithm + nodecapacity >= 2 || throw(ArgumentError("`nodecapacity` must be at least 2, got $nodecapacity")) + isnothing(iterate(data)) && throw(ArgumentError("cannot build an `RTree` from an empty collection")) + E = typeof(GI.extent(first(data))) + extents = E[GI.extent(x) for x in data] + perm = loadorder(algorithm, extents, nodecapacity) + leaves = extents[perm] + levels = _pack_levels(leaves, nodecapacity) + total = reduce(Extents.union, levels[1]) + return RTree{A, E}(algorithm, nodecapacity, total, levels, perm) +end + +Extents.extent(tree::RTree) = tree.extent + +function Base.show(io::IO, tree::RTree{A}) where A + print(io, "RTree{", nameof(A), "}(", length(tree.indices), " leaves, ", + length(tree.levels), " levels, capacity ", tree.nodecapacity, ")") +end +Base.show(io::IO, ::MIME"text/plain", tree::RTree) = Base.show(io, tree) diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl index 080da277e4..03fb069ed6 100644 --- a/src/utils/NaturalIndexing.jl +++ b/src/utils/NaturalIndexing.jl @@ -155,33 +155,43 @@ Extents.extent(node::NaturalIndexNode) = node.extent SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndex}) = true SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndexNode}) = true -function SpatialTreeInterface.nchild(node::NaturalIndexNode) +# A node's children all live in one extents vector at the next level. +# Resolve that vector once per node and index into it, instead of following +# `parent_index -> levels -> level -> extents` again for every child. +@inline _child_extents(node::NaturalIndexNode) = node.parent_index.levels[node.level + 1].extents + +@inline function _child_range(node::NaturalIndexNode, child_extents) start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1 - stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents)) - return stop_idx - start_idx + 1 + stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(child_extents)) + return start_idx:stop_idx end +SpatialTreeInterface.nchild(node::NaturalIndexNode) = length(_child_range(node, _child_extents(node))) + function SpatialTreeInterface.getchild(node::NaturalIndexNode, i::Int) child_index = (node.index - 1) * node.parent_index.nodecapacity + i return NaturalIndexNode( - node.parent_index, + node.parent_index, node.level + 1, # increment level by 1 child_index, # index of this particular child - node.parent_index.levels[node.level+1].extents[child_index] # the extent of this child + _child_extents(node)[child_index] # the extent of this child ) end # Get all children of a node function SpatialTreeInterface.getchild(node::NaturalIndexNode) - return (SpatialTreeInterface.getchild(node, i) for i in 1:SpatialTreeInterface.nchild(node)) + extents = _child_extents(node) + parent, childlevel = node.parent_index, node.level + 1 + range = _child_range(node, extents) + return (NaturalIndexNode(parent, childlevel, ci, @inbounds extents[ci]) for ci in range) end SpatialTreeInterface.isleaf(node::NaturalIndexNode) = node.level == length(node.parent_index.levels) - 1 function SpatialTreeInterface.child_indices_extents(node::NaturalIndexNode) - start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1 - stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents)) - return ((i, node.parent_index.levels[node.level+1].extents[i]) for i in start_idx:stop_idx) + extents = _child_extents(node) + range = _child_range(node, extents) + return ((i, @inbounds extents[i]) for i in range) end # implementation for "root node" / top level tree diff --git a/test/runtests.jl b/test/runtests.jl index 86c92d3002..f67aceee5f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,7 @@ end @safetestset "Utils" begin include("utils/utils.jl") end @safetestset "LoopStateMachine" begin include("utils/LoopStateMachine.jl") end @safetestset "SpatialTreeInterface" begin include("utils/SpatialTreeInterface.jl") end +@safetestset "FlexibleRTrees" begin include("utils/FlexibleRTrees.jl") end @safetestset "UnitSpherical" begin include("utils/unitspherical.jl") end @safetestset "RobustCrossProduct" begin include("utils/robustcrossproduct.jl") end # Methods diff --git a/test/utils/FlexibleRTrees.jl b/test/utils/FlexibleRTrees.jl new file mode 100644 index 0000000000..58a57d17c9 --- /dev/null +++ b/test/utils/FlexibleRTrees.jl @@ -0,0 +1,86 @@ +using Test + +import GeoInterface as GI +import GeometryOps as GO +import GeometryOps.FlexibleRTrees as FRT +import GeometryOps.FlexibleRTrees: RTree, STR, HPR, Unsorted, query, hilbert_key +import GeometryOps.SpatialTreeInterface as STI +import Extents +using Random: Xoshiro + +# Random boxes with side lengths ~5% of the unit cube, in N dims. +function random_extents(rng, n, N) + dims = (:X, :Y, :Z, :M)[1:N] + return [begin + lo = ntuple(_ -> rand(rng), N) + hi = lo .+ 0.05 .* ntuple(_ -> rand(rng), N) + Extents.Extent(NamedTuple{dims}(tuple.(lo, hi))) + end for _ in 1:n] +end + +brute_force(ext, extents) = findall(e -> Extents.intersects(ext, e), extents) + +grow(ext, d) = Extents.buffer(ext, NamedTuple{keys(ext)}(ntuple(_ -> d, length(keys(ext))))) + +@testset "query ≡ brute force ($(N)D, $alg, n = $n)" for + N in (2, 3), + alg in (STR(), HPR(), Unsorted()), + n in (1, 5, 16, 17, 100, 256, 1000) + rng = Xoshiro(hash((N, n))) + extents = random_extents(rng, n, N) + tree = RTree(alg, extents; nodecapacity = 8) + queries = vcat( + random_extents(rng, 20, N), # small probes + [reduce(Extents.union, extents)], # everything + [grow(reduce(Extents.union, extents), 10.0)], # superset + [grow(e, 3.0) for e in random_extents(Xoshiro(0), 5, N)], # big probes + ) + for q in queries + @test query(tree, q) == brute_force(q, extents) + end + # A query far outside everything returns nothing. + faraway = Extents.Extent(NamedTuple{((:X, :Y, :Z, :M)[1:N])}(ntuple(_ -> (99.0, 100.0), N))) + @test isempty(query(tree, faraway)) +end + +@testset "construction and type stability" begin + rng = Xoshiro(7) + extents = random_extents(rng, 300, 2) + tree = @inferred RTree(STR(), extents) + @test tree isa RTree{STR, eltype(extents)} + @inferred RTree(HPR(), extents) + @inferred RTree(Unsorted(), extents) + # The query path is inferrable too (depth_first_search returns Vector{Int}). + q = Extents.Extent(X = (0.2, 0.4), Y = (0.2, 0.4)) + @inferred query(tree, q) + # Deep and shallow trees have the SAME concrete type — the point of the flat layout. + tiny = RTree(STR(), extents[1:3]) + @test typeof(tiny) === typeof(tree) + + # Geometries as input work through GI.extent. + lines = GO.to_edgelist(GI.LinearRing([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 0.0)])) + gtree = RTree(HPR(), lines) + # Edge 2 is the right side; edge 3 (the closing diagonal) has a bbox + # covering the whole square, so extent-intersects finds it too. + @test query(gtree, Extents.Extent(X = (0.9, 1.1), Y = (0.4, 0.6))) == [2, 3] + + @test_throws ArgumentError RTree(STR(), Extents.Extent{(:X, :Y)}[]) + @test_throws ArgumentError RTree(STR(), random_extents(rng, 5, 2); nodecapacity = 1) + @test occursin("RTree{HPR}", sprint(show, gtree)) +end + +@testset "Hilbert curve properties" begin + # Order-1 2D curve: the classic U through the four quadrants. + keys1 = [hilbert_key((UInt32(x), UInt32(y)), 1) for (x, y) in ((0, 0), (0, 1), (1, 1), (1, 0))] + @test keys1 == [0, 1, 2, 3] + # In any dimension: the curve visits every grid cell exactly once + # (bijectivity), and consecutive cells are adjacent (unit Manhattan step). + for (N, bits) in ((2, 4), (3, 2)) + side = 2^bits + cells = vec(collect(Iterators.product(ntuple(_ -> 0:(side - 1), N)...))) + keys = [hilbert_key(UInt32.(c), bits) for c in cells] + @test allunique(keys) + path = cells[sortperm(keys)] + @test all(sum(abs.(path[i + 1] .- path[i])) == 1 for i in 1:(length(path) - 1)) + end +end From d2084f4963f001f2974953b1a4b8082c9d4c5741 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 15:18:40 -0700 Subject: [PATCH 08/20] Bump version from 0.1.41 to 0.1.42 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index fde30baf29..90ae4f8c41 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeometryOps" uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" authors = ["Anshul Singhvi ", "Rafael Schouten ", "Skylar Gering ", "and contributors"] -version = "0.1.41" +version = "0.1.42" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" From 74d634a13aa9802ffc7da1f85d687fcca4a80782 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 16:06:32 -0700 Subject: [PATCH 09/20] Allow spherical caps to intersect 3D extents (#432) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Allow `SphericalCap`s to intersect 3D `Extents.Extent`s Compare the squared Euclidean distance from the cap center to the box against the squared chord radius — the same comparison as S2's `S2Cap::Contains(S2Point)`, extended to a box by clamping the center into it. The chord radius is `2sin(min(radius, π)/2)` following S2's `S1ChordAngle`, which keeps precision at tiny radii where `2 - 2cos(radius)` cancels to zero. Conservative for spatial tree pruning: never a false negative for geometry on the sphere. Co-Authored-By: Claude Fable 5 * improve docstring Co-authored-by: Anshul Singhvi --------- Co-authored-by: Claude Fable 5 --- src/utils/UnitSpherical/UnitSpherical.jl | 1 + src/utils/UnitSpherical/cap.jl | 28 +++++++++++++++ test/utils/unitspherical.jl | 44 +++++++++++++++++++++++- 3 files changed, 72 insertions(+), 1 deletion(-) diff --git a/src/utils/UnitSpherical/UnitSpherical.jl b/src/utils/UnitSpherical/UnitSpherical.jl index c6293b64c9..4b4a535164 100644 --- a/src/utils/UnitSpherical/UnitSpherical.jl +++ b/src/utils/UnitSpherical/UnitSpherical.jl @@ -3,6 +3,7 @@ module UnitSpherical using CoordinateTransformations using StaticArrays, LinearAlgebra import GeoInterface as GI, GeoFormatTypes as GFT +import Extents import Random diff --git a/src/utils/UnitSpherical/cap.jl b/src/utils/UnitSpherical/cap.jl index 483ddebdbf..8956a583e0 100644 --- a/src/utils/UnitSpherical/cap.jl +++ b/src/utils/UnitSpherical/cap.jl @@ -113,6 +113,34 @@ function _contains(cap::SphericalCap, point::UnitSphericalPoint) spherical_distance(cap.point, point) <= cap.radius end +# ## Cap–extent intersection + +""" + Extents.intersects(cap::SphericalCap, ext::Extents.Extent{(:X, :Y, :Z)}) + Extents.intersects(ext::Extents.Extent{(:X, :Y, :Z)}, cap::SphericalCap) + +Whether `cap` intersects the part of the unit sphere covered by the 3D +Cartesian bounding box `ext` (also in unit-spherical space). + +A point on the sphere lies in the cap iff its Euclidean (chord) distance to +the cap's center is at most ``2 * sin(radius/2)``, so this tests whether the box +comes within that distance of the center. + +This is a fail-safe comparison, so may have false positives but never false negatives. +""" +function Extents.intersects(cap::SphericalCap, ext::Extents.Extent{(:X, :Y, :Z)}) + c = cap.point + dx = c.x - clamp(c.x, ext.X[1], ext.X[2]) + dy = c.y - clamp(c.y, ext.Y[1], ext.Y[2]) + dz = c.z - clamp(c.z, ext.Z[1], ext.Z[2]) + # the chord radius as S2's `S1ChordAngle(S1Angle)` computes it: accurate + # for small radii, where `2 - 2 * radiuslike` cancels to 0 + chord = 2 * sin(0.5 * min(cap.radius, π)) + return dx^2 + dy^2 + dz^2 <= chord^2 +end +Extents.intersects(ext::Extents.Extent{(:X, :Y, :Z)}, cap::SphericalCap) = + Extents.intersects(cap, ext) + #Comment by asinghvi: this could be transformed to GO.union function _merge(x::SphericalCap, y::SphericalCap) diff --git a/test/utils/unitspherical.jl b/test/utils/unitspherical.jl index 1bba6b8815..9a0fe41f15 100644 --- a/test/utils/unitspherical.jl +++ b/test/utils/unitspherical.jl @@ -3,7 +3,9 @@ using LinearAlgebra using GeometryOps.UnitSpherical +import GeometryOps as GO import GeoInterface as GI +import Extents @testset "spherical_distance" begin @testset "Basic correctness" begin @@ -544,4 +546,44 @@ end result = spherical_arc_intersection(a7, b7, a8, b8) @test result.type == arc_overlap @test length(result.points) == 2 -end \ No newline at end of file +end + +@testset "Cap–extent intersection" begin + cap = SphericalCap(UnitSphericalPoint(1.0, 0.0, 0.0), π/4) + around_center = Extents.Extent(X = (0.9, 1.1), Y = (-0.1, 0.1), Z = (-0.1, 0.1)) + antipodal = Extents.Extent(X = (-1.1, -0.9), Y = (-0.1, 0.1), Z = (-0.1, 0.1)) + @test Extents.intersects(cap, around_center) + @test Extents.intersects(around_center, cap) + @test !Extents.intersects(cap, antipodal) + @test !Extents.intersects(antipodal, cap) + + # Degenerate (point) extents straddling the cap boundary + just_inside = UnitSphericalPoint(cos(π/4 - 1e-3), sin(π/4 - 1e-3), 0.0) + just_outside = UnitSphericalPoint(cos(π/4 + 1e-3), sin(π/4 + 1e-3), 0.0) + @test Extents.intersects(cap, GI.extent(just_inside)) + @test !Extents.intersects(cap, GI.extent(just_outside)) + + # A whole-sphere cap reaches the antipodal box; radii past π clamp to full + @test Extents.intersects(SphericalCap(UnitSphericalPoint(1.0, 0.0, 0.0), π), antipodal) + @test Extents.intersects(SphericalCap(UnitSphericalPoint(1.0, 0.0, 0.0), 3π/2), antipodal) + + # Tiny radii keep full precision (the chord radius is computed as + # `2sin(radius/2)`, not from `radiuslike = cos(radius)`, which rounds + # to 1 below radius ≈ 1.5e-8 and would empty the cap) + tinycap = SphericalCap(UnitSphericalPoint(1.0, 0.0, 0.0), 1e-9) + @test Extents.intersects(tinycap, GI.extent(UnitSphericalPoint(cos(0.9e-9), sin(0.9e-9), 0.0))) + @test !Extents.intersects(tinycap, GI.extent(UnitSphericalPoint(cos(1.1e-9), sin(1.1e-9), 0.0))) + + @testset "Spatial tree pruning against brute force" begin + using Random: Xoshiro + rng = Xoshiro(1312) + points = rand(rng, UnitSphericalPoint{Float64}, 1000) + tree = GO.FlexibleRTrees.RTree(GO.FlexibleRTrees.HPR(), points) + for radius in (0.05, 0.5, 2.0, π) + querycap = SphericalCap(rand(rng, UnitSphericalPoint{Float64}), radius) + hits = sort!(GO.SpatialTreeInterface.depth_first_search( + Base.Fix1(Extents.intersects, querycap), tree)) + @test hits == findall(p -> spherical_distance(querycap.point, p) <= querycap.radius, points) + end + end +end From 7145c7cb53457981451b7f287b01a5c09ccffa78 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 16:27:21 -0700 Subject: [PATCH 10/20] Add `spherical_arc_extent` and spherical manifold edge decomposition MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A great-circle arc bulges out of its endpoints' bounding box, so edge extents built from endpoints alone silently miss queries that touch only the bulge. `spherical_arc_extent(a, b)` computes the true 3D Cartesian extent of the shorter arc: each coordinate along the arc is a sinusoid, so an interior extremum exists iff the coordinate rises at one endpoint and falls at the other, where it attains the sinusoid's amplitude — no trigonometric calls, with tangents from `robust_cross_product` so nearly-degenerate arcs stay stable. Bounds are padded by a few ulps to guarantee containment, as S2's `S2LatLngRectBounder` does. `eachedge`, `to_edgelist`, and `lazy_edgelist` gain manifold-first methods: `Spherical()` yields edges as `UnitSphericalPoint` pairs (converting geographic input) whose `GI.Line`s carry arc extents, and rings that are already `UnitSphericalPoint`s take that path with no manifold argument. Previously `_lineedge` threw a `MethodError` for such rings. Spatial trees built over these edges (`NaturalIndex`, `RTree`) therefore index spherical edges in 3D, correctly. Co-Authored-By: Claude Fable 5 --- src/utils/UnitSpherical/UnitSpherical.jl | 2 + src/utils/UnitSpherical/arc_extent.jl | 79 ++++++++++++++++++++++++ src/utils/utils.jl | 30 ++++++++- test/utils/unitspherical.jl | 47 ++++++++++++++ test/utils/utils.jl | 37 +++++++++++ 5 files changed, 194 insertions(+), 1 deletion(-) create mode 100644 src/utils/UnitSpherical/arc_extent.jl diff --git a/src/utils/UnitSpherical/UnitSpherical.jl b/src/utils/UnitSpherical/UnitSpherical.jl index 4b4a535164..cdc70b62e5 100644 --- a/src/utils/UnitSpherical/UnitSpherical.jl +++ b/src/utils/UnitSpherical/UnitSpherical.jl @@ -21,11 +21,13 @@ include("slerp.jl") include("cap.jl") include("predicates.jl") include("arc_intersection.jl") +include("arc_extent.jl") export UnitSphericalPoint, UnitSphereFromGeographic, GeographicFromUnitSphere, slerp, SphericalCap, spherical_distance, spherical_orient, point_on_spherical_arc, spherical_arc_intersection, ArcIntersectionResult, arc_cross, arc_hinge, arc_overlap, arc_disjoint, + spherical_arc_extent, to_unit_spherical_points """ diff --git a/src/utils/UnitSpherical/arc_extent.jl b/src/utils/UnitSpherical/arc_extent.jl new file mode 100644 index 0000000000..c571b58c85 --- /dev/null +++ b/src/utils/UnitSpherical/arc_extent.jl @@ -0,0 +1,79 @@ +# # Spherical arc extents + +#= +```@docs; canonical=false +spherical_arc_extent +``` + +## Why not the extent of the endpoints? + +A great-circle arc bulges away from the chord between its endpoints, so the +axis-aligned bounding box of the endpoints does not, in general, contain the +arc. The classic case is two points at the same latitude: the arc between +them passes closer to the pole than either endpoint, e.g. two points at +`z = 0.9` on either side of the prime meridian are joined by an arc whose +midpoint has `z = 0.9 / cos(θ/2) > 0.9`. A spatial index built on endpoint +boxes would silently miss queries that touch only the bulge. + +## How the extremum is found + +With `t̂ₐ` the unit tangent at `a` pointing along the arc, the arc is +`p(φ) = a cos(φ) + t̂ₐ sin(φ)` for `φ ∈ [0, θ]`, so each Cartesian +coordinate is a sinusoid `pᵢ(φ) = Rᵢ cos(φ - φᵢ)` with amplitude +`Rᵢ = hypot(aᵢ, t̂ₐᵢ)`. Since `θ ≤ π`, at most one interior maximum and one +interior minimum exist per axis, and the endpoint derivatives decide: an +interior maximum exists iff `pᵢ` is increasing at `a` and decreasing at `b` +(`t̂ₐᵢ > 0 > t̂ᵦᵢ`), where it attains `Rᵢ`; likewise a minimum attains `-Rᵢ`. +No trigonometric calls are needed, and the tangents come from +[`robust_cross_product`](@ref), so nearly-degenerate and nearly-antipodal +arcs stay stable. + +Bounds are padded by a few ulps so that the extent is guaranteed to contain +the arc despite floating point error, in the same spirit as S2's +`S2LatLngRectBounder`, which widens its bounds by their maximum error. +=# + +""" + spherical_arc_extent(a, b)::Extents.Extent{(:X, :Y, :Z)} + +The 3D Cartesian extent of the shorter great-circle arc between `a` and `b` +on the unit sphere. Accepts `UnitSphericalPoint`s, or any GeoInterface +point (interpreted geographically, as longitude/latitude, like the +`UnitSphericalPoint` constructor itself). + +The extent is exact up to floating point error and padded by a few ulps, so +it always contains the arc — unlike the extent of the endpoints, which the +arc bulges out of wherever a coordinate attains its extremum between them. +For antipodal endpoints the arc's plane is ambiguous; the one chosen by +[`robust_cross_product`](@ref) is used. + +## Example + +```jldoctest +julia> using GeometryOps.UnitSpherical + +julia> ext = spherical_arc_extent(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(0, 1, 0)); + +julia> ext.X[2] ≈ 1 && ext.Y[2] ≈ 1 +true +``` +""" +spherical_arc_extent(a, b) = spherical_arc_extent(UnitSphericalPoint(a), UnitSphericalPoint(b)) +function spherical_arc_extent(a::UnitSphericalPoint{T1}, b::UnitSphericalPoint{T2}) where {T1, T2} + F = float(promote_type(T1, T2)) + pad = 4 * eps(F) + if a == b + bounds = ntuple(i -> (F(a[i]) - pad, F(a[i]) + pad), 3) + else + n = robust_cross_product(a, b) + ta = normalize(cross(n, a)) # unit tangent at `a`, pointing along the arc + tb = normalize(cross(n, b)) # unit tangent at `b`, pointing along the arc + bounds = ntuple(3) do i + lo, hi = minmax(F(a[i]), F(b[i])) + ta[i] > 0 > tb[i] && (hi = hypot(F(a[i]), F(ta[i]))) + ta[i] < 0 < tb[i] && (lo = -hypot(F(a[i]), F(ta[i]))) + (lo - pad, hi + pad) + end + end + return Extents.Extent(X = bounds[1], Y = bounds[2], Z = bounds[3]) +end diff --git a/src/utils/utils.jl b/src/utils/utils.jl index ddefac1a48..7aae0831b2 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -159,6 +159,7 @@ Currently they only work on linear rings. """ eachedge(geom, [::Type{T}]) + eachedge(m::Manifold, geom, [::Type{T}]) Decompose a geometry into a list of edges. Currently only works for LineString and LinearRing. @@ -166,11 +167,19 @@ Currently only works for LineString and LinearRing. Returns some iterator, which yields tuples of points. Each tuple is an edge. It goes `(p1, p2), (p2, p3), (p3, p4), ...` etc. + +On `Planar()` (the manifold-less default) points are 2D coordinate tuples. +On `Spherical()` they are `UnitSphericalPoint`s: geographic (longitude, +latitude) input is converted, and `UnitSphericalPoint`s pass through as-is. """ eachedge(geom) = eachedge(GI.trait(geom), geom, Float64) function eachedge(geom, ::Type{T}) where T eachedge(GI.trait(geom), geom, T) end +eachedge(::Planar, geom, ::Type{T} = Float64) where T = eachedge(geom, T) +function eachedge(::Spherical, geom, ::Type{T} = Float64) where T + return (map(UnitSpherical.UnitSphericalPoint, ps) for ps in eachedge(geom, T)) +end # implementation for LineString and LinearRing function eachedge(trait::GI.AbstractCurveTrait, geom, ::Type{T}) where T return (_tuple_point.((GI.getpoint(geom, i), GI.getpoint(geom, i+1)), T) for i in 1:GI.npoint(geom)-1) @@ -188,11 +197,20 @@ end """ to_edgelist(geom, [::Type{T}]) + to_edgelist(m::Manifold, geom, [::Type{T}]) Convert a geometry into a vector of `GI.Line` objects with attached extents. + +On `Spherical()` — or whenever the geometry's points are already +`UnitSphericalPoint`s — each edge is a great-circle arc of +`UnitSphericalPoint`s carrying its 3D [`UnitSpherical.spherical_arc_extent`](@ref), +so spatial indices built over the edges (e.g. `NaturalIndex`, `RTree`) +bound the arcs correctly. """ -to_edgelist(geom, ::Type{T} = Float64) where T = +to_edgelist(geom, ::Type{T} = Float64) where T = [_lineedge(ps, T) for ps in eachedge(geom, T)] +to_edgelist(m::Manifold, geom, ::Type{T} = Float64) where T = + [_lineedge(ps, T) for ps in eachedge(m, geom, T)] """ to_edgelist(ext::E, geom, [::Type{T}])::(::Vector{GI.Line}, ::Vector{Int}) @@ -222,15 +240,25 @@ function _lineedge(ps::Tuple, ::Type{T}) where T e = GI.extent(l) return GI.Line(l.geom; extent=e) end +# On the unit sphere, an edge is a great-circle arc: its extent is 3D and +# must cover the arc's bulge past the endpoints, not just the endpoints. +function _lineedge(ps::Tuple{<:UnitSpherical.UnitSphericalPoint, <:UnitSpherical.UnitSphericalPoint}, ::Type{T}) where T + a, b = UnitSpherical.UnitSphericalPoint{T}.(ps) + return GI.Line(StaticArrays.SVector((a, b)); extent = UnitSpherical.spherical_arc_extent(a, b)) +end """ lazy_edgelist(geom, [::Type{T}]) + lazy_edgelist(m::Manifold, geom, [::Type{T}]) Return an iterator over `GI.Line` objects with attached extents. """ function lazy_edgelist(geom, ::Type{T} = Float64) where T (_lineedge(ps, T) for ps in eachedge(geom, T)) end +function lazy_edgelist(m::Manifold, geom, ::Type{T} = Float64) where T + (_lineedge(ps, T) for ps in eachedge(m, geom, T)) +end """ edge_extents(geom, [::Type{T}]) diff --git a/test/utils/unitspherical.jl b/test/utils/unitspherical.jl index 9a0fe41f15..a390cf957a 100644 --- a/test/utils/unitspherical.jl +++ b/test/utils/unitspherical.jl @@ -587,3 +587,50 @@ end end end end + +@testset "spherical_arc_extent" begin + # Quarter arc along the equator: x and y peak at the endpoints, z is flat + ext = spherical_arc_extent(UnitSphericalPoint(1.0, 0.0, 0.0), UnitSphericalPoint(0.0, 1.0, 0.0)) + @test ext isa Extents.Extent{(:X, :Y, :Z)} + @test ext.X[1] ≈ 0 atol = 1e-14 + @test ext.X[2] ≈ 1 atol = 1e-14 + @test ext.Y[1] ≈ 0 atol = 1e-14 + @test ext.Y[2] ≈ 1 atol = 1e-14 + @test ext.Z[1] ≈ 0 atol = 1e-14 + @test ext.Z[2] ≈ 0 atol = 1e-14 + + # The arc bulges out of the endpoints' box: two points at z = 0.9 on + # either side of the x = 0 great circle, joined over the pole + z = 0.9; s = sqrt(1 - z^2) + bulging = spherical_arc_extent(UnitSphericalPoint(0.0, -s, z), UnitSphericalPoint(0.0, s, z)) + @test bulging.Z[2] ≈ 1 atol = 1e-14 # the pole, not the endpoints' z = 0.9 + @test bulging.Z[2] > z + @test bulging.Z[1] ≈ z atol = 1e-14 + @test bulging.Y[1] ≈ -s atol = 1e-14 + @test bulging.Y[2] ≈ s atol = 1e-14 + + # Geographic (lon, lat) input takes the same path + @test spherical_arc_extent((0.0, 0.0), (90.0, 0.0)) == + spherical_arc_extent(UnitSphericalPoint(1.0, 0.0, 0.0), UnitSphericalPoint(0.0, 1.0, 0.0)) + + # Degenerate arc: a point + p = UnitSphericalPoint(1.0, 0.0, 0.0) + degenerate = spherical_arc_extent(p, p) + @test degenerate.X[1] <= 1 <= degenerate.X[2] + @test degenerate.X[2] - degenerate.X[1] < 1e-14 + + @testset "Containment of dense arc samples" begin + using Random: Xoshiro + rng = Xoshiro(999) + inext(q, e) = e.X[1] <= q[1] <= e.X[2] && e.Y[1] <= q[2] <= e.Y[2] && e.Z[1] <= q[3] <= e.Z[2] + for _ in 1:100 + a, b = rand(rng, UnitSphericalPoint{Float64}), rand(rng, UnitSphericalPoint{Float64}) + e = spherical_arc_extent(a, b) + @test all(t -> inext(slerp(a, b, t), e), range(0.0, 1.0, length = 101)) + # nearly-degenerate arcs stay stable through robust_cross_product + b2 = UnitSphericalPoint(normalize(a + 1e-12 * rand(rng, UnitSphericalPoint{Float64}))) + e2 = spherical_arc_extent(a, b2) + @test all(t -> inext(slerp(a, b2, t), e2), range(0.0, 1.0, length = 11)) + end + end +end diff --git a/test/utils/utils.jl b/test/utils/utils.jl index 79ddf4e301..de1648dc18 100644 --- a/test/utils/utils.jl +++ b/test/utils/utils.jl @@ -2,6 +2,7 @@ using Test import GeometryOps as GO, GeoInterface as GI import Extents +using GeometryOps.UnitSpherical point = GI.Point(1.0, 1.0) linestring = GI.LineString([(1.0, 1.0), (2.0, 2.0)]) @@ -190,3 +191,39 @@ end @test_throws ArgumentError collect(GO.lazy_edge_extents(GI.MultiPoint([(1.0, 1.0), (2.0, 2.0)]))) end + +@testset "eachedge and to_edgelist on the sphere" begin + lonlat_ring = GI.LinearRing([(0.0, 0.0), (90.0, 0.0), (0.0, 90.0), (0.0, 0.0)]) + + edges = collect(GO.eachedge(GO.Spherical(), lonlat_ring, Float64)) + @test length(edges) == 3 + @test all(ps -> ps isa Tuple{UnitSphericalPoint{Float64}, UnitSphericalPoint{Float64}}, edges) + # Planar() is the manifold-less behavior + @test collect(GO.eachedge(GO.Planar(), lonlat_ring, Float64)) == collect(GO.eachedge(lonlat_ring, Float64)) + + el = GO.to_edgelist(GO.Spherical(), lonlat_ring) + @test length(el) == 3 + @test all(l -> GI.extent(l) isa Extents.Extent{(:X, :Y, :Z)}, el) + # rings of UnitSphericalPoints take the spherical path with no manifold argument + usp_ring = GI.LinearRing(to_unit_spherical_points(lonlat_ring)) + @test GI.extent.(GO.to_edgelist(usp_ring)) == GI.extent.(el) + # the lazy variant agrees + @test GI.extent.(collect(GO.lazy_edgelist(GO.Spherical(), lonlat_ring))) == GI.extent.(el) + + # spatial trees over spherical edges index in 3D + ni = GO.NaturalIndexing.NaturalIndex(el) + rt = GO.FlexibleRTrees.RTree(GO.FlexibleRTrees.STR(), el) + @test Extents.extent(ni) isa Extents.Extent{(:X, :Y, :Z)} + @test Extents.extent(rt) isa Extents.Extent{(:X, :Y, :Z)} + + # The arc's bulge is indexed: an edge between two points at z = 0.9 arcs + # over the pole, so a query box touching only the polar region must find + # it, even though the endpoints' own box tops out at z = 0.9 + z = 0.9; s = sqrt(1 - z^2) + seg = GI.LineString([UnitSphericalPoint(0.0, -s, z), UnitSphericalPoint(0.0, s, z)]) + tree = GO.FlexibleRTrees.RTree(GO.FlexibleRTrees.STR(), GO.to_edgelist(seg)) + bulge_only = Extents.Extent(X = (-0.01, 0.01), Y = (-0.01, 0.01), Z = (0.95, 1.05)) + @test GO.FlexibleRTrees.query(tree, bulge_only) == [1] + endpoint_box = Extents.union(GI.extent(UnitSphericalPoint(0.0, -s, z)), GI.extent(UnitSphericalPoint(0.0, s, z))) + @test !Extents.intersects(endpoint_box, bulge_only) +end From f75d1de250218f82a7b40c70bc2f3f0301c230de Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 16:32:29 -0700 Subject: [PATCH 11/20] Trim comments and docstrings Co-Authored-By: Claude Fable 5 --- src/utils/UnitSpherical/arc_extent.jl | 8 +++----- src/utils/utils.jl | 14 +++++--------- test/utils/unitspherical.jl | 9 ++++----- test/utils/utils.jl | 7 +++---- 4 files changed, 15 insertions(+), 23 deletions(-) diff --git a/src/utils/UnitSpherical/arc_extent.jl b/src/utils/UnitSpherical/arc_extent.jl index c571b58c85..d59fadf41e 100644 --- a/src/utils/UnitSpherical/arc_extent.jl +++ b/src/utils/UnitSpherical/arc_extent.jl @@ -41,11 +41,9 @@ on the unit sphere. Accepts `UnitSphericalPoint`s, or any GeoInterface point (interpreted geographically, as longitude/latitude, like the `UnitSphericalPoint` constructor itself). -The extent is exact up to floating point error and padded by a few ulps, so -it always contains the arc — unlike the extent of the endpoints, which the -arc bulges out of wherever a coordinate attains its extremum between them. -For antipodal endpoints the arc's plane is ambiguous; the one chosen by -[`robust_cross_product`](@ref) is used. +The extent is exact up to floating point error, padded by a few ulps so it +always contains the arc. For antipodal endpoints the arc's plane is +ambiguous; the one chosen by [`robust_cross_product`](@ref) is used. ## Example diff --git a/src/utils/utils.jl b/src/utils/utils.jl index 7aae0831b2..d270d1d9ee 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -168,9 +168,9 @@ Returns some iterator, which yields tuples of points. Each tuple is an edge. It goes `(p1, p2), (p2, p3), (p3, p4), ...` etc. -On `Planar()` (the manifold-less default) points are 2D coordinate tuples. -On `Spherical()` they are `UnitSphericalPoint`s: geographic (longitude, -latitude) input is converted, and `UnitSphericalPoint`s pass through as-is. +On `Planar()` (the default) points are 2D coordinate tuples. On +`Spherical()` they are `UnitSphericalPoint`s: geographic (longitude, +latitude) input is converted, `UnitSphericalPoint`s pass through. """ eachedge(geom) = eachedge(GI.trait(geom), geom, Float64) function eachedge(geom, ::Type{T}) where T @@ -202,10 +202,8 @@ end Convert a geometry into a vector of `GI.Line` objects with attached extents. On `Spherical()` — or whenever the geometry's points are already -`UnitSphericalPoint`s — each edge is a great-circle arc of -`UnitSphericalPoint`s carrying its 3D [`UnitSpherical.spherical_arc_extent`](@ref), -so spatial indices built over the edges (e.g. `NaturalIndex`, `RTree`) -bound the arcs correctly. +`UnitSphericalPoint`s — each edge carries the 3D +[`UnitSpherical.spherical_arc_extent`](@ref) of its great-circle arc. """ to_edgelist(geom, ::Type{T} = Float64) where T = [_lineedge(ps, T) for ps in eachedge(geom, T)] @@ -240,8 +238,6 @@ function _lineedge(ps::Tuple, ::Type{T}) where T e = GI.extent(l) return GI.Line(l.geom; extent=e) end -# On the unit sphere, an edge is a great-circle arc: its extent is 3D and -# must cover the arc's bulge past the endpoints, not just the endpoints. function _lineedge(ps::Tuple{<:UnitSpherical.UnitSphericalPoint, <:UnitSpherical.UnitSphericalPoint}, ::Type{T}) where T a, b = UnitSpherical.UnitSphericalPoint{T}.(ps) return GI.Line(StaticArrays.SVector((a, b)); extent = UnitSpherical.spherical_arc_extent(a, b)) diff --git a/test/utils/unitspherical.jl b/test/utils/unitspherical.jl index a390cf957a..890b8ebf4b 100644 --- a/test/utils/unitspherical.jl +++ b/test/utils/unitspherical.jl @@ -599,17 +599,16 @@ end @test ext.Z[1] ≈ 0 atol = 1e-14 @test ext.Z[2] ≈ 0 atol = 1e-14 - # The arc bulges out of the endpoints' box: two points at z = 0.9 on - # either side of the x = 0 great circle, joined over the pole + # Two points at z = 0.9 joined over the pole: the arc reaches z = 1 z = 0.9; s = sqrt(1 - z^2) bulging = spherical_arc_extent(UnitSphericalPoint(0.0, -s, z), UnitSphericalPoint(0.0, s, z)) - @test bulging.Z[2] ≈ 1 atol = 1e-14 # the pole, not the endpoints' z = 0.9 + @test bulging.Z[2] ≈ 1 atol = 1e-14 @test bulging.Z[2] > z @test bulging.Z[1] ≈ z atol = 1e-14 @test bulging.Y[1] ≈ -s atol = 1e-14 @test bulging.Y[2] ≈ s atol = 1e-14 - # Geographic (lon, lat) input takes the same path + # Geographic (lon, lat) input @test spherical_arc_extent((0.0, 0.0), (90.0, 0.0)) == spherical_arc_extent(UnitSphericalPoint(1.0, 0.0, 0.0), UnitSphericalPoint(0.0, 1.0, 0.0)) @@ -627,7 +626,7 @@ end a, b = rand(rng, UnitSphericalPoint{Float64}), rand(rng, UnitSphericalPoint{Float64}) e = spherical_arc_extent(a, b) @test all(t -> inext(slerp(a, b, t), e), range(0.0, 1.0, length = 101)) - # nearly-degenerate arcs stay stable through robust_cross_product + # nearly-degenerate arc b2 = UnitSphericalPoint(normalize(a + 1e-12 * rand(rng, UnitSphericalPoint{Float64}))) e2 = spherical_arc_extent(a, b2) @test all(t -> inext(slerp(a, b2, t), e2), range(0.0, 1.0, length = 11)) diff --git a/test/utils/utils.jl b/test/utils/utils.jl index de1648dc18..cf3d9f02d2 100644 --- a/test/utils/utils.jl +++ b/test/utils/utils.jl @@ -198,7 +198,7 @@ end edges = collect(GO.eachedge(GO.Spherical(), lonlat_ring, Float64)) @test length(edges) == 3 @test all(ps -> ps isa Tuple{UnitSphericalPoint{Float64}, UnitSphericalPoint{Float64}}, edges) - # Planar() is the manifold-less behavior + # Planar() matches the manifold-less form @test collect(GO.eachedge(GO.Planar(), lonlat_ring, Float64)) == collect(GO.eachedge(lonlat_ring, Float64)) el = GO.to_edgelist(GO.Spherical(), lonlat_ring) @@ -216,9 +216,8 @@ end @test Extents.extent(ni) isa Extents.Extent{(:X, :Y, :Z)} @test Extents.extent(rt) isa Extents.Extent{(:X, :Y, :Z)} - # The arc's bulge is indexed: an edge between two points at z = 0.9 arcs - # over the pole, so a query box touching only the polar region must find - # it, even though the endpoints' own box tops out at z = 0.9 + # An edge between two points at z = 0.9 arcs over the pole; a query box + # touching only the polar region must still find it z = 0.9; s = sqrt(1 - z^2) seg = GI.LineString([UnitSphericalPoint(0.0, -s, z), UnitSphericalPoint(0.0, s, z)]) tree = GO.FlexibleRTrees.RTree(GO.FlexibleRTrees.STR(), GO.to_edgelist(seg)) From 34c156eca8742f4fcd3b80e610374d3ead2e9cc1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 16:34:47 -0700 Subject: [PATCH 12/20] Trim the arc extent literate header Co-Authored-By: Claude Fable 5 --- src/utils/UnitSpherical/arc_extent.jl | 30 +++++++++------------------ 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/utils/UnitSpherical/arc_extent.jl b/src/utils/UnitSpherical/arc_extent.jl index d59fadf41e..9577f616f5 100644 --- a/src/utils/UnitSpherical/arc_extent.jl +++ b/src/utils/UnitSpherical/arc_extent.jl @@ -5,32 +5,22 @@ spherical_arc_extent ``` -## Why not the extent of the endpoints? - A great-circle arc bulges away from the chord between its endpoints, so the -axis-aligned bounding box of the endpoints does not, in general, contain the -arc. The classic case is two points at the same latitude: the arc between -them passes closer to the pole than either endpoint, e.g. two points at +endpoints' bounding box does not in general contain the arc: two points at `z = 0.9` on either side of the prime meridian are joined by an arc whose -midpoint has `z = 0.9 / cos(θ/2) > 0.9`. A spatial index built on endpoint -boxes would silently miss queries that touch only the bulge. - -## How the extremum is found +midpoint has `z = 0.9 / cos(θ/2) > 0.9`. [`spherical_arc_extent`](@ref) +computes a box that contains the whole arc. With `t̂ₐ` the unit tangent at `a` pointing along the arc, the arc is `p(φ) = a cos(φ) + t̂ₐ sin(φ)` for `φ ∈ [0, θ]`, so each Cartesian coordinate is a sinusoid `pᵢ(φ) = Rᵢ cos(φ - φᵢ)` with amplitude -`Rᵢ = hypot(aᵢ, t̂ₐᵢ)`. Since `θ ≤ π`, at most one interior maximum and one -interior minimum exist per axis, and the endpoint derivatives decide: an -interior maximum exists iff `pᵢ` is increasing at `a` and decreasing at `b` -(`t̂ₐᵢ > 0 > t̂ᵦᵢ`), where it attains `Rᵢ`; likewise a minimum attains `-Rᵢ`. -No trigonometric calls are needed, and the tangents come from -[`robust_cross_product`](@ref), so nearly-degenerate and nearly-antipodal -arcs stay stable. - -Bounds are padded by a few ulps so that the extent is guaranteed to contain -the arc despite floating point error, in the same spirit as S2's -`S2LatLngRectBounder`, which widens its bounds by their maximum error. +`Rᵢ = hypot(aᵢ, t̂ₐᵢ)`. Since `θ ≤ π` there is at most one interior maximum +and one interior minimum per axis: a maximum exists iff `pᵢ` increases at +`a` and decreases at `b` (`t̂ₐᵢ > 0 > t̂ᵦᵢ`), where it attains `Rᵢ`; minima +mirror. The tangents come from [`robust_cross_product`](@ref), which keeps +nearly-degenerate and nearly-antipodal arcs stable. Bounds are padded by a +few ulps to absorb floating point error, as S2's `S2LatLngRectBounder` pads +by its maximum error. =# """ From 528ac9e1dd783e3f68bfb820154d75d49b8ba229 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 6 Jul 2026 19:54:01 -0700 Subject: [PATCH 13/20] Add manifold-aware `Extents.extent` with spherical region extents MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `Extents.extent(m::Manifold, geom, [T])` computes extents on a manifold: `Planar()` delegates to `GI.extent`; `Spherical()` returns 3D Cartesian extents on the unit sphere. Curves are covered by the union of their edges' arc extents. Rings and polygons are regions under S2's loop convention (CCW, interior on the left): an extremum over a region lies on the boundary or at one of the six axis points `±eᵢ`, so the boundary extent is widened per axis by an enclosure check from the ring's winding number and signed area. GO does not export `extent`; the methods extend `Extents.extent`, reachable as `GO.extent`. Also fix `to_unit_spherical_points` to pass `UnitSphericalPoint`s through unchanged instead of reinterpreting their first two Cartesian coordinates as geographic longitude/latitude. Co-Authored-By: Claude Fable 5 --- src/GeometryOps.jl | 1 + src/methods/extent.jl | 125 +++++++++++++++++++++++ src/utils/UnitSpherical/UnitSpherical.jl | 10 +- test/methods/extent.jl | 117 +++++++++++++++++++++ test/runtests.jl | 1 + 5 files changed, 249 insertions(+), 5 deletions(-) create mode 100644 src/methods/extent.jl create mode 100644 test/methods/extent.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 58d380eb34..02fcd11934 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -69,6 +69,7 @@ include("methods/centroid.jl") include("methods/convex_hull.jl") include("methods/distance.jl") include("methods/equals.jl") +include("methods/extent.jl") include("methods/perimeter.jl") include("methods/clipping/predicates.jl") include("methods/clipping/clipping_processor.jl") diff --git a/src/methods/extent.jl b/src/methods/extent.jl new file mode 100644 index 0000000000..b5d8c9db06 --- /dev/null +++ b/src/methods/extent.jl @@ -0,0 +1,125 @@ +# # Extent + +#= +```@docs; canonical=false +Extents.extent(::Manifold, ::Any) +``` + +`Extents.extent(m::Manifold, geom)` computes the extent of a geometry *on a +manifold*. On `Planar()` it delegates to `GI.extent`. On `Spherical()` it +returns 3D Cartesian `Extent{(:X, :Y, :Z)}`s on the unit sphere — the boxes +that spatial indices over spherical geometry prune with. + +On the sphere, an extremum of a coordinate over a *region* is attained +either on the boundary or at one of the six on-sphere critical points of +that coordinate: `(±1,0,0)`, `(0,±1,0)`, `(0,0,±1)`. A polygon that +strictly encloses one of them (a cell over a pole, say) therefore extends +past every boundary edge's extent, so a region's box is the union of its +edges' [`UnitSpherical.spherical_arc_extent`](@ref)s plus an enclosure +check per axis point. + +Enclosure follows S2's loop convention (`s2loop.h`): "All loops are defined +to have a CCW orientation, i.e. the interior of the loop is on the left +side of the edges. This implies that a clockwise loop enclosing a small +area is interpreted to be a CCW loop enclosing a very large area." A ring +encloses exactly one of `±eᵢ` iff its winding number about that axis is +`±1`, and the sign picks which; both are enclosed only when the interior is +the larger side of the ring (negative signed area, or area above `2π`), in +which case the axis is extended to `[-1, 1]`. + +Two documented approximations, both conservative-safe for meshes: the +winding accumulates wrapped angle deltas, so a single edge must not sweep +more than a half turn about an axis (an edge passing closer to `±eᵢ` than +roughly its own length) — such an edge's own arc extent already reaches +within `(distance)²/2` of `±1`. And a region whose interior strictly +contains an antipodal pair away from the axis points (a thin tube pole to +pole) is beyond the winding test; its boundary extents again come within +`(distance)²/2` of the truth. +=# + +""" + extent(m::Manifold, geom, [::Type{T} = Float64])::Extents.Extent + +The extent of `geom` on the manifold `m` — this method lives on, and +returns an, `Extents.Extent`. + +On `Planar()`, `GI.extent(geom)`. On `Spherical()`, the 3D Cartesian +extent of the geometry on the unit sphere, with geographic (longitude, +latitude) input converted like `UnitSphericalPoint`: curves are covered by +the union of their edges' great-circle arc extents, and rings and polygons +are treated as regions — wound CCW with the interior on the left, per S2's +loop convention — whose extent also covers any enclosed pole or other +on-sphere axis extreme. + +## Example + +```jldoctest +julia> import GeometryOps as GO, GeoInterface as GI + +julia> cap = GI.Polygon([[(lon, 60.0) for lon in 0.0:30.0:360.0]]); # around the north pole + +julia> GO.extent(GO.Spherical(), cap).Z[2] +1.0 +``` +""" +function Extents.extent(m::Manifold, geom, ::Type{T} = Float64) where T + return _extent(m, GI.trait(geom), geom, T) +end + +_extent(::Planar, trait, geom, ::Type{T}) where T = GI.extent(geom) + +_extent(::Spherical, ::GI.PointTrait, geom, ::Type{T}) where T = + GI.extent(UnitSpherical.UnitSphericalPoint(geom)) +_extent(m::Spherical, ::Union{GI.LineTrait, GI.LineStringTrait}, geom, ::Type{T}) where T = + mapreduce(GI.extent, Extents.union, lazy_edgelist(m, geom, T)) +_extent(m::Spherical, ::GI.LinearRingTrait, geom, ::Type{T}) where T = + _spherical_region_extent(UnitSpherical.to_unit_spherical_points(geom)) +_extent(m::Spherical, ::GI.PolygonTrait, geom, ::Type{T}) where T = + _extent(m, GI.LinearRingTrait(), GI.getexterior(geom), T) +# multi-geometries and collections; holes never extend a polygon's extent, +# so only the exterior ring above matters +_extent(m::Spherical, ::GI.AbstractGeometryTrait, geom, ::Type{T}) where T = + mapreduce(g -> Extents.extent(m, g, T), Extents.union, GI.getgeom(geom)) + +function _spherical_region_extent(pts::Vector{<:UnitSpherical.UnitSphericalPoint}) + n = length(pts) + n > 1 && pts[end] == pts[1] && (n -= 1) + ext = mapreduce(Extents.union, 1:n) do i + UnitSpherical.spherical_arc_extent(pts[i], pts[mod1(i + 1, n)]) + end + n < 3 && return ext + + # winding about each axis from wrapped angle deltas in the plane + # perpendicular to it; a vertex on an axis makes that winding + # meaningless, but also puts ±1 into the edge extents above, so skip it + winding = zeros(MVector{3, Float64}) + onaxis = MVector(false, false, false) + angles(p) = (atan(p.z, p.y), atan(p.x, p.z), atan(p.y, p.x)) + prev = angles(pts[n]) + for i in 1:n + p = pts[i] + onaxis[1] |= p.y == 0 && p.z == 0 + onaxis[2] |= p.z == 0 && p.x == 0 + onaxis[3] |= p.x == 0 && p.y == 0 + cur = angles(p) + winding .+= rem.(cur .- prev, 2π, RoundNearest) + prev = cur + end + + ringarea = sum(i -> _spherical_triangle_area(Girard(), pts[1], pts[i], pts[i + 1]), 2:(n - 1); init = 0.0) + bigregion = ringarea < 0 || ringarea > 2π + + axis_bounds = values(ext) + bounds = ntuple(3) do i + lo, hi = axis_bounds[i] + if !onaxis[i] + winding[i] > π && (hi = one(hi)) + winding[i] < -π && (lo = -one(lo)) + if bigregion && abs(winding[i]) <= π + lo, hi = -one(lo), one(hi) + end + end + (lo, hi) + end + return Extents.Extent(X = bounds[1], Y = bounds[2], Z = bounds[3]) +end diff --git a/src/utils/UnitSpherical/UnitSpherical.jl b/src/utils/UnitSpherical/UnitSpherical.jl index cdc70b62e5..eba53275c2 100644 --- a/src/utils/UnitSpherical/UnitSpherical.jl +++ b/src/utils/UnitSpherical/UnitSpherical.jl @@ -31,14 +31,14 @@ export UnitSphericalPoint, UnitSphereFromGeographic, GeographicFromUnitSphere, to_unit_spherical_points """ - to_unit_spherical_points(ring) -> Vector{UnitSphericalPoint{Float64}} + to_unit_spherical_points(ring) -> Vector{<:UnitSphericalPoint} -Convert a ring (linear ring or any GeoInterface point iterator) to a vector of UnitSphericalPoints. -Uses UnitSphereFromGeographic which is a no-op for already-converted points. +Convert a ring (linear ring or any GeoInterface point iterator) to a vector of +UnitSphericalPoints, treating geographic input as (longitude, latitude). +`UnitSphericalPoint`s pass through unchanged. """ function to_unit_spherical_points(ring) - transform = UnitSphereFromGeographic() - return [transform((GI.x(p), GI.y(p))) for p in GI.getpoint(ring)] + return [UnitSphericalPoint(p) for p in GI.getpoint(ring)] end end \ No newline at end of file diff --git a/test/methods/extent.jl b/test/methods/extent.jl new file mode 100644 index 0000000000..65fd2022d4 --- /dev/null +++ b/test/methods/extent.jl @@ -0,0 +1,117 @@ +using Test +using LinearAlgebra + +import GeometryOps as GO, GeoInterface as GI +import Extents +using GeometryOps.UnitSpherical +using Random: Xoshiro + +# k vertices of the z = z₀ circle, CCW seen from +z (the S2 interior-on-left +# convention: the enclosed region is the cap containing the north pole) +polar_ring(z, k) = [UnitSphericalPoint(sqrt(1 - z^2) * cos(t), sqrt(1 - z^2) * sin(t), z) + for t in range(0, 2π; length = k + 1)[1:(end - 1)]] + +@testset "extent(Planar(), ...)" begin + poly = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]]) + @test GO.extent(GO.Planar(), poly) == GI.extent(poly) + @test GO.extent(GO.Planar(), GI.Point(1.0, 2.0)) == GI.extent(GI.Point(1.0, 2.0)) +end + +@testset "extent(Spherical(), ...)" begin + m = GO.Spherical() + z = 0.9; s = sqrt(1 - z^2) + + @testset "CCW polar cap ring encloses the pole" begin + ext = GO.extent(m, GI.LinearRing(polar_ring(z, 8))) + @test ext isa Extents.Extent{(:X, :Y, :Z)} + @test ext.Z[2] == 1 + @test ext.Z[1] ≈ z atol = 1e-12 + @test -1 < ext.X[1] && ext.X[2] < 1 + @test -1 < ext.Y[1] && ext.Y[2] < 1 + end + + @testset "CW ring is the complement (S2 convention)" begin + ext = GO.extent(m, GI.LinearRing(reverse(polar_ring(z, 8)))) + @test ext.Z[1] == -1 + @test ext.Z[2] < 1 # sup z of the complement is on the boundary + @test ext.X == (-1.0, 1.0) # complement contains ±eₓ and ±e_y + @test ext.Y == (-1.0, 1.0) + end + + @testset "Geographic polygon around the pole" begin + cap = GI.Polygon([[(lon, 60.0) for lon in 0.0:30.0:360.0]]) + ext = GO.extent(m, cap) + @test ext.Z[2] == 1 + @test ext.Z[1] ≈ sind(60) atol = 1e-12 + end + + @testset "No enclosure: region extent equals curve extent" begin + pts = [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)] + @test GO.extent(m, GI.Polygon([pts])) == GO.extent(m, GI.LineString(pts)) + end + + @testset "Vertex exactly at the pole" begin + ring = GI.LinearRing([UnitSphericalPoint(0.0, 0.0, 1.0), + UnitSphericalPoint(s, 0.0, z), + UnitSphericalPoint(0.0, s, z), + UnitSphericalPoint(0.0, 0.0, 1.0)]) + ext = GO.extent(m, ring) + @test ext.Z[2] >= 1 + @test all(b -> all(isfinite, b), values(ext)) + end + + @testset "Cap larger than a hemisphere" begin + ext = GO.extent(m, GI.LinearRing(polar_ring(cosd(100), 16))) + @test ext.Z[2] == 1 + @test ext.X == (-1.0, 1.0) # ±eₓ and ±e_y are interior + @test ext.Y == (-1.0, 1.0) + @test ext.Z[1] < cosd(100) # arcs bulge below the vertex circle + end + + @testset "Points, multis, and collections" begin + p = GI.Point(0.0, 0.0) # lon/lat → (1, 0, 0) + ep = GO.extent(m, p) + @test ep.X[1] == ep.X[2] == 1.0 + emp = GO.extent(m, GI.MultiPoint([(0.0, 0.0), (90.0, 0.0)])) + @test emp.X == (0.0, 1.0) && emp.Y == (0.0, 1.0) + + north = GI.Polygon([[(lon, 60.0) for lon in 0.0:30.0:360.0]]) + south = GI.Polygon([[(lon, -60.0) for lon in 360.0:-30.0:0.0]]) # CCW around the south pole + ext = GO.extent(m, GI.MultiPolygon([north, south])) + @test ext.Z == (-1.0, 1.0) + end + + @testset "LineStrings are curves, not regions" begin + ext = GO.extent(m, GI.LineString(vcat(polar_ring(z, 8), [polar_ring(z, 8)[1]]))) + @test ext.Z[2] < 1 # no interior, no pole + end + + @testset "Random cells contain their samples and enclosed axis points" begin + rng = Xoshiro(2026) + inext(q, e) = e.X[1] <= q[1] <= e.X[2] && e.Y[1] <= q[2] <= e.Y[2] && e.Z[1] <= q[3] <= e.Z[2] + axispoints = [UnitSphericalPoint(v...) for v in + ((1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1))] + for _ in 1:50 + c = rand(rng, UnitSphericalPoint{Float64}) + r = 0.05 + 0.95 * rand(rng) + u = normalize(cross(c, abs(c.z) < 0.9 ? UnitSphericalPoint(0.0, 0.0, 1.0) : UnitSphericalPoint(1.0, 0.0, 0.0))) + v = cross(c, u) + k = rand(rng, 3:8) + ring = [UnitSphericalPoint(cos(r) * c + sin(r) * (cos(t) * u + sin(t) * v)) + for t in range(0, 2π; length = k + 1)[1:(end - 1)]] + ext = GO.extent(m, GI.LinearRing(ring)) + # boundary and interior samples lie inside + samples = [slerp(ring[i], ring[mod1(i + 1, k)], t) + for i in 1:k for t in range(0.0, 1.0; length = 33)] + @test all(q -> inext(q, ext), samples) + @test inext(c, ext) + # the polygon is star-shaped around c, so it contains the cap + # whose radius is the boundary's least distance to c — any axis + # point in that cap must be covered + r_in = 0.98 * minimum(q -> spherical_distance(c, q), samples) + for a in axispoints + spherical_distance(c, a) < r_in && @test inext(a, ext) + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index f67aceee5f..8540d6b741 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -34,6 +34,7 @@ end @safetestset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end @safetestset "Distance" begin include("methods/distance.jl") end @safetestset "Equals" begin include("methods/equals.jl") end +@safetestset "Extent" begin include("methods/extent.jl") end @safetestset "Minimum Bounding Circle" begin include("methods/minimum_bounding_circle.jl") end @safetestset "Polygonize" begin include("methods/polygonize.jl") end # Clipping From 5ee72bec316d1bf91855e35bef2923da5ccaa10b Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 7 Jul 2026 17:34:40 -0700 Subject: [PATCH 14/20] Replace the winding enclosure heuristic with S2-style crossing parity Axis-point enclosure in `_spherical_region_extent` is now decided the way `S2Loop::InitBound` decides pole containment: the side of an anchor edge gives the departure side of the arc from the query point to the anchor's midpoint, and the parity of transversal boundary crossings flips it. Degenerate configurations retry with the next anchor and fall back to a conservative extension, never an under-covering box. This removes both documented approximations of the winding approach: a region containing both poles at small area (dumbbell) and edges sweeping near an axis are now handled by construction. Co-Authored-By: Claude Fable 5 --- src/methods/extent.jl | 136 ++++++++++++++++++++++++++++------------- test/methods/extent.jl | 43 ++++++++++++- 2 files changed, 136 insertions(+), 43 deletions(-) diff --git a/src/methods/extent.jl b/src/methods/extent.jl index b5d8c9db06..9e3f110dd1 100644 --- a/src/methods/extent.jl +++ b/src/methods/extent.jl @@ -21,20 +21,24 @@ check per axis point. Enclosure follows S2's loop convention (`s2loop.h`): "All loops are defined to have a CCW orientation, i.e. the interior of the loop is on the left side of the edges. This implies that a clockwise loop enclosing a small -area is interpreted to be a CCW loop enclosing a very large area." A ring -encloses exactly one of `±eᵢ` iff its winding number about that axis is -`±1`, and the sign picks which; both are enclosed only when the interior is -the larger side of the ring (negative signed area, or area above `2π`), in -which case the axis is extended to `[-1, 1]`. - -Two documented approximations, both conservative-safe for meshes: the -winding accumulates wrapped angle deltas, so a single edge must not sweep -more than a half turn about an axis (an edge passing closer to `±eᵢ` than -roughly its own length) — such an edge's own arc extent already reaches -within `(distance)²/2` of `±1`. And a region whose interior strictly -contains an antipodal pair away from the axis points (a thin tube pole to -pole) is beyond the winding test; its boundary extents again come within -`(distance)²/2` of the truth. +area is interpreted to be a CCW loop enclosing a very large area." + +The enclosure test is crossing parity, the way `S2Loop::InitBound` decides +pole containment (`s2loop.cc`). Pick an anchor edge whose great circle +does not pass through the query point `q`: which side of that edge `q` +lies on says whether the arc from the edge's midpoint to `q` *departs* +into the interior (the left side) or the exterior, and each transversal +boundary crossing along the arc flips that. The departure side is exactly +`q`'s side because the arc can meet the anchor's great circle again only +at the midpoint's antipode, which an arc shorter than a half turn never +reaches. + +Where S2 resolves degenerate configurations with exact predicates and +symbolic perturbation, this test detects them — a vertex within +[`UnitSpherical.spherical_orient`](@ref)'s tolerance of a test arc's great +circle, a crossing too close to an arc endpoint to call — and retries with +the next edge as anchor. If every anchor is degenerate the axis is +extended to `±1`, so the box can come out loose but never under-covers. =# """ @@ -89,37 +93,85 @@ function _spherical_region_extent(pts::Vector{<:UnitSpherical.UnitSphericalPoint end n < 3 && return ext - # winding about each axis from wrapped angle deltas in the plane - # perpendicular to it; a vertex on an axis makes that winding - # meaningless, but also puts ±1 into the edge extents above, so skip it - winding = zeros(MVector{3, Float64}) - onaxis = MVector(false, false, false) - angles(p) = (atan(p.z, p.y), atan(p.x, p.z), atan(p.y, p.x)) - prev = angles(pts[n]) - for i in 1:n - p = pts[i] - onaxis[1] |= p.y == 0 && p.z == 0 - onaxis[2] |= p.z == 0 && p.x == 0 - onaxis[3] |= p.x == 0 && p.y == 0 - cur = angles(p) - winding .+= rem.(cur .- prev, 2π, RoundNearest) - prev = cur + lo = MVector(ext.X[1], ext.Y[1], ext.Z[1]) + hi = MVector(ext.X[2], ext.Y[2], ext.Z[2]) + for i in 1:3, s in (1.0, -1.0) + q = UnitSpherical.UnitSphericalPoint(ntuple(j -> j == i ? s : 0.0, 3)) + inside = _spherical_ring_contains(pts, n, q) + # nothing = undecidable: extend anyway so the box never under-covers + if inside === nothing || inside + s > 0 ? (hi[i] = one(hi[i])) : (lo[i] = -one(lo[i])) + end end + return Extents.Extent(X = (lo[1], hi[1]), Y = (lo[2], hi[2]), Z = (lo[3], hi[3])) +end - ringarea = sum(i -> _spherical_triangle_area(Girard(), pts[1], pts[i], pts[i + 1]), 2:(n - 1); init = 0.0) - bigregion = ringarea < 0 || ringarea > 2π - - axis_bounds = values(ext) - bounds = ntuple(3) do i - lo, hi = axis_bounds[i] - if !onaxis[i] - winding[i] > π && (hi = one(hi)) - winding[i] < -π && (lo = -one(lo)) - if bigregion && abs(winding[i]) <= π - lo, hi = -one(lo), one(hi) +# Crossing-parity containment of `q` in the closed region left of the ring, +# after S2Loop::Contains/InitBound. Returns `nothing` when every anchor edge +# is degenerate with respect to `q`. +function _spherical_ring_contains(pts, n, q) + for j in 1:n + UnitSpherical.point_on_spherical_arc(q, pts[j], pts[mod1(j + 1, n)]) && return true + end + for j in 1:n + a, b = pts[j], pts[mod1(j + 1, n)] + a == b && continue + side = UnitSpherical.spherical_orient(a, b, q) + side == 0 && continue + mid = a + b + norm(mid) < 1e-9 && continue # near-antipodal edge, midpoint unstable + m = UnitSpherical.UnitSphericalPoint(normalize(mid)) + dot(q, m) < -1 + 1e-9 && continue # test arc q → m would span a half turn + crossings = 0 + ok = true + for k in 1:n + k == j && continue + c = _arc_crossing_parity(q, m, pts[k], pts[mod1(k + 1, n)]) + if c == -1 + ok = false + break end + crossings += c end - (lo, hi) + ok || continue + # walking from `m` toward `q` departs onto `q`'s side of the anchor + # edge (the arc meets that great circle again only at `-m`); positive + # side is the interior, and each crossing flips it + return isodd(crossings) ? side < 0 : side > 0 end - return Extents.Extent(X = bounds[1], Y = bounds[2], Z = bounds[3]) + return nothing +end + +# Crossing parity of the test arc q → m against ring edge a → b: 1 for a +# transversal crossing, 0 for none, -1 for too close to degenerate to call. +function _arc_crossing_parity(q, m, a, b) + # a vertex exactly antipodal to `q` lies on every great circle through + # `q`, but its edges can reach the test arc only at `q` itself, which + # the on-boundary check has already excluded + (a == -q || b == -q) && return 0 + sa = UnitSpherical.spherical_orient(q, m, a) + sb = UnitSpherical.spherical_orient(q, m, b) + (sa == 0 || sb == 0) && return -1 + sa == sb && return 0 + # `q` on this edge's great circle (but not on the edge — checked + # upfront): the two circles meet only at `±q`, and the test arc reaches + # neither, so the edge cannot cross it. This is systematic, not rare — + # a lonlat grid's meridian edges pass through `±eₓ`/`±e_y` exactly — + # and no anchor changes it, so it must resolve rather than retry. + sq = UnitSpherical.spherical_orient(a, b, q) + sq == 0 && return 0 + sm = UnitSpherical.spherical_orient(a, b, m) + sm == 0 && return -1 + sq == sm && return 0 + # each arc now crosses the other's great circle exactly once, at one of + # the two antipodal circle intersections; the arcs cross iff those are + # the same point, i.e. iff the intersection direction `x` points into + # both arcs' hemispheres + x = cross(normalize(UnitSpherical.robust_cross_product(q, m)), + normalize(UnitSpherical.robust_cross_product(a, b))) + d1 = dot(x, q + m) + d2 = dot(x, a + b) + tol = 16 * eps(Float64) * norm(x) + (abs(d1) <= tol || abs(d2) <= tol) && return -1 + return (d1 > 0) == (d2 > 0) ? 1 : 0 end diff --git a/test/methods/extent.jl b/test/methods/extent.jl index 65fd2022d4..aa0f02a0eb 100644 --- a/test/methods/extent.jl +++ b/test/methods/extent.jl @@ -46,10 +46,18 @@ end end @testset "No enclosure: region extent equals curve extent" begin - pts = [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)] + pts = [(5.0, 5.0), (15.0, 5.0), (15.0, 15.0), (5.0, 15.0), (5.0, 5.0)] @test GO.extent(m, GI.Polygon([pts])) == GO.extent(m, GI.LineString(pts)) end + @testset "Axis point on the boundary clamps, and extends nothing else" begin + pts = [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)] # corner at +eₓ + ext = GO.extent(m, GI.Polygon([pts])) + @test ext.X[2] == 1.0 + @test ext.X[1] > 0.9 && ext.Y[2] < 0.2 && ext.Z[2] < 0.2 + @test ext.Y[1] > -0.1 && ext.Z[1] > -0.1 + end + @testset "Vertex exactly at the pole" begin ring = GI.LinearRing([UnitSphericalPoint(0.0, 0.0, 1.0), UnitSphericalPoint(s, 0.0, z), @@ -68,6 +76,30 @@ end @test ext.Z[1] < cosd(100) # arcs bulge below the vertex circle end + @testset "Dumbbell: both poles through a thin corridor" begin + # two polar caps (above ±85°) joined by a corridor over lon ∈ (355°, 5°); + # small area, no winding about any axis — only a containment test sees + # that both poles and (1, 0, 0) are interior + north = [(lon, 85.0) for lon in range(5.0, 355.0; length = 15)] # eastward, pole on the left + south = [(lon, -85.0) for lon in range(355.0, 5.0; length = 15)] # westward, pole on the left + ext = GO.extent(m, GI.Polygon([vcat(north, south, [north[1]])])) + @test ext.Z == (-1.0, 1.0) + @test ext.X[2] == 1.0 # (1, 0, 0) is inside the corridor + @test ext.X[1] > -1 # (-1, 0, 0) is outside + @test -0.2 < ext.Y[1] && ext.Y[2] < 0.2 + end + + @testset "Lonlat polar cells: exact pole vertex, no far-pole leak" begin + for lon0 in 0.0:30.0:330.0 + cell = GI.Polygon([[(lon0, 80.0), (lon0 + 30.0, 80.0), + (lon0 + 30.0, 90.0), (lon0, 90.0), (lon0, 80.0)]]) + ext = GO.extent(m, cell) + @test ext.Z[2] == 1.0 # the pole is a vertex + @test ext.Z[1] > 0.9 # the south pole must not leak in + @test ext.X[1] > -1 && ext.X[2] < 1 && ext.Y[1] > -1 && ext.Y[2] < 1 + end + end + @testset "Points, multis, and collections" begin p = GI.Point(0.0, 0.0) # lon/lat → (1, 0, 0) ep = GO.extent(m, p) @@ -112,6 +144,15 @@ end for a in axispoints spherical_distance(c, a) < r_in && @test inext(a, ext) end + # the reversed ring bounds the complement (S2 convention): it + # shares the boundary, covers axis points outside the cap, and + # between them the two boxes cover every axis point + rext = GO.extent(m, GI.LinearRing(reverse(ring))) + @test all(q -> inext(q, rext), samples) + for a in axispoints + spherical_distance(c, a) > r && @test inext(a, rext) + @test inext(a, ext) || inext(a, rext) + end end end end From 51d6cd66301aecb8c5d7ee931a59f3fb475cd083 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 7 Jul 2026 17:50:44 -0700 Subject: [PATCH 15/20] Trim comments in the region extent parity code Co-Authored-By: Claude Fable 5 --- src/methods/extent.jl | 31 +++++++++++++++---------------- test/methods/extent.jl | 6 +++--- 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/src/methods/extent.jl b/src/methods/extent.jl index 9e3f110dd1..66ad5c2da4 100644 --- a/src/methods/extent.jl +++ b/src/methods/extent.jl @@ -24,14 +24,13 @@ side of the edges. This implies that a clockwise loop enclosing a small area is interpreted to be a CCW loop enclosing a very large area." The enclosure test is crossing parity, the way `S2Loop::InitBound` decides -pole containment (`s2loop.cc`). Pick an anchor edge whose great circle -does not pass through the query point `q`: which side of that edge `q` -lies on says whether the arc from the edge's midpoint to `q` *departs* -into the interior (the left side) or the exterior, and each transversal -boundary crossing along the arc flips that. The departure side is exactly -`q`'s side because the arc can meet the anchor's great circle again only -at the midpoint's antipode, which an arc shorter than a half turn never -reaches. +pole containment (`s2loop.cc`). For an anchor edge whose great circle +misses the query point `q`, the side of that edge `q` falls on is the side +the arc from the edge's midpoint to `q` departs into — left is the +interior — and each transversal boundary crossing along the arc flips it. +The departure side equals `q`'s side because the arc can meet the anchor's +great circle again only at the midpoint's antipode, which an arc shorter +than a half turn never reaches. Where S2 resolves degenerate configurations with exact predicates and symbolic perturbation, this test detects them — a vertex within @@ -98,7 +97,7 @@ function _spherical_region_extent(pts::Vector{<:UnitSpherical.UnitSphericalPoint for i in 1:3, s in (1.0, -1.0) q = UnitSpherical.UnitSphericalPoint(ntuple(j -> j == i ? s : 0.0, 3)) inside = _spherical_ring_contains(pts, n, q) - # nothing = undecidable: extend anyway so the box never under-covers + # nothing (undecidable) extends too; the box must never under-cover if inside === nothing || inside s > 0 ? (hi[i] = one(hi[i])) : (lo[i] = -one(lo[i])) end @@ -146,18 +145,18 @@ end # transversal crossing, 0 for none, -1 for too close to degenerate to call. function _arc_crossing_parity(q, m, a, b) # a vertex exactly antipodal to `q` lies on every great circle through - # `q`, but its edges can reach the test arc only at `q` itself, which - # the on-boundary check has already excluded + # `q`; its edges can reach the test arc only at `q` itself, excluded by + # the on-boundary check (a == -q || b == -q) && return 0 sa = UnitSpherical.spherical_orient(q, m, a) sb = UnitSpherical.spherical_orient(q, m, b) (sa == 0 || sb == 0) && return -1 sa == sb && return 0 - # `q` on this edge's great circle (but not on the edge — checked - # upfront): the two circles meet only at `±q`, and the test arc reaches - # neither, so the edge cannot cross it. This is systematic, not rare — - # a lonlat grid's meridian edges pass through `±eₓ`/`±e_y` exactly — - # and no anchor changes it, so it must resolve rather than retry. + # `q` on this edge's great circle but off the edge (checked upfront): + # the circles meet only at `±q`, both out of the test arc's reach — no + # crossing. This degeneracy is anchor-independent (a lonlat grid's + # meridian edges hold `±eₓ`/`±e_y` exactly), so it resolves instead of + # returning -1. sq = UnitSpherical.spherical_orient(a, b, q) sq == 0 && return 0 sm = UnitSpherical.spherical_orient(a, b, m) diff --git a/test/methods/extent.jl b/test/methods/extent.jl index aa0f02a0eb..8868ea2cd4 100644 --- a/test/methods/extent.jl +++ b/test/methods/extent.jl @@ -77,9 +77,9 @@ end end @testset "Dumbbell: both poles through a thin corridor" begin - # two polar caps (above ±85°) joined by a corridor over lon ∈ (355°, 5°); - # small area, no winding about any axis — only a containment test sees - # that both poles and (1, 0, 0) are interior + # two polar caps (above ±85°) joined by a corridor over lon ∈ (355°, 5°): + # small area, zero net winding about every axis, and both poles and + # (1, 0, 0) interior north = [(lon, 85.0) for lon in range(5.0, 355.0; length = 15)] # eastward, pole on the left south = [(lon, -85.0) for lon in range(355.0, 5.0; length = 15)] # westward, pole on the left ext = GO.extent(m, GI.Polygon([vcat(north, south, [north[1]])])) From c21f4c3162bb59aa7eb09ba2201e66f4547c86b8 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 11 Jul 2026 15:46:46 -0400 Subject: [PATCH 16/20] Condense spherical extent prose Tighten the literate headers, docstrings, and comments of `spherical_arc_extent` and `Extents.extent(m::Manifold, geom)`; fix an overstated dumbbell test comment (zero net winding holds about the polar axis, not every axis). Co-Authored-By: Claude Fable 5 --- src/methods/extent.jl | 73 ++++++++++----------------- src/utils/UnitSpherical/arc_extent.jl | 28 +++++----- test/methods/extent.jl | 2 +- 3 files changed, 41 insertions(+), 62 deletions(-) diff --git a/src/methods/extent.jl b/src/methods/extent.jl index 66ad5c2da4..47b35c092c 100644 --- a/src/methods/extent.jl +++ b/src/methods/extent.jl @@ -10,49 +10,35 @@ manifold*. On `Planar()` it delegates to `GI.extent`. On `Spherical()` it returns 3D Cartesian `Extent{(:X, :Y, :Z)}`s on the unit sphere — the boxes that spatial indices over spherical geometry prune with. -On the sphere, an extremum of a coordinate over a *region* is attained -either on the boundary or at one of the six on-sphere critical points of -that coordinate: `(±1,0,0)`, `(0,±1,0)`, `(0,0,±1)`. A polygon that -strictly encloses one of them (a cell over a pole, say) therefore extends -past every boundary edge's extent, so a region's box is the union of its -edges' [`UnitSpherical.spherical_arc_extent`](@ref)s plus an enclosure -check per axis point. - -Enclosure follows S2's loop convention (`s2loop.h`): "All loops are defined -to have a CCW orientation, i.e. the interior of the loop is on the left -side of the edges. This implies that a clockwise loop enclosing a small -area is interpreted to be a CCW loop enclosing a very large area." - -The enclosure test is crossing parity, the way `S2Loop::InitBound` decides -pole containment (`s2loop.cc`). For an anchor edge whose great circle -misses the query point `q`, the side of that edge `q` falls on is the side -the arc from the edge's midpoint to `q` departs into — left is the -interior — and each transversal boundary crossing along the arc flips it. -The departure side equals `q`'s side because the arc can meet the anchor's -great circle again only at the midpoint's antipode, which an arc shorter -than a half turn never reaches. - -Where S2 resolves degenerate configurations with exact predicates and -symbolic perturbation, this test detects them — a vertex within -[`UnitSpherical.spherical_orient`](@ref)'s tolerance of a test arc's great -circle, a crossing too close to an arc endpoint to call — and retries with -the next edge as anchor. If every anchor is degenerate the axis is -extended to `±1`, so the box can come out loose but never under-covers. +An extremum of a coordinate over a *region* on the sphere lies either on +the boundary or at an on-sphere critical point of that coordinate — across +the three coordinates, the six axis points `(±1,0,0)`, `(0,±1,0)`, +`(0,0,±1)`. A region's box is therefore the union of its edges' +[`UnitSpherical.spherical_arc_extent`](@ref)s plus an enclosure check per +axis point, decided by crossing parity as in `S2Loop::InitBound` +(`s2loop.cc`). + +Rings follow S2's loop convention (`s2loop.h`): CCW, interior on the left, +so a clockwise ring encloses the complement. Configurations too close to +degenerate for [`UnitSpherical.spherical_orient`](@ref) to call are retried +with the next edge as anchor; if every anchor fails, the axis is extended +to `±1`, so the box can come out loose but never under-covers. =# """ extent(m::Manifold, geom, [::Type{T} = Float64])::Extents.Extent -The extent of `geom` on the manifold `m` — this method lives on, and -returns an, `Extents.Extent`. +The extent of `geom` on the manifold `m`, as an `Extents.Extent`. The +method extends `Extents.extent` (GeometryOps does not export `extent`), so +call it as `GO.extent(m, geom)`. On `Planar()`, `GI.extent(geom)`. On `Spherical()`, the 3D Cartesian extent of the geometry on the unit sphere, with geographic (longitude, latitude) input converted like `UnitSphericalPoint`: curves are covered by -the union of their edges' great-circle arc extents, and rings and polygons -are treated as regions — wound CCW with the interior on the left, per S2's -loop convention — whose extent also covers any enclosed pole or other -on-sphere axis extreme. +the union of their edges' great-circle arc extents; rings and polygons are +regions — wound CCW with the interior on the left, per S2's loop +convention — whose extent also covers any enclosed axis point (a pole, +say). ## Example @@ -144,28 +130,25 @@ end # Crossing parity of the test arc q → m against ring edge a → b: 1 for a # transversal crossing, 0 for none, -1 for too close to degenerate to call. function _arc_crossing_parity(q, m, a, b) - # a vertex exactly antipodal to `q` lies on every great circle through - # `q`; its edges can reach the test arc only at `q` itself, excluded by - # the on-boundary check + # a vertex at `-q` lies on every great circle through `q`; its edges can + # reach the test arc only at `q` itself, excluded by the on-boundary check (a == -q || b == -q) && return 0 sa = UnitSpherical.spherical_orient(q, m, a) sb = UnitSpherical.spherical_orient(q, m, b) (sa == 0 || sb == 0) && return -1 sa == sb && return 0 - # `q` on this edge's great circle but off the edge (checked upfront): - # the circles meet only at `±q`, both out of the test arc's reach — no - # crossing. This degeneracy is anchor-independent (a lonlat grid's - # meridian edges hold `±eₓ`/`±e_y` exactly), so it resolves instead of - # returning -1. + # `q` on this edge's great circle but off the edge (checked upfront): the + # circles meet only at `±q`, out of the test arc's reach — no crossing. + # Anchor-independent (lonlat meridian edges hold `±eₓ`/`±e_y` exactly), + # so resolve instead of returning -1. sq = UnitSpherical.spherical_orient(a, b, q) sq == 0 && return 0 sm = UnitSpherical.spherical_orient(a, b, m) sm == 0 && return -1 sq == sm && return 0 # each arc now crosses the other's great circle exactly once, at one of - # the two antipodal circle intersections; the arcs cross iff those are - # the same point, i.e. iff the intersection direction `x` points into - # both arcs' hemispheres + # the two antipodal circle intersections; the arcs cross iff it is the + # same one, i.e. iff `x` points into both arcs' hemispheres x = cross(normalize(UnitSpherical.robust_cross_product(q, m)), normalize(UnitSpherical.robust_cross_product(a, b))) d1 = dot(x, q + m) diff --git a/src/utils/UnitSpherical/arc_extent.jl b/src/utils/UnitSpherical/arc_extent.jl index 9577f616f5..63302c5afe 100644 --- a/src/utils/UnitSpherical/arc_extent.jl +++ b/src/utils/UnitSpherical/arc_extent.jl @@ -5,22 +5,18 @@ spherical_arc_extent ``` -A great-circle arc bulges away from the chord between its endpoints, so the -endpoints' bounding box does not in general contain the arc: two points at -`z = 0.9` on either side of the prime meridian are joined by an arc whose -midpoint has `z = 0.9 / cos(θ/2) > 0.9`. [`spherical_arc_extent`](@ref) -computes a box that contains the whole arc. - -With `t̂ₐ` the unit tangent at `a` pointing along the arc, the arc is -`p(φ) = a cos(φ) + t̂ₐ sin(φ)` for `φ ∈ [0, θ]`, so each Cartesian -coordinate is a sinusoid `pᵢ(φ) = Rᵢ cos(φ - φᵢ)` with amplitude -`Rᵢ = hypot(aᵢ, t̂ₐᵢ)`. Since `θ ≤ π` there is at most one interior maximum -and one interior minimum per axis: a maximum exists iff `pᵢ` increases at -`a` and decreases at `b` (`t̂ₐᵢ > 0 > t̂ᵦᵢ`), where it attains `Rᵢ`; minima -mirror. The tangents come from [`robust_cross_product`](@ref), which keeps -nearly-degenerate and nearly-antipodal arcs stable. Bounds are padded by a -few ulps to absorb floating point error, as S2's `S2LatLngRectBounder` pads -by its maximum error. +A great-circle arc bulges away from the chord between its endpoints (two +points at `z = 0.9` joined over the pole reach `z = 1`), so the endpoints' +bounding box does not contain the arc. [`spherical_arc_extent`](@ref) +computes a box that does. + +Each Cartesian coordinate along the arc is the sinusoid +`pᵢ(φ) = aᵢ cos φ + t̂ₐᵢ sin φ` (`t̂ₐ` the unit tangent at `a`), of +amplitude `hypot(aᵢ, t̂ₐᵢ)`; the arc spans at most a half turn, so an +interior extremum on axis `i` exists iff `pᵢ` rises at `a` and falls at +`b`. Tangents come from [`robust_cross_product`](@ref), keeping +nearly-degenerate and nearly-antipodal arcs stable; bounds are padded by a +few ulps, as S2's `S2LatLngRectBounder` pads by its maximum error. =# """ diff --git a/test/methods/extent.jl b/test/methods/extent.jl index 8868ea2cd4..2ab76cfbaa 100644 --- a/test/methods/extent.jl +++ b/test/methods/extent.jl @@ -78,7 +78,7 @@ end @testset "Dumbbell: both poles through a thin corridor" begin # two polar caps (above ±85°) joined by a corridor over lon ∈ (355°, 5°): - # small area, zero net winding about every axis, and both poles and + # small area, zero net winding about the polar axis, and both poles and # (1, 0, 0) interior north = [(lon, 85.0) for lon in range(5.0, 355.0; length = 15)] # eastward, pole on the left south = [(lon, -85.0) for lon in range(355.0, 5.0; length = 15)] # westward, pole on the left From 84e02ceaf0c154cf3bfa8ee7a5afba96ee9c4d4c Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 11 Jul 2026 16:53:45 -0400 Subject: [PATCH 17/20] Delegate spherical relate bounds to the shared extent substrate `rk_interaction_bounds(::Spherical)` now composes `spherical_arc_extent` and `_spherical_region_extent` over kernel-converted points instead of carrying a private extent stack; `arc_extent`, `_point_box`, `_arcs_extent`, `_sph_bounds`, and `_widen_area_axes` are deleted. Area widening thereby gains the parity test's super-hemisphere and degeneracy handling. Rings keep dim-1 curve bounds (a CW hole must not become a complement region), and boxes keep a few ulps of padding. The antipodal-edge throw moves from the extent computation to explicit ingest validation (`_validate_relate_edges`, run in the extent-cache pass), since `spherical_arc_extent` deliberately picks a stable plane instead of throwing. Stored-extent reuse on `Spherical` now requires a 3D `(X, Y, Z)` extent (a user's lon/lat box must not be compared against unit-sphere boxes), and `_union_stored_extents` computes missing child extents through `rk_interaction_bounds` instead of `GI.extent`. Co-Authored-By: Claude Fable 5 --- .../relateng/edge_intersector.jl | 4 +- src/methods/geom_relations/relateng/kernel.jl | 6 + .../relateng/kernel_spherical.jl | 143 ++++++++---------- .../relateng/relate_geometry.jl | 30 ++-- .../correction/antipodal_edge_split.jl | 7 +- .../relateng/kernel_conformance_spherical.jl | 18 +-- test/methods/relateng/spherical_end_to_end.jl | 8 +- 7 files changed, 102 insertions(+), 114 deletions(-) diff --git a/src/methods/geom_relations/relateng/edge_intersector.jl b/src/methods/geom_relations/relateng/edge_intersector.jl index 673f2b0157..03ed7a6441 100644 --- a/src/methods/geom_relations/relateng/edge_intersector.jl +++ b/src/methods/geom_relations/relateng/edge_intersector.jl @@ -287,7 +287,7 @@ _segment_envs_disjoint(::Planar, a0, a1, b0, b1) = # extent, which is exact in xyz and has no antimeridian pathology — so this # prune is valid (and tighter than no pruning). _segment_envs_disjoint(::Spherical, a0, a1, b0, b1) = - !Extents.intersects(arc_extent(a0, a1), arc_extent(b0, b1)) + !Extents.intersects(spherical_arc_extent(a0, a1), spherical_arc_extent(b0, b1)) _segment_envs_disjoint(::Manifold, a0, a1, b0, b1) = false # The spatial index built over per-segment extents for the tree-accelerated @@ -427,7 +427,7 @@ end # the endpoints, or the bulge-aware 3D great-circle arc extent on the sphere # (the endpoint box would miss a long arc's bulge and wrongly prune crossings). _segment_extent(::Planar, p, q) = Extents.Extent(X = minmax(p[1], q[1]), Y = minmax(p[2], q[2])) -_segment_extent(::Spherical, p, q) = arc_extent(p, q) +_segment_extent(::Spherical, p, q) = spherical_arc_extent(p, q) _segment_extent_type(::Planar) = Extents.Extent{(:X, :Y), NTuple{2, NTuple{2, Float64}}} _segment_extent_type(::Spherical) = Extents.Extent{(:X, :Y, :Z), NTuple{3, NTuple{2, Float64}}} diff --git a/src/methods/geom_relations/relateng/kernel.jl b/src/methods/geom_relations/relateng/kernel.jl index a6cf3ab01f..f5f77f4da9 100644 --- a/src/methods/geom_relations/relateng/kernel.jl +++ b/src/methods/geom_relations/relateng/kernel.jl @@ -261,6 +261,12 @@ function _to_kernel_points(m::Manifold, geom) return pts end +# Manifold hook for edge validation at ingest (the `RelateGeometry` +# extent-cache pass): a manifold may reject edges it cannot represent. Planar +# edges are always fine; `Spherical` throws on exactly-antipodal vertex pairs +# (see kernel_spherical.jl). +_validate_relate_edges(::Manifold, curve) = nothing + # Degenerate interaction box of a single *kernel* point, matching the # dimensionality of `rk_interaction_bounds` (2D for planar tuples, 3D for 3D # kernel points such as `UnitSphericalPoint`). Used by the point-locator line diff --git a/src/methods/geom_relations/relateng/kernel_spherical.jl b/src/methods/geom_relations/relateng/kernel_spherical.jl index 1f71c96a1a..86c5082d2f 100644 --- a/src/methods/geom_relations/relateng/kernel_spherical.jl +++ b/src/methods/geom_relations/relateng/kernel_spherical.jl @@ -91,12 +91,6 @@ end _kernel_point_type(::Spherical) = UnitSphericalPoint{Float64} @inline _to_kernel_point(::Spherical, p) = _spherical_kernel_point(p) -# 3D AABB of a minor great-circle arc (spike-proven, 0/102k fuzz escapes). A -# great-circle arc bulges outside the coordinate box of its endpoints; the -# extremum of coordinate i on the circle with normal n is at ±w, the normalized -# projection of axis eᵢ onto the circle's plane. The box is extended by whichever -# of ±w lies on the minor arc; a few ulps of widening absorb the roundoff in w. -@inline _on_minor_arc(w, a, b, n) = (cross(a, w) ⋅ n) >= 0.0 && (cross(w, b) ⋅ n) >= 0.0 @inline _widen(lo, hi) = (prevfloat(lo, 4), nextfloat(hi, 4)) @noinline _throw_antipodal_edge(a, b) = throw(ArgumentError( @@ -104,42 +98,27 @@ _kernel_point_type(::Spherical) = UnitSphericalPoint{Float64} "unique great-circle arc; densify it first with the `AntipodalEdgeSplit` " * "correction (it inserts the lon/lat midpoint)")) -function arc_extent(a, b) - n = cross(a, b) - n2 = n ⋅ n - #-- a vanishing normal with the endpoints pointing opposite ways means the - #-- vertices are exactly antipodal: infinitely many great circles pass - #-- through them, so the edge has no well-defined arc. (A vanishing normal - #-- with a·b > 0 is a zero-length/repeated vertex, which is fine.) - n2 == 0.0 && (a ⋅ b) < 0.0 && _throw_antipodal_edge(a, b) - xlo, xhi = minmax(a[1], b[1]); ylo, yhi = minmax(a[2], b[2]); zlo, zhi = minmax(a[3], b[3]) - if n2 > 0.0 - invn2 = inv(n2) - @inbounds for i in 1:3 - ei_n = n[i] - wx = (i == 1) - ei_n * n[1] * invn2 - wy = (i == 2) - ei_n * n[2] * invn2 - wz = (i == 3) - ei_n * n[3] * invn2 - wnorm = sqrt(wx^2 + wy^2 + wz^2) - wnorm == 0.0 && continue # axis ⟂ plane: coordinate constant 0, endpoints cover it - w = UnitSphericalPoint(wx / wnorm, wy / wnorm, wz / wnorm) - for ww in (w, -w) - if _on_minor_arc(ww, a, b, n) - ci = ww[i] - if i == 1; xlo = min(xlo, ci); xhi = max(xhi, ci) - elseif i == 2; ylo = min(ylo, ci); yhi = max(yhi, ci) - else; zlo = min(zlo, ci); zhi = max(zhi, ci) - end - end - end - end +# Ingest validation, run once per curve at `RelateGeometry` construction (the +# extent-cache pass). A vanishing cross product with the endpoints pointing +# opposite ways means the vertices are exactly antipodal: infinitely many +# great circles pass through them, so the edge has no well-defined arc and +# relate's contract is to throw informatively. (`spherical_arc_extent` +# deliberately does not throw — it picks a stable plane — so the guard lives +# here, not in the extent.) A vanishing cross with `a ⋅ b > 0` is a +# zero-length/repeated vertex, which is fine. +function _validate_relate_edges(::Spherical, curve) + n = GI.npoint(curve) + n < 2 && return nothing + prev = _spherical_kernel_point(GI.getpoint(curve, 1)) + for i in 2:n + cur = _spherical_kernel_point(GI.getpoint(curve, i)) + c = cross(prev, cur) + (c ⋅ c) == 0.0 && (prev ⋅ cur) < 0.0 && _throw_antipodal_edge(prev, cur) + prev = cur end - return Extents.Extent(X = _widen(xlo, xhi), Y = _widen(ylo, yhi), Z = _widen(zlo, zhi)) + return nothing end -@inline _point_box(u) = - Extents.Extent(X = _widen(GI.x(u), GI.x(u)), Y = _widen(GI.y(u), GI.y(u)), Z = _widen(GI.z(u), GI.z(u))) - # ## rk_classify_intersection # # The two great circles meet at ±d, d = (a0×a1)×(b0×b1). `SS_PROPER` iff one of @@ -351,57 +330,59 @@ function rk_point_in_ring(m::Spherical, p, ring; exact) return (isodd(crossings) ⊻ pole_inside) ? LOC_INTERIOR : LOC_EXTERIOR end -# Interaction bounds on the sphere: a 3D `Extent{(:X,:Y,:Z)}` in unit-sphere xyz -# (the engine works in xyz after ingest), as the union of `arc_extent` over the -# geometry's edges, plus — for area elements — an axis-point extension so the box -# covers the interior, not just the boundary slab. -rk_interaction_bounds(m::Spherical, geom) = _sph_bounds(m, GI.trait(geom), geom) - -# Converted (lon/lat → unit xyz) vertices of a ring/curve. -_ring_usp(ring) = [_spherical_kernel_point(p) for p in GI.getpoint(ring)] - -# Union of arc_extents over consecutive vertices (a single point box if n == 1). -function _arcs_extent(usp) - length(usp) == 1 && return _point_box(usp[1]) - ext = arc_extent(usp[1], usp[2]) - for i in 2:length(usp)-1 - ext = Extents.union(ext, arc_extent(usp[i], usp[i+1])) +# Interaction bounds on the sphere, built from the shared substrate +# (`spherical_arc_extent` per edge, `_spherical_region_extent` for area +# interiors) over kernel-converted points, so the box and the ingested +# vertices agree bit-for-bit. Relate-specific glue: rings are a polygon's +# linework / dim-1 curve elements (JTS semantics), not S2 regions — their +# edges are bounded directly, so a CW hole cannot become a complement region — +# and every box gets a few ulps of padding so a kernel point that differs +# from another conversion path by sub-ulp noise still prunes as interacting. +rk_interaction_bounds(m::Spherical, geom) = + _pad_bounds(_sph_interaction_extent(m, GI.trait(geom), geom)) + +function _sph_interaction_extent(m::Spherical, ::GI.AbstractPointTrait, geom) + p = _spherical_kernel_point(geom) + return Extents.Extent(X = (p[1], p[1]), Y = (p[2], p[2]), Z = (p[3], p[3])) +end +function _sph_interaction_extent(m::Spherical, ::GI.AbstractCurveTrait, geom) + n = GI.npoint(geom) + prev = _spherical_kernel_point(GI.getpoint(geom, 1)) + ext = spherical_arc_extent(prev, prev) + for i in 2:n + cur = _spherical_kernel_point(GI.getpoint(geom, i)) + ext = Extents.union(ext, spherical_arc_extent(prev, cur)) + prev = cur end return ext end - -function _sph_bounds(::Spherical, ::GI.AbstractPointTrait, geom) - return _point_box(_spherical_kernel_point(geom)) -end -_sph_bounds(::Spherical, ::GI.AbstractCurveTrait, geom) = _arcs_extent(_ring_usp(geom)) -function _sph_bounds(::Spherical, ::GI.AbstractPolygonTrait, geom) - exterior = _ring_usp(GI.getexterior(geom)) - ext = _arcs_extent(exterior) +function _sph_interaction_extent(m::Spherical, ::GI.AbstractPolygonTrait, geom) + # region box of the exterior ring: edge arc extents plus enclosed-axis + # widening, on the same converted points the engine ingests + ext = _spherical_region_extent(_ring_usp(GI.getexterior(geom))) + # a valid polygon's holes lie inside that region — but JTS's element + # envelope also covers a stray hole outside the shell, and extraction + # relies on that to keep the element alive + # (see `_extract_segment_strings_from_atomic!`) for hole in GI.gethole(geom) - ext = Extents.union(ext, _arcs_extent(_ring_usp(hole))) + GI.isempty(hole) && continue + ext = Extents.union(ext, _sph_interaction_extent(m, GI.trait(hole), hole)) end - # An area interior reaches beyond its boundary slab (e.g. a ring around a - # pole has a thin boundary band but its interior reaches z = 1). Widen each - # axis whose ±eᵢ is interior to the exterior ring out to ±1 (conservative — - # over-covering can only under-prune, never miss an interaction). - return _widen_area_axes(ext, exterior) -end - -function _widen_area_axes(ext, pts) - xlo, xhi = ext.X; ylo, yhi = ext.Y; zlo, zhi = ext.Z - _ring_contains_dir(pts, UnitSphericalPoint(1.0, 0.0, 0.0)) && (xhi = 1.0) - _ring_contains_dir(pts, UnitSphericalPoint(-1.0, 0.0, 0.0)) && (xlo = -1.0) - _ring_contains_dir(pts, UnitSphericalPoint(0.0, 1.0, 0.0)) && (yhi = 1.0) - _ring_contains_dir(pts, UnitSphericalPoint(0.0, -1.0, 0.0)) && (ylo = -1.0) - _ring_contains_dir(pts, UnitSphericalPoint(0.0, 0.0, 1.0)) && (zhi = 1.0) - _ring_contains_dir(pts, UnitSphericalPoint(0.0, 0.0, -1.0)) && (zlo = -1.0) - return Extents.Extent(X = (xlo, xhi), Y = (ylo, yhi), Z = (zlo, zhi)) + return ext end -function _sph_bounds(m::Spherical, ::GI.AbstractGeometryTrait, geom) +function _sph_interaction_extent(m::Spherical, ::GI.AbstractGeometryTrait, geom) ext = nothing for g in GI.getgeom(geom) - e = _sph_bounds(m, GI.trait(g), g) + GI.isempty(g) && continue + e = _sph_interaction_extent(m, GI.trait(g), g) ext = ext === nothing ? e : Extents.union(ext, e) end return ext end + +# Converted (kernel-ingest: unit, signed-zero) vertices of a ring/curve. +_ring_usp(ring) = [_spherical_kernel_point(p) for p in GI.getpoint(ring)] + +_pad_bounds(::Nothing) = nothing +_pad_bounds(ext) = Extents.Extent( + X = _widen(ext.X...), Y = _widen(ext.Y...), Z = _widen(ext.Z...)) diff --git a/src/methods/geom_relations/relateng/relate_geometry.jl b/src/methods/geom_relations/relateng/relate_geometry.jl index ce5577f1cb..b6541e8b82 100644 --- a/src/methods/geom_relations/relateng/relate_geometry.jl +++ b/src/methods/geom_relations/relateng/relate_geometry.jl @@ -133,25 +133,36 @@ _has_stored_extent(geom) = geom isa GI.Wrappers.WrapperGeometry && hasproperty(geom, :extent) && geom.extent isa Extents.Extent +# A stored extent is reusable as interaction bounds only if it lives in the +# space the kernel prunes in: any stored extent on `Planar`, but only a 3D +# `(X, Y, Z)` extent on `Spherical` — user inputs typically carry lon/lat +# boxes, which must not be compared against unit-sphere boxes. Our own cache +# pass always stores `(X, Y, Z)`; a user storing one is trusted to mean it. +_reusable_stored_extent(::Manifold, geom) = _has_stored_extent(geom) +_reusable_stored_extent(::Spherical, geom) = + _has_stored_extent(geom) && geom.extent isa Extents.Extent{(:X, :Y, :Z)} + _relate_cache_extents(m::Manifold, geom) = _relate_cache_extents(m, GI.trait(geom), geom) #-- point elements: their extent is themselves, nothing to cache _relate_cache_extents(::Manifold, ::Union{GI.AbstractPointTrait, GI.AbstractMultiPointTrait}, geom) = geom -#-- linework leaves: lines and rings (the only level where coordinates are read) +#-- linework leaves: lines and rings (the only level where coordinates are +#-- read, and hence where the manifold's edge validation runs) function _relate_cache_extents(m::Manifold, trait::GI.AbstractCurveTrait, line) - (GI.isempty(line) || _has_stored_extent(line)) && return line + (GI.isempty(line) || _reusable_stored_extent(m, line)) && return line + _validate_relate_edges(m, line) return GI.geointerface_geomtype(trait)(line; extent = rk_interaction_bounds(m, line), crs = GI.crs(line)) end function _relate_cache_extents(m::Manifold, trait::GI.AbstractPolygonTrait, poly) GI.isempty(poly) && return poly - if _has_stored_extent(poly) && all(r -> GI.isempty(r) || _has_stored_extent(r), GI.getring(poly)) + if _reusable_stored_extent(m, poly) && all(r -> GI.isempty(r) || _reusable_stored_extent(m, r), GI.getring(poly)) return poly end rings = [_relate_cache_extents(m, GI.trait(r), r) for r in GI.getring(poly)] - ext = _union_stored_extents(rings) + ext = _union_stored_extents(m, rings) ext === nothing && return poly return GI.geointerface_geomtype(trait)(rings; extent = ext, crs = GI.crs(poly)) end @@ -159,7 +170,7 @@ end #-- collections (covers Multi* types too): recurse, union the child extents function _relate_cache_extents(m::Manifold, trait::GI.AbstractGeometryCollectionTrait, geom) children = [_relate_cache_extents(m, GI.trait(g), g) for g in GI.getgeom(geom)] - ext = _union_stored_extents(children) + ext = _union_stored_extents(m, children) ext === nothing && return geom return GI.geointerface_geomtype(trait)(children; extent = ext, crs = GI.crs(geom)) end @@ -169,16 +180,17 @@ _relate_cache_extents(::Manifold, ::GI.AbstractTrait, geom) = geom # Union of the children's extents, reading stored ones and computing only # for non-empty children that have none (e.g. point members of a GC); -# `nothing` when no child contributes one. -function _union_stored_extents(children) +# `nothing` when no child contributes one. Computed extents go through +# `rk_interaction_bounds` so they stay in the manifold's coordinate space. +function _union_stored_extents(m::Manifold, children) ext = nothing for c in children - ce = if _has_stored_extent(c) + ce = if _reusable_stored_extent(m, c) c.extent elseif GI.isempty(c) nothing else - GI.extent(c; fallback = true) + rk_interaction_bounds(m, c) end ce === nothing && continue ext = ext === nothing ? ce : Extents.union(ext, ce) diff --git a/src/transformations/correction/antipodal_edge_split.jl b/src/transformations/correction/antipodal_edge_split.jl index 1ca9f11737..2465dd80fd 100644 --- a/src/transformations/correction/antipodal_edge_split.jl +++ b/src/transformations/correction/antipodal_edge_split.jl @@ -5,8 +5,8 @@ export AntipodalEdgeSplit #= On the sphere an edge between two *antipodal* vertices (a point and its antipode) has no unique great-circle arc — infinitely many great circles pass -through both — so the spherical RelateNG kernel refuses such an edge (see -[`arc_extent`](@ref) / `relate` with a `Spherical` manifold). This correction +through both — so `relate` with a `Spherical` manifold refuses such an edge +at ingest (the kernel's edge validation). This correction is the documented remedy: it splits every antipodal edge by inserting the lon/lat midpoint of its endpoints, replacing one ambiguous edge with two well-defined (roughly quarter-circle) arcs. @@ -52,7 +52,8 @@ end (::AntipodalEdgeSplit)(::GI.AbstractCurveTrait, curve) = _split_antipodal_curve(curve) # Whether the lon/lat points `p`, `q` map to exactly-antipodal unit vectors — -# the same condition `arc_extent` throws on (vanishing cross, negative dot). +# the same condition the kernel's edge validation throws on (vanishing cross, +# negative dot). function _is_antipodal_lonlat(p, q) up = _spherical_kernel_point(p) uq = _spherical_kernel_point(q) diff --git a/test/methods/relateng/kernel_conformance_spherical.jl b/test/methods/relateng/kernel_conformance_spherical.jl index 7547405959..abb582ffc0 100644 --- a/test/methods/relateng/kernel_conformance_spherical.jl +++ b/test/methods/relateng/kernel_conformance_spherical.jl @@ -53,20 +53,8 @@ function kernel_conformance_suite_spherical(m; exact) @test !GO.rk_point_on_segment(m, _usp(0, 0, 1), a, b; exact) # pole, off the circle @test !GO.rk_point_on_segment(m, _usp(-1, 1, 0), a, b; exact) # on circle, outside the minor-arc span end - @testset "arc_extent contains the arc (bulge captured)" begin - rng2 = Random.MersenneTwister(42) - for _ in 1:500 - p = GO.rk_normalize_usp(_usp(randn(rng2), randn(rng2), randn(rng2))) - q = GO.rk_normalize_usp(_usp(randn(rng2), randn(rng2), randn(rng2))) - e = GO.arc_extent(p, q) - for s in 0:0.05:1 - u = slerp(p, q, s) - @test e.X[1] <= GI.x(u) <= e.X[2] - @test e.Y[1] <= GI.y(u) <= e.Y[2] - @test e.Z[1] <= GI.z(u) <= e.Z[2] - end - end - end + # arc containment (bulge capture) is the shared `spherical_arc_extent`'s + # contract, tested exhaustively in test/utils/unitspherical.jl @testset "rk_interaction_bounds is 3D and contains the converted vertices" begin ring = GI.LinearRing([(0.,0.), (10.,0.), (10.,10.), (0.,10.), (0.,0.)]) @@ -274,7 +262,7 @@ function kernel_conformance_suite_spherical(m; exact) verts = [_usp(2,0,1), _usp(0,2,1), _usp(-2,0,1), _usp(0,-2,1), _usp(2,0,1)] poly = GI.Polygon([GI.LinearRing(verts)]) e = GO.rk_interaction_bounds(m, poly) - @test e.Z[2] == 1.0 # interior reaches the enclosed north pole + @test e.Z[2] >= 1.0 # interior reaches the enclosed north pole # the boundary ring (curve) bound tops out well below the pole eb = GO.rk_interaction_bounds(m, GI.LinearRing(verts)) @test eb.Z[2] < 0.99 diff --git a/test/methods/relateng/spherical_end_to_end.jl b/test/methods/relateng/spherical_end_to_end.jl index 3d16297d35..6e59edc61f 100644 --- a/test/methods/relateng/spherical_end_to_end.jl +++ b/test/methods/relateng/spherical_end_to_end.jl @@ -32,7 +32,7 @@ end # 2D coordinate box. A long near-equatorial arc (lon 0 → 170) bulges to y ≈ 1 # at lon 90 while its endpoint box has y ∈ [0, 0.17]; a polygon straddling the # equator at lon 90 crosses it there. A 2D endpoint box prunes that pair away -# (wrong DE-9IM); the bulge-aware `arc_extent` keeps it. +# (wrong DE-9IM); the bulge-aware `spherical_arc_extent` keeps it. @testset "spherical tree accelerator agrees with NestedLoop (arc bulge)" begin A = GI.Polygon([GI.LinearRing([(0., 0.), (170., 0.), (85., 40.), (0., 0.)])]) B = GI.Polygon([GI.LinearRing([(88., -2.), (92., -2.), (92., 2.), (88., 2.), (88., -2.)])]) @@ -78,12 +78,12 @@ end p0 = GO._to_kernel_point(Spherical(), (0., 0.)) # (1, 0, 0) p180 = GO._to_kernel_point(Spherical(), (180., 0.)) # (-1, 0, 0): antipodal p90 = GO._to_kernel_point(Spherical(), (90., 0.)) # (0, 1, 0) - err = try; GO.arc_extent(p0, p180); nothing; catch e; e; end + err = try; GO._validate_relate_edges(Spherical(), GI.LineString([p0, p180])); nothing; catch e; e; end @test err isa ArgumentError @test occursin("AntipodalEdgeSplit", err.msg) #-- a normal edge and a repeated (zero-length) vertex do NOT throw - @test (GO.arc_extent(p0, p90); true) - @test (GO.arc_extent(p0, p0); true) + @test (GO._validate_relate_edges(Spherical(), GI.LineString([p0, p90])); true) + @test (GO._validate_relate_edges(Spherical(), GI.LineString([p0, p0])); true) #-- the whole relate rejects a polygon carrying an antipodal edge at ingest bad = GI.Polygon([GI.LinearRing([(0., 0.), (180., 0.), (90., 80.), (0., 0.)])]) @test_throws ArgumentError GO.relate(alg, bad, GI.Point(10., 10.)) From ea69854c3c57e0375a2343146b926948e4eab3d6 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 11 Jul 2026 17:01:36 -0400 Subject: [PATCH 18/20] Promote ring-containment parity to UnitSpherical, graft into `rk_point_in_ring` `spherical_ring_contains` moves from the extent internals to UnitSpherical as documented public API (the extent code and the relate kernel are its two consumers), with its three geometric predicates injectable: `orient`, `on_arc`, and `proper_crossing` default to the tolerance-based substrate predicates, so the extent path is unchanged. `rk_point_in_ring(::Spherical)` keeps its exact boundary classification and delegates the interior decision to the shared anchor-retry parity, injecting `rk_orient` and `_arcs_cross_properly` so the answer is as exact as the kernel's predicates. This replaces the fixed north-pole reference walk, whose documented degeneracies (query at/antipodal to the pole, vertex on the query meridian) and super-hemisphere winding limitation no longer apply; an all-anchors-degenerate configuration now throws instead of answering wrong. The scale-dependent guards in the shared function are normalized so the conformance suite's exact integer (non-unit) rings work unchanged. Co-Authored-By: Claude Fable 5 --- src/methods/extent.jl | 69 +---------- .../relateng/kernel_spherical.jl | 61 ++++------ src/utils/UnitSpherical/UnitSpherical.jl | 1 + src/utils/UnitSpherical/predicates.jl | 110 ++++++++++++++++++ 4 files changed, 136 insertions(+), 105 deletions(-) diff --git a/src/methods/extent.jl b/src/methods/extent.jl index 47b35c092c..5e07c677df 100644 --- a/src/methods/extent.jl +++ b/src/methods/extent.jl @@ -82,7 +82,7 @@ function _spherical_region_extent(pts::Vector{<:UnitSpherical.UnitSphericalPoint hi = MVector(ext.X[2], ext.Y[2], ext.Z[2]) for i in 1:3, s in (1.0, -1.0) q = UnitSpherical.UnitSphericalPoint(ntuple(j -> j == i ? s : 0.0, 3)) - inside = _spherical_ring_contains(pts, n, q) + inside = UnitSpherical.spherical_ring_contains(pts, n, q) # nothing (undecidable) extends too; the box must never under-cover if inside === nothing || inside s > 0 ? (hi[i] = one(hi[i])) : (lo[i] = -one(lo[i])) @@ -90,70 +90,3 @@ function _spherical_region_extent(pts::Vector{<:UnitSpherical.UnitSphericalPoint end return Extents.Extent(X = (lo[1], hi[1]), Y = (lo[2], hi[2]), Z = (lo[3], hi[3])) end - -# Crossing-parity containment of `q` in the closed region left of the ring, -# after S2Loop::Contains/InitBound. Returns `nothing` when every anchor edge -# is degenerate with respect to `q`. -function _spherical_ring_contains(pts, n, q) - for j in 1:n - UnitSpherical.point_on_spherical_arc(q, pts[j], pts[mod1(j + 1, n)]) && return true - end - for j in 1:n - a, b = pts[j], pts[mod1(j + 1, n)] - a == b && continue - side = UnitSpherical.spherical_orient(a, b, q) - side == 0 && continue - mid = a + b - norm(mid) < 1e-9 && continue # near-antipodal edge, midpoint unstable - m = UnitSpherical.UnitSphericalPoint(normalize(mid)) - dot(q, m) < -1 + 1e-9 && continue # test arc q → m would span a half turn - crossings = 0 - ok = true - for k in 1:n - k == j && continue - c = _arc_crossing_parity(q, m, pts[k], pts[mod1(k + 1, n)]) - if c == -1 - ok = false - break - end - crossings += c - end - ok || continue - # walking from `m` toward `q` departs onto `q`'s side of the anchor - # edge (the arc meets that great circle again only at `-m`); positive - # side is the interior, and each crossing flips it - return isodd(crossings) ? side < 0 : side > 0 - end - return nothing -end - -# Crossing parity of the test arc q → m against ring edge a → b: 1 for a -# transversal crossing, 0 for none, -1 for too close to degenerate to call. -function _arc_crossing_parity(q, m, a, b) - # a vertex at `-q` lies on every great circle through `q`; its edges can - # reach the test arc only at `q` itself, excluded by the on-boundary check - (a == -q || b == -q) && return 0 - sa = UnitSpherical.spherical_orient(q, m, a) - sb = UnitSpherical.spherical_orient(q, m, b) - (sa == 0 || sb == 0) && return -1 - sa == sb && return 0 - # `q` on this edge's great circle but off the edge (checked upfront): the - # circles meet only at `±q`, out of the test arc's reach — no crossing. - # Anchor-independent (lonlat meridian edges hold `±eₓ`/`±e_y` exactly), - # so resolve instead of returning -1. - sq = UnitSpherical.spherical_orient(a, b, q) - sq == 0 && return 0 - sm = UnitSpherical.spherical_orient(a, b, m) - sm == 0 && return -1 - sq == sm && return 0 - # each arc now crosses the other's great circle exactly once, at one of - # the two antipodal circle intersections; the arcs cross iff it is the - # same one, i.e. iff `x` points into both arcs' hemispheres - x = cross(normalize(UnitSpherical.robust_cross_product(q, m)), - normalize(UnitSpherical.robust_cross_product(a, b))) - d1 = dot(x, q + m) - d2 = dot(x, a + b) - tol = 16 * eps(Float64) * norm(x) - (abs(d1) <= tol || abs(d2) <= tol) && return -1 - return (d1 > 0) == (d2 > 0) ? 1 : 0 -end diff --git a/src/methods/geom_relations/relateng/kernel_spherical.jl b/src/methods/geom_relations/relateng/kernel_spherical.jl index 86c5082d2f..2ae7ebf1fd 100644 --- a/src/methods/geom_relations/relateng/kernel_spherical.jl +++ b/src/methods/geom_relations/relateng/kernel_spherical.jl @@ -262,9 +262,7 @@ function rk_nodes_coincide(::Spherical, k1::NodeKey, k2::NodeKey; exact) return _iszero3(_cross3(d1, d2)) && _dot3(d1, d2) > 0 end -# ## rk_point_in_ring (meridian-crossing parity, S2 convention) - -const _NPOLE = UnitSphericalPoint(0.0, 0.0, 1.0) +# ## rk_point_in_ring (anchor-retry crossing parity, S2 convention) # Whether the two minor arcs (p0,p1) and (q0,q1) cross properly (interior to # both). The great circles meet at ±d, d = (p0×p1)×(q0×q1); a proper crossing is @@ -279,34 +277,6 @@ function _arcs_cross_properly(bt, p0, p1, q0, q1) return _strictly_in_arc3(nd, P0, P1, na) && _strictly_in_arc3(nd, Q0, Q1, nb) end -# Whether direction `N` lies in the ring's interior (S2 interior-on-left), -# decided by the ring's winding around `N`: project each edge onto the plane ⊥ N -# and sum the signed turn (each < π in magnitude, so there is no atan2 branch -# ambiguity — unlike a per-triangle solid-angle sum). A single interior-on-left -# enclosure sweeps +2π. Robust for rings smaller than a hemisphere (the common -# geographic case); a super-hemisphere interior would under-report a -# non-encircled interior axis — a documented limitation, see the design doc. -function _ring_contains_dir(pts, N) - θ = 0.0 - for i in 1:length(pts)-1 - a = pts[i]; b = pts[i+1] - ua = a - (a ⋅ N) * N - ub = b - (b ⋅ N) * N - θ += atan(N ⋅ cross(ua, ub), ua ⋅ ub) - end - return θ > π -end -@inline _ring_contains_pole(pts) = _ring_contains_dir(pts, _NPOLE) - -# Location of `p` relative to the area enclosed by `ring`. Boundary first (exact -# arc membership), then parity of proper crossings of the minor arc p→NPOLE with -# the ring edges — "same region as the pole" — resolved to interior/exterior by -# the pole's own containment. -# -# NOTE: the north-pole reference is degenerate when `p` is at/antipodal to the -# pole, or a ring vertex sits exactly on the p→pole meridian. The engine inputs -# in this work avoid those; the deterministic-perturbation / other-pole -# treatment (mirroring planar RayCrossingCounter) is a follow-up. # Ring vertices as spherical kernel points. A 3D ring is already in kernel # coordinates (e.g. the conformance suite's exact integer USP rings) and is read # verbatim — renormalizing would perturb the exact orient the boundary test @@ -316,20 +286,37 @@ _ring_kernel_pts(ring) = _ring_kernel_pts(booltype(GI.is3d(GI.getpoint(ring, 1)) _ring_kernel_pts(::True, ring) = _node_points(ring) _ring_kernel_pts(::False, ring) = _ring_usp(ring) +# Location of `p` relative to the area enclosed by `ring` (S2 convention: +# interior on the left). Boundary first (exact arc membership), then the +# shared anchor-retry crossing parity — `spherical_ring_contains` with this +# kernel's predicates injected: `rk_orient` for sides and +# `_arcs_cross_properly` for transversality, so the interior decision is as +# exact as the predicates. Unlike a fixed reference point (the previous +# north-pole walk), the anchor retry has no preferred direction to go +# degenerate at, and handles super-hemisphere interiors. If every anchor is +# degenerate with respect to `p` — unreachable for a non-degenerate ring and +# an off-boundary point — the configuration is refused, not answered wrong. function rk_point_in_ring(m::Spherical, p, ring; exact) pts = _ring_kernel_pts(ring) @inbounds for i in 1:length(pts)-1 rk_point_on_segment(m, p, pts[i], pts[i+1]; exact) && return LOC_BOUNDARY end bt = booltype(exact) - crossings = 0 - @inbounds for i in 1:length(pts)-1 - _arcs_cross_properly(bt, p, _NPOLE, pts[i], pts[i+1]) && (crossings += 1) - end - pole_inside = _ring_contains_pole(pts) - return (isodd(crossings) ⊻ pole_inside) ? LOC_INTERIOR : LOC_EXTERIOR + n = length(pts) + n > 1 && pts[end] == pts[1] && (n -= 1) + inside = spherical_ring_contains(pts, n, p; + orient = (a, b, c) -> rk_orient(m, a, b, c; exact), + on_arc = Returns(false), # boundary classified exactly above + proper_crossing = (q, mid, a, b) -> _arcs_cross_properly(bt, q, mid, a, b) ? 1 : 0) + inside === nothing && _throw_degenerate_point_in_ring(p) + return inside ? LOC_INTERIOR : LOC_EXTERIOR end +@noinline _throw_degenerate_point_in_ring(p) = throw(ArgumentError( + "rk_point_in_ring: every anchor edge of the ring is degenerate with " * + "respect to the query point $(_tup3(p)) — the ring is degenerate at " * + "this point")) + # Interaction bounds on the sphere, built from the shared substrate # (`spherical_arc_extent` per edge, `_spherical_region_extent` for area # interiors) over kernel-converted points, so the box and the ingested diff --git a/src/utils/UnitSpherical/UnitSpherical.jl b/src/utils/UnitSpherical/UnitSpherical.jl index eba53275c2..fd92b3e498 100644 --- a/src/utils/UnitSpherical/UnitSpherical.jl +++ b/src/utils/UnitSpherical/UnitSpherical.jl @@ -25,6 +25,7 @@ include("arc_extent.jl") export UnitSphericalPoint, UnitSphereFromGeographic, GeographicFromUnitSphere, slerp, SphericalCap, spherical_distance, spherical_orient, point_on_spherical_arc, + spherical_ring_contains, spherical_arc_intersection, ArcIntersectionResult, arc_cross, arc_hinge, arc_overlap, arc_disjoint, spherical_arc_extent, diff --git a/src/utils/UnitSpherical/predicates.jl b/src/utils/UnitSpherical/predicates.jl index e0e9d36d74..1449bde36e 100644 --- a/src/utils/UnitSpherical/predicates.jl +++ b/src/utils/UnitSpherical/predicates.jl @@ -92,3 +92,113 @@ function point_on_spherical_arc(p::UnitSphericalPoint, a::UnitSphericalPoint, b: # (in terms of angle, so larger dot product) return (ap ≥ ab - tol) && (bp ≥ ab - tol) end + +""" + spherical_ring_contains(pts, n, q; orient, on_arc, proper_crossing) -> Union{Bool, Nothing} + +Whether `q` lies in the closed region on the left of the ring `pts[1:n]` +(S2 loop convention: counterclockwise winding, interior on the left — a +clockwise ring therefore contains the complement of the region it visually +encloses). The closing edge `pts[n] → pts[1]` is implied; points on the +boundary count as contained. Returns `nothing` when every anchor edge is +degenerate with respect to `q` — callers must treat that conservatively. + +Containment is decided by crossing parity, the way `S2Loop::Contains` / +`InitBound` decide pole containment: which side of an anchor edge `q` falls +on, flipped once per transversal crossing of the arc from the anchor's +midpoint to `q` with the other edges. Anchors degenerate with respect to +`q` are skipped and the next edge tried. + +The geometric predicates are injectable, so callers with stricter +requirements can substitute their own: + +- `orient(a, b, c)`: sign-valued orientation of `c` against the oriented + great circle through `a, b`; default [`spherical_orient`](@ref). +- `on_arc(q, a, b)::Bool`: boundary membership; default + [`point_on_spherical_arc`](@ref). Pass `Returns(false)` when boundary + points have already been classified. +- `proper_crossing(q, m, a, b)::Int`: `1` if the minor arcs `(q, m)` and + `(a, b)` cross transversally in both interiors, `0` if not, `-1` for too + close to degenerate to call; only consulted once `orient` places both + endpoint pairs strictly transversally. The default decides through + `robust_cross_product` with a small tolerance band and assumes unit + input. + +The injected predicates receive the input points untouched (which may be +non-unit for scale-invariant predicates); only the constructed reference +midpoint is normalized. +""" +function spherical_ring_contains(pts, n, q; + orient = spherical_orient, + on_arc = point_on_spherical_arc, + proper_crossing = _hemisphere_proper_crossing) + for j in 1:n + on_arc(q, pts[j], pts[mod1(j + 1, n)]) && return true + end + nq = norm(q) + for j in 1:n + a, b = pts[j], pts[mod1(j + 1, n)] + a == b && continue + side = orient(a, b, q) + side == 0 && continue + mid = a + b + # near-antipodal edge: the midpoint direction is unstable + norm(mid) < 1e-9 * (norm(a) + norm(b)) && continue + m = UnitSphericalPoint(normalize(mid)) + # test arc q → m would span a half turn + dot(q, m) < (-1 + 1e-9) * nq && continue + crossings = 0 + ok = true + for k in 1:n + k == j && continue + c = _arc_crossing_parity(q, m, pts[k], pts[mod1(k + 1, n)]; orient, proper_crossing) + if c == -1 + ok = false + break + end + crossings += c + end + ok || continue + # walking from `m` toward `q` departs onto `q`'s side of the anchor + # edge (the arc meets that great circle again only at `-m`); positive + # side is the interior, and each crossing flips it + return isodd(crossings) ? side < 0 : side > 0 + end + return nothing +end + +# Crossing parity of the test arc q → m against ring edge a → b: 1 for a +# transversal crossing, 0 for none, -1 for too close to degenerate to call +# (with an exact `orient`, only exact incidences return -1). +function _arc_crossing_parity(q, m, a, b; orient, proper_crossing) + # a vertex at `-q` lies on every great circle through `q`; its edges can + # reach the test arc only at `q` itself, excluded by the on-boundary check + (a == -q || b == -q) && return 0 + sa = orient(q, m, a) + sb = orient(q, m, b) + (sa == 0 || sb == 0) && return -1 + (sa > 0) == (sb > 0) && return 0 + # `q` on this edge's great circle but off the edge (checked upfront): the + # circles meet only at `±q`, out of the test arc's reach — no crossing. + # Anchor-independent (lonlat meridian edges hold `±eₓ`/`±e_y` exactly), + # so resolve instead of returning -1. + sq = orient(a, b, q) + sq == 0 && return 0 + sm = orient(a, b, m) + sm == 0 && return -1 + (sq > 0) == (sm > 0) && return 0 + return proper_crossing(q, m, a, b) +end + +# Default transversality decision: the circles' intersection direction `x` +# must point into both arcs' hemispheres (each arc holds exactly one of `±x` +# once the endpoint sides are strict). Tolerance-banded; assumes unit input. +function _hemisphere_proper_crossing(q, m, a, b) + x = cross(normalize(robust_cross_product(q, m)), + normalize(robust_cross_product(a, b))) + d1 = dot(x, q + m) + d2 = dot(x, a + b) + tol = 16 * eps(Float64) * norm(x) + (abs(d1) <= tol || abs(d2) <= tol) && return -1 + return (d1 > 0) == (d2 > 0) ? 1 : 0 +end From 06a77c2631ce7d6c020d48d4c4ecdf349d49ecbe Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 11 Jul 2026 17:04:14 -0400 Subject: [PATCH 19/20] Remove `rk_bounds_disjoint` and `rk_bounds_covers` MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `rk_bounds_disjoint` was `!Extents.intersects` under a kernel name — its call sites now say that directly. `rk_bounds_covers` had no callers (the covers short-circuit uses `ext_covers` = `Extents.covers`, which is dimension-generic over shared keys). The kernel contract now documents `rk_interaction_bounds` as the shared manifold extent plus relate glue, and the `_validate_relate_edges` ingest hook. Co-Authored-By: Claude Fable 5 --- src/methods/geom_relations/relateng/kernel.jl | 40 +++++-------------- .../geom_relations/relateng/point_locator.jl | 2 +- .../relateng/relate_geometry.jl | 4 +- .../geom_relations/relateng/relate_ng.jl | 2 +- test/methods/relateng/kernel.jl | 5 ++- .../relateng/kernel_conformance_spherical.jl | 12 ------ 6 files changed, 18 insertions(+), 47 deletions(-) diff --git a/src/methods/geom_relations/relateng/kernel.jl b/src/methods/geom_relations/relateng/kernel.jl index f5f77f4da9..b4e5d101fa 100644 --- a/src/methods/geom_relations/relateng/kernel.jl +++ b/src/methods/geom_relations/relateng/kernel.jl @@ -36,19 +36,19 @@ repeated last point): one of `LOC_INTERIOR`, `LOC_BOUNDARY`, `LOC_EXTERIOR`. rk_interaction_bounds(m, geom)::Extents.Extent -The bounding region within which `geom` can interact with another geometry. -On the plane this is the ordinary extent; other manifolds may need to widen -it (e.g. great-circle edges bulge outside the coordinate box of their -endpoints). +The bounding region within which `geom` can interact with another geometry: +the shared manifold extent (`Extents.extent(m, geom)`) in kernel +coordinates, plus any relate-specific conservatism (ulp padding, dim-1 +curve semantics for rings — see kernel_spherical.jl). Extent tests on +these boxes use `Extents.intersects`/`Extents.covers` directly (via the +`nothing`-tolerant `ext_intersects`/`ext_covers` where empty geometries +can occur). - rk_bounds_disjoint(extA, extB)::Bool - rk_bounds_covers(extA, extB)::Bool + _validate_relate_edges(m, curve) -Conservative interaction-bounds tests used for short-circuiting: -`rk_bounds_disjoint` must only return `true` when no interaction is possible; -`rk_bounds_covers` must only return `true` when `extA` covers `extB` in the -X/Y dimensions. These operate on the extents produced by -`rk_interaction_bounds` and are manifold-generic. +Ingest validation hook, run once per curve at `RelateGeometry` +construction: throw on edges the manifold cannot represent. A no-op on +`Planar`; `Spherical` rejects exactly-antipodal vertex pairs. rk_classify_intersection(m, a0, a1, b0, b1; exact)::SegSegClass @@ -153,24 +153,6 @@ end # Manifold-generic helpers -# Conservative interaction-bounds tests (contract above): manifold-generic, -# operating on the extents produced by `rk_interaction_bounds`. -# `rk_bounds_disjoint` is dimension-generic already (`Extents.intersects`). -# `rk_bounds_covers` checks X/Y on the plane and additionally Z for the 3D -# extents the spherical kernel produces; the X/Y path is byte-identical to the -# original planar definition. `hasproperty(_, :Z)` folds to a compile-time -# constant for a concretely-typed `Extent`, so both paths stay allocation-free. -rk_bounds_disjoint(extA, extB) = !Extents.intersects(extA, extB) - -function rk_bounds_covers(extA, extB) - (extA.X[1] <= extB.X[1] && extB.X[2] <= extA.X[2]) && - (extA.Y[1] <= extB.Y[1] && extB.Y[2] <= extA.Y[2]) && - _bounds_covers_z(extA, extB) -end -@inline _bounds_covers_z(extA, extB) = - (!hasproperty(extA, :Z) || !hasproperty(extB, :Z)) || - (extA.Z[1] <= extB.Z[1] && extB.Z[2] <= extA.Z[2]) - # Symbolic node identity (design D2). One concrete isbits key type for both # node kinds so Dict{NodeKey{P}, ...} is type-stable. Equality and hashing # are the default bit-pattern (egal) semantics for isbits structs; this is diff --git a/src/methods/geom_relations/relateng/point_locator.jl b/src/methods/geom_relations/relateng/point_locator.jl index 700244764c..7c1ef96edb 100644 --- a/src/methods/geom_relations/relateng/point_locator.jl +++ b/src/methods/geom_relations/relateng/point_locator.jl @@ -460,7 +460,7 @@ end function locate_on_line(loc::RelatePointLocator, p, is_node::Bool, line) #-- Java: lineEnv.intersects(p) short-circuit (p is already a kernel point) pt_ext = _kernel_point_box(p) - if rk_bounds_disjoint(rk_interaction_bounds(loc.m, line), pt_ext) + if !Extents.intersects(rk_interaction_bounds(loc.m, line), pt_ext) return LOC_EXTERIOR end #-- Java: PointLocation.isOnLine over the coordinate sequence diff --git a/src/methods/geom_relations/relateng/relate_geometry.jl b/src/methods/geom_relations/relateng/relate_geometry.jl index b6541e8b82..8058a6d486 100644 --- a/src/methods/geom_relations/relateng/relate_geometry.jl +++ b/src/methods/geom_relations/relateng/relate_geometry.jl @@ -560,7 +560,7 @@ function _extract_segment_strings_from_atomic!(rg::RelateGeometry, is_a::Bool, g if elem_ext === missing elem_ext = rk_interaction_bounds(rg.m, geom) end - rk_bounds_disjoint(ext_filter, elem_ext) && return nothing + !Extents.intersects(ext_filter, elem_ext) && return nothing end rg.element_id += Int32(1) @@ -590,7 +590,7 @@ function _extract_ring_to_segment_string!(rg::RelateGeometry, is_a::Bool, ring, if ring_ext === missing ring_ext = rk_interaction_bounds(rg.m, ring) end - rk_bounds_disjoint(ext_filter, ring_ext) && return nothing + !Extents.intersects(ext_filter, ring_ext) && return nothing end #-- orient the points if required diff --git a/src/methods/geom_relations/relateng/relate_ng.jl b/src/methods/geom_relations/relateng/relate_ng.jl index deed6d2d4c..7f3de205cb 100644 --- a/src/methods/geom_relations/relateng/relate_ng.jl +++ b/src/methods/geom_relations/relateng/relate_ng.jl @@ -769,7 +769,7 @@ end # (only reachable when the target is empty but the predicate still requires # exterior checks). _elem_env_disjoint(m::Manifold, elem, target_ext) = - target_ext === nothing || rk_bounds_disjoint(rk_interaction_bounds(m, elem), target_ext) + target_ext === nothing || !Extents.intersects(rk_interaction_bounds(m, elem), target_ext) # Java LineString.isClosed: false for an empty line, otherwise exact 2D # coordinate equality of the endpoints. diff --git a/test/methods/relateng/kernel.jl b/test/methods/relateng/kernel.jl index 9639779b75..7665c8f10f 100644 --- a/test/methods/relateng/kernel.jl +++ b/test/methods/relateng/kernel.jl @@ -6,6 +6,7 @@ import GeometryOps as GO import GeometryOps: Planar, True, False import GeometryOps.UnitSpherical: UnitSphericalPoint import GeoInterface as GI +import Extents const PT = Tuple{Float64, Float64} m = Planar() @@ -43,8 +44,8 @@ end @testset "bounds" begin pa = GI.Polygon([[(0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,0.0)]]) ea = GO.rk_interaction_bounds(m, pa) - @test !GO.rk_bounds_disjoint(ea, ea) - @test GO.rk_bounds_covers(ea, ea) + @test Extents.intersects(ea, ea) + @test Extents.covers(ea, ea) end @testset "rk_classify_intersection" begin diff --git a/test/methods/relateng/kernel_conformance_spherical.jl b/test/methods/relateng/kernel_conformance_spherical.jl index abb582ffc0..d531b05489 100644 --- a/test/methods/relateng/kernel_conformance_spherical.jl +++ b/test/methods/relateng/kernel_conformance_spherical.jl @@ -67,18 +67,6 @@ function kernel_conformance_suite_spherical(m; exact) @test e.Z[1] <= GI.z(u) <= e.Z[2] end end - @testset "rk_bounds_covers respects Z" begin - big = Extents.Extent(X = (0., 2.), Y = (0., 2.), Z = (0., 2.)) - inside = Extents.Extent(X = (0.5, 1.), Y = (0.5, 1.), Z = (0.5, 1.)) - outsideZ = Extents.Extent(X = (0.5, 1.), Y = (0.5, 1.), Z = (0.5, 3.)) - @test GO.rk_bounds_covers(big, inside) - @test !GO.rk_bounds_covers(big, outsideZ) - @test !GO.rk_bounds_disjoint(big, inside) - # the 2D covering relation is unchanged - big2 = Extents.Extent(X = (0., 2.), Y = (0., 2.)) - @test GO.rk_bounds_covers(big2, Extents.Extent(X = (0.5, 1.), Y = (0.5, 1.))) - @test !GO.rk_bounds_covers(big2, Extents.Extent(X = (0.5, 3.), Y = (0.5, 1.))) - end @testset "rk_classify_intersection: symmetry and incidence consistency" begin n_proper = 0; n_touch = 0; n_collinear = 0 for _ in 1:2000 From 8d5ea7e550e5f98118fb598eeb470dab07bc5333 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 11 Jul 2026 17:44:06 -0400 Subject: [PATCH 20/20] Tighten substrate-unification comments and helpers Share `_exactly_antipodal` between the ingest validation and `AntipodalEdgeSplit` (one spelling of the condition), route the spherical point interaction box through `GI.extent`, drop narrative from the new kernel comments, and update the extent-cache comment for `_reusable_stored_extent`. Co-Authored-By: Claude Fable 5 --- .../relateng/kernel_spherical.jl | 52 ++++++++----------- .../relateng/relate_geometry.jl | 6 +-- .../correction/antipodal_edge_split.jl | 11 ++-- src/utils/UnitSpherical/predicates.jl | 33 ++++++------ 4 files changed, 44 insertions(+), 58 deletions(-) diff --git a/src/methods/geom_relations/relateng/kernel_spherical.jl b/src/methods/geom_relations/relateng/kernel_spherical.jl index 2ae7ebf1fd..bcc092dabc 100644 --- a/src/methods/geom_relations/relateng/kernel_spherical.jl +++ b/src/methods/geom_relations/relateng/kernel_spherical.jl @@ -98,22 +98,20 @@ _kernel_point_type(::Spherical) = UnitSphericalPoint{Float64} "unique great-circle arc; densify it first with the `AntipodalEdgeSplit` " * "correction (it inserts the lon/lat midpoint)")) -# Ingest validation, run once per curve at `RelateGeometry` construction (the -# extent-cache pass). A vanishing cross product with the endpoints pointing -# opposite ways means the vertices are exactly antipodal: infinitely many -# great circles pass through them, so the edge has no well-defined arc and -# relate's contract is to throw informatively. (`spherical_arc_extent` -# deliberately does not throw — it picks a stable plane — so the guard lives -# here, not in the extent.) A vanishing cross with `a ⋅ b > 0` is a -# zero-length/repeated vertex, which is fine. +# Exactly antipodal pair: vanishing cross product, opposed directions. A +# vanishing cross with `u ⋅ v > 0` is a zero-length/repeated vertex, fine. +_exactly_antipodal(u, v) = iszero(cross(u, v)) && (u ⋅ v) < 0.0 + +# Ingest validation, once per curve at `RelateGeometry` construction: an +# exactly-antipodal edge has no unique great-circle arc, so throw rather +# than pick one (`spherical_arc_extent` picks a stable plane, never throws). function _validate_relate_edges(::Spherical, curve) n = GI.npoint(curve) n < 2 && return nothing prev = _spherical_kernel_point(GI.getpoint(curve, 1)) for i in 2:n cur = _spherical_kernel_point(GI.getpoint(curve, i)) - c = cross(prev, cur) - (c ⋅ c) == 0.0 && (prev ⋅ cur) < 0.0 && _throw_antipodal_edge(prev, cur) + _exactly_antipodal(prev, cur) && _throw_antipodal_edge(prev, cur) prev = cur end return nothing @@ -288,14 +286,11 @@ _ring_kernel_pts(::False, ring) = _ring_usp(ring) # Location of `p` relative to the area enclosed by `ring` (S2 convention: # interior on the left). Boundary first (exact arc membership), then the -# shared anchor-retry crossing parity — `spherical_ring_contains` with this -# kernel's predicates injected: `rk_orient` for sides and -# `_arcs_cross_properly` for transversality, so the interior decision is as -# exact as the predicates. Unlike a fixed reference point (the previous -# north-pole walk), the anchor retry has no preferred direction to go -# degenerate at, and handles super-hemisphere interiors. If every anchor is -# degenerate with respect to `p` — unreachable for a non-degenerate ring and -# an off-boundary point — the configuration is refused, not answered wrong. +# shared `spherical_ring_contains` with this kernel's predicates injected — +# `rk_orient` for sides, `_arcs_cross_properly` for transversality — so the +# interior decision is as exact as the predicates. All anchors degenerate +# (unreachable for a non-degenerate ring and an off-boundary point) is +# refused, not answered wrong. function rk_point_in_ring(m::Spherical, p, ring; exact) pts = _ring_kernel_pts(ring) @inbounds for i in 1:length(pts)-1 @@ -317,24 +312,23 @@ end "respect to the query point $(_tup3(p)) — the ring is degenerate at " * "this point")) -# Interaction bounds on the sphere, built from the shared substrate +# Interaction bounds on the sphere: the shared substrate # (`spherical_arc_extent` per edge, `_spherical_region_extent` for area -# interiors) over kernel-converted points, so the box and the ingested -# vertices agree bit-for-bit. Relate-specific glue: rings are a polygon's -# linework / dim-1 curve elements (JTS semantics), not S2 regions — their -# edges are bounded directly, so a CW hole cannot become a complement region — -# and every box gets a few ulps of padding so a kernel point that differs -# from another conversion path by sub-ulp noise still prunes as interacting. +# interiors) over kernel-converted points, so box and ingested vertices +# agree bit-for-bit. Rings are dim-1 linework here (JTS semantics), not S2 +# regions — a CW hole must not become a complement region. Boxes get a few +# ulps of padding so a vertex from another conversion path still prunes as +# interacting. rk_interaction_bounds(m::Spherical, geom) = _pad_bounds(_sph_interaction_extent(m, GI.trait(geom), geom)) -function _sph_interaction_extent(m::Spherical, ::GI.AbstractPointTrait, geom) - p = _spherical_kernel_point(geom) - return Extents.Extent(X = (p[1], p[1]), Y = (p[2], p[2]), Z = (p[3], p[3])) -end +_sph_interaction_extent(m::Spherical, ::GI.AbstractPointTrait, geom) = + GI.extent(_spherical_kernel_point(geom)) function _sph_interaction_extent(m::Spherical, ::GI.AbstractCurveTrait, geom) n = GI.npoint(geom) prev = _spherical_kernel_point(GI.getpoint(geom, 1)) + # seeding with pts[1]'s box covers the degenerate n == 1 curve; it is + # absorbed by the first edge box otherwise ext = spherical_arc_extent(prev, prev) for i in 2:n cur = _spherical_kernel_point(GI.getpoint(geom, i)) diff --git a/src/methods/geom_relations/relateng/relate_geometry.jl b/src/methods/geom_relations/relateng/relate_geometry.jl index 8058a6d486..d3971b79e7 100644 --- a/src/methods/geom_relations/relateng/relate_geometry.jl +++ b/src/methods/geom_relations/relateng/relate_geometry.jl @@ -124,9 +124,9 @@ bounds embedded at every level, in one coordinate pass. Wrappers share the original linework objects (a `GI.LinearRing(ring; extent)` around a same-trait geometry stores `ring`'s coordinate backing, copying nothing), so this costs O(#elements) small allocations plus the one extent scan the -constructor performed anyway. Levels that already carry a stored extent -are reused as-is, so re-wrapping an already-cached tree (or user inputs -built with `extent = ...`) does no coordinate work. +constructor performed anyway. Levels whose stored extent is usable as +interaction bounds (`_reusable_stored_extent`) are reused as-is, so +re-wrapping an already-cached tree does no coordinate work. ==========================================================================# _has_stored_extent(geom) = diff --git a/src/transformations/correction/antipodal_edge_split.jl b/src/transformations/correction/antipodal_edge_split.jl index 2465dd80fd..7f8f15a93a 100644 --- a/src/transformations/correction/antipodal_edge_split.jl +++ b/src/transformations/correction/antipodal_edge_split.jl @@ -52,14 +52,9 @@ end (::AntipodalEdgeSplit)(::GI.AbstractCurveTrait, curve) = _split_antipodal_curve(curve) # Whether the lon/lat points `p`, `q` map to exactly-antipodal unit vectors — -# the same condition the kernel's edge validation throws on (vanishing cross, -# negative dot). -function _is_antipodal_lonlat(p, q) - up = _spherical_kernel_point(p) - uq = _spherical_kernel_point(q) - n = cross(up, uq) - return iszero(n[1]) && iszero(n[2]) && iszero(n[3]) && (up ⋅ uq) < 0 -end +# the condition the kernel's edge validation throws on. +_is_antipodal_lonlat(p, q) = + _exactly_antipodal(_spherical_kernel_point(p), _spherical_kernel_point(q)) # Insert the lon/lat midpoint into every antipodal edge of a ring/line, # returning the input unchanged (no copy) when there is nothing to split. diff --git a/src/utils/UnitSpherical/predicates.jl b/src/utils/UnitSpherical/predicates.jl index 1449bde36e..3e55f14419 100644 --- a/src/utils/UnitSpherical/predicates.jl +++ b/src/utils/UnitSpherical/predicates.jl @@ -97,36 +97,33 @@ end spherical_ring_contains(pts, n, q; orient, on_arc, proper_crossing) -> Union{Bool, Nothing} Whether `q` lies in the closed region on the left of the ring `pts[1:n]` -(S2 loop convention: counterclockwise winding, interior on the left — a -clockwise ring therefore contains the complement of the region it visually -encloses). The closing edge `pts[n] → pts[1]` is implied; points on the -boundary count as contained. Returns `nothing` when every anchor edge is -degenerate with respect to `q` — callers must treat that conservatively. +(S2 loop convention: counterclockwise winding, interior on the left, so a +clockwise ring contains the complement). The closing edge `pts[n] → pts[1]` +is implied; boundary points count as contained. Returns `nothing` when +every anchor edge is degenerate with respect to `q` — callers must treat +that conservatively. Containment is decided by crossing parity, the way `S2Loop::Contains` / `InitBound` decide pole containment: which side of an anchor edge `q` falls on, flipped once per transversal crossing of the arc from the anchor's -midpoint to `q` with the other edges. Anchors degenerate with respect to -`q` are skipped and the next edge tried. +midpoint to `q` with the other edges; degenerate anchors are skipped and +the next edge tried. -The geometric predicates are injectable, so callers with stricter -requirements can substitute their own: +The geometric predicates are injectable, for callers with stricter +requirements. They receive the input points untouched (which may be +non-unit for scale-invariant predicates — the defaults assume unit input); +only the constructed reference midpoint is normalized. - `orient(a, b, c)`: sign-valued orientation of `c` against the oriented great circle through `a, b`; default [`spherical_orient`](@ref). - `on_arc(q, a, b)::Bool`: boundary membership; default [`point_on_spherical_arc`](@ref). Pass `Returns(false)` when boundary - points have already been classified. + points are already classified. - `proper_crossing(q, m, a, b)::Int`: `1` if the minor arcs `(q, m)` and `(a, b)` cross transversally in both interiors, `0` if not, `-1` for too - close to degenerate to call; only consulted once `orient` places both - endpoint pairs strictly transversally. The default decides through - `robust_cross_product` with a small tolerance band and assumes unit - input. - -The injected predicates receive the input points untouched (which may be -non-unit for scale-invariant predicates); only the constructed reference -midpoint is normalized. + close to call; consulted once `orient` places both endpoint pairs + strictly transversally. The default uses `robust_cross_product` with a + small tolerance band. """ function spherical_ring_contains(pts, n, q; orient = spherical_orient,