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: diff --git a/Project.toml b/Project.toml index 35c043dfa2..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.40" +version = "0.1.42" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -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/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)) 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/GeometryOps.jl b/src/GeometryOps.jl index b13c1eb681..e6a2f15288 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 @@ -67,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/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/src/methods/extent.jl b/src/methods/extent.jl new file mode 100644 index 0000000000..5e07c677df --- /dev/null +++ b/src/methods/extent.jl @@ -0,0 +1,92 @@ +# # 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. + +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`, 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; 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 + +```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 + + 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 = 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])) + end + end + return Extents.Extent(X = (lo[1], hi[1]), Y = (lo[2], hi[2]), Z = (lo[3], hi[3])) +end 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..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 @@ -261,6 +243,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..bcc092dabc 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,25 @@ _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 +# 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)) + _exactly_antipodal(prev, cur) && _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 @@ -283,9 +260,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 @@ -300,34 +275,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 @@ -337,71 +284,86 @@ _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 `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 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 -# 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])) +@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: the shared substrate +# (`spherical_arc_extent` per edge, `_spherical_region_extent` for area +# 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)) + +_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)) + 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/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 ce5577f1cb..d3971b79e7 100644 --- a/src/methods/geom_relations/relateng/relate_geometry.jl +++ b/src/methods/geom_relations/relateng/relate_geometry.jl @@ -124,34 +124,45 @@ 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) = 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) @@ -548,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) @@ -578,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/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/src/transformations/correction/antipodal_edge_split.jl b/src/transformations/correction/antipodal_edge_split.jl index 1ca9f11737..7f8f15a93a 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,13 +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 `arc_extent` 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/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/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/src/utils/UnitSpherical/UnitSpherical.jl b/src/utils/UnitSpherical/UnitSpherical.jl index c6293b64c9..fd92b3e498 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 @@ -20,22 +21,25 @@ 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_ring_contains, spherical_arc_intersection, ArcIntersectionResult, arc_cross, arc_hinge, arc_overlap, arc_disjoint, + spherical_arc_extent, 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/src/utils/UnitSpherical/arc_extent.jl b/src/utils/UnitSpherical/arc_extent.jl new file mode 100644 index 0000000000..63302c5afe --- /dev/null +++ b/src/utils/UnitSpherical/arc_extent.jl @@ -0,0 +1,63 @@ +# # Spherical arc extents + +#= +```@docs; canonical=false +spherical_arc_extent +``` + +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. +=# + +""" + 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, 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 + +```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/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/src/utils/UnitSpherical/predicates.jl b/src/utils/UnitSpherical/predicates.jl index e0e9d36d74..3e55f14419 100644 --- a/src/utils/UnitSpherical/predicates.jl +++ b/src/utils/UnitSpherical/predicates.jl @@ -92,3 +92,110 @@ 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, 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; degenerate anchors are skipped and +the next edge tried. + +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 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 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, + 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 diff --git a/src/utils/utils.jl b/src/utils/utils.jl index ddefac1a48..d270d1d9ee 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 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 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,18 @@ 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 carries the 3D +[`UnitSpherical.spherical_arc_extent`](@ref) of its great-circle arc. """ -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 +238,23 @@ function _lineedge(ps::Tuple, ::Type{T}) where T e = GI.extent(l) return GI.Line(l.geom; extent=e) end +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/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 diff --git a/test/methods/extent.jl b/test/methods/extent.jl new file mode 100644 index 0000000000..2ab76cfbaa --- /dev/null +++ b/test/methods/extent.jl @@ -0,0 +1,158 @@ +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 = [(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), + 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 "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 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 + 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) + @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 + # 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 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/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 7547405959..d531b05489 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.)]) @@ -79,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 @@ -274,7 +250,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.)) diff --git a/test/runtests.jl b/test/runtests.jl index da3b50e888..bee2c88e94 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 @@ -34,7 +35,9 @@ end @safetestset "RelateNG" begin include("methods/relateng/runtests.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 @safetestset "Coverage" begin include("methods/clipping/coverage.jl") end @safetestset "Cut" begin include("methods/clipping/cut.jl") end 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]) 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 diff --git a/test/utils/unitspherical.jl b/test/utils/unitspherical.jl index 1bba6b8815..890b8ebf4b 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,90 @@ 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 + +@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 + + # 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 + @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 + @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 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)) + end + end +end diff --git a/test/utils/utils.jl b/test/utils/utils.jl index 79ddf4e301..cf3d9f02d2 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,38 @@ 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() 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) + @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)} + + # 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)) + 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