diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ff1e526d6..3f8c2b28c 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -18,7 +18,7 @@ jobs: fail-fast: false matrix: version: - - '1.9' + - '1.10' - '1' - 'nightly' os: @@ -33,8 +33,8 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - uses: julia-actions/cache@v1 - - name: Dev GeometryOpsCore` - run: julia --project=. -e 'using Pkg; Pkg.develop(; path = joinpath(".", "GeometryOpsCore"))' + - name: Dev GeometryOpsCore and add other packages + run: julia --project=. -e 'using Pkg; Pkg.develop(; path = joinpath(".", "GeometryOpsCore"));' - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 @@ -60,13 +60,28 @@ jobs: with: version: '1' - name: Build and add versions - run: julia --project=docs -e 'using Pkg; Pkg.develop([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore"))]); Pkg.add([PackageSpec(name = "GeoMakie", rev = "master"), PackageSpec(name="GeoInterface", rev="bugfix_vars")])' + run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])' - uses: julia-actions/julia-docdeploy@v1 with: install-package: false env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + doctests: + name: Doctests + runs-on: ubuntu-latest + permissions: + contents: write + statuses: write + actions: write + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/cache@v1 + - uses: julia-actions/setup-julia@v1 + with: + version: '1' + - name: Build and add versions + run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])' - run: | julia --project=docs -e ' using Documenter: DocMeta, doctest diff --git a/.gitignore b/.gitignore index 37d130b1a..cce699c12 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,6 @@ /docs/build/ /docs/src/source/ .vscode/ -.DS_Store \ No newline at end of file +.DS_Store + +benchmarks/Manifest.toml \ No newline at end of file diff --git a/GeometryOpsCore/Project.toml b/GeometryOpsCore/Project.toml index e3bc961d8..4922403ac 100644 --- a/GeometryOpsCore/Project.toml +++ b/GeometryOpsCore/Project.toml @@ -1,15 +1,17 @@ name = "GeometryOpsCore" uuid = "05efe853-fabf-41c8-927e-7063c8b9f013" authors = ["Anshul Singhvi ", "Rafael Schouten ", "Skylar Gering ", "and contributors"] -version = "0.1.2" +version = "0.1.3" [deps] DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" +StableTasks = "91464d47-22a1-43fe-8b7f-2d57ee82463f" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [compat] -julia = "1.9" DataAPI = "1" GeoInterface = "1.2" +StableTasks = "0.1.5" Tables = "1" +julia = "1.9" diff --git a/GeometryOpsCore/README.md b/GeometryOpsCore/README.md index adf5a85a6..67ab0d0bb 100644 --- a/GeometryOpsCore/README.md +++ b/GeometryOpsCore/README.md @@ -2,7 +2,13 @@ This is a "core" package for [GeometryOps.jl](https://github.com/JuliaGeo/GeometryOps.jl), that defines some basic primitive functions and types for GeometryOps. -Generally, you would depend on this to use either the GeometryOps types (like `Linear`, `Spherical`, etc) or the primitive functions like `apply`, `applyreduce`, `flatten`, etc. +It defines, all in all: +- Manifolds and the manifold interface +- The Algorithm type and the algorithm interface +- Low level functions like apply, applyreduce, flatten, etc. +- Common methods that should work across all geometries! + +Generally, you would depend on this to use either the GeometryOps types (like `Planar`, `Spherical`, etc) or the primitive functions like `apply`, `applyreduce`, `flatten`, etc. All of these are also accessible from GeometryOps, so it's preferable that you use GeometryOps directly. Tests are in the main GeometryOps tests, we don't have separate tests for GeometryOpsCore since it's in a monorepo structure. \ No newline at end of file diff --git a/GeometryOpsCore/src/GeometryOpsCore.jl b/GeometryOpsCore/src/GeometryOpsCore.jl index bdd8fdf76..4dd993739 100644 --- a/GeometryOpsCore/src/GeometryOpsCore.jl +++ b/GeometryOpsCore/src/GeometryOpsCore.jl @@ -6,7 +6,7 @@ import GeoInterface import GeoInterface as GI import GeoInterface: Extents -# Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`. +# Import all exported names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`. for name in names(GeoInterface) @eval using GeoInterface: $name end @@ -16,9 +16,17 @@ end using Tables using DataAPI +import StableTasks include("keyword_docs.jl") -include("types.jl") +include("constants.jl") + +include("types/manifold.jl") +include("types/algorithm.jl") +include("types/operation.jl") +include("types/exceptions.jl") +include("types/booltypes.jl") +include("types/traittarget.jl") include("apply.jl") include("applyreduce.jl") diff --git a/GeometryOpsCore/src/apply.jl b/GeometryOpsCore/src/apply.jl index 311edbc2f..a465cc83c 100644 --- a/GeometryOpsCore/src/apply.jl +++ b/GeometryOpsCore/src/apply.jl @@ -30,18 +30,11 @@ Functions like [`flip`](@ref), [`reproject`](@ref), [`transform`](@ref), even [` using the `apply` framework. Similarly, [`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the [`applyreduce`](@ref) framework. -## Docstrings - -### Functions - ```@docs; collapse=true, canonical=false apply -applyreduce ``` -=# -#= ## What is `apply`? `apply` applies some function to every geometry matching the `Target` @@ -69,7 +62,7 @@ Be careful making a union across "levels" of nesting, e.g. `Union{FeatureTrait,PolygonTrait}`, as `_apply` will just never reach `PolygonTrait` when all the polygons are wrapped in a `FeatureTrait` object. -## Embedding: +### Embedding `extent` and `crs` can be embedded in all geometries, features, and feature collections as part of `apply`. Geometries deeper than `Target` @@ -78,7 +71,7 @@ will of course not have new `extent` or `crs` embedded. - `calc_extent` signals to recalculate an `Extent` and embed it. - `crs` will be embedded as-is -## Threading +### Threading Threading is used at the outermost level possible - over an array, feature collection, or e.g. a MultiPolygonTrait where @@ -86,6 +79,22 @@ each `PolygonTrait` sub-geometry may be calculated on a different thread. Currently, threading defaults to `false` for all objects, but can be turned on by passing the keyword argument `threaded=true` to `apply`. + +Threading uses [StableTasks.jl](https://github.com/JuliaFolds2/StableTasks.jl) to provide +type-stable tasks (base Julia `Threads.@spawn` is not type stable). This is completely cost-free +and saves some allocations when running multithreaded. + +The current strategy is to launch 2 tasks for each CPU thread, to provide load balancing. We +assume Julia will manage these tasks efficiently, and we don't want to run too many tasks +since each task does have some overhead when it's created. This may need revisiting in the future, +but it's a pretty easy heuristic to use. + +## Implementation + +Literate.jl source code is below. + +*** + =# """ @@ -319,6 +328,18 @@ end using Base.Threads: nthreads, @threads, @spawn +#= +Here we used to use the compiler directive `@assume_effects :foldable` to force the compiler +to lookup through the closure. This alone makes e.g. `flip` 2.5x faster! + +But it caused inference to fail, so we've removed it. No effect on runtime so far as we can tell, +at least in Julia 1.11. +=# +@inline function _maptasks(f::F, taskrange, threaded::False)::Vector where F + map(f, taskrange) +end + + # Threading utility, modified Mason Protters threading PSA # run `f` over ntasks, where f receives an AbstractArray/range # of linear indices @@ -333,7 +354,7 @@ using Base.Threads: nthreads, @threads, @spawn # Map over the chunks tasks = map(task_chunks) do chunk # Spawn a task to process this chunk - @spawn begin + StableTasks.@spawn begin # Where we map `f` over the chunk indices map(f, chunk) end @@ -341,14 +362,4 @@ using Base.Threads: nthreads, @threads, @spawn # Finally we join the results into a new vector return mapreduce(fetch, vcat, tasks) -end -#= -Here we used to use the compiler directive `@assume_effects :foldable` to force the compiler -to lookup through the closure. This alone makes e.g. `flip` 2.5x faster! - -But it caused inference to fail, so we've removed it. No effect on runtime so far as we can tell, -at least in Julia 1.11. -=# -@inline function _maptasks(f::F, taskrange, threaded::False)::Vector where F - map(f, taskrange) end \ No newline at end of file diff --git a/GeometryOpsCore/src/applyreduce.jl b/GeometryOpsCore/src/applyreduce.jl index c6a833797..231dda070 100644 --- a/GeometryOpsCore/src/applyreduce.jl +++ b/GeometryOpsCore/src/applyreduce.jl @@ -18,6 +18,36 @@ and perform some operation on it. [`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the [`applyreduce`](@ref) framework. + +```@docs +applyreduce +``` + + +### Threading + +Threading is used at the outermost level possible - over +an array, feature collection, or e.g. a MultiPolygonTrait where +each `PolygonTrait` sub-geometry may be calculated on a different thread. + +Currently, threading defaults to `false` for all objects, but can be turned on +by passing the keyword argument `threaded=true` to `apply`. + +Threading uses [StableTasks.jl](https://github.com/JuliaFolds2/StableTasks.jl) to provide +type-stable tasks (base Julia `Threads.@spawn` is not type stable). This is completely cost-free +and saves some allocations when running multithreaded. + +The current strategy is to launch 2 tasks for each CPU thread, to provide load balancing. We +assume Julia will manage these tasks efficiently, and we don't want to run too many tasks +since each task does have some overhead when it's created. This may need revisiting in the future, +but it's a pretty easy heuristic to use. + +## Implementation + +Literate.jl source code is below. + +*** + =# """ @@ -135,7 +165,7 @@ import Base.Threads: nthreads, @threads, @spawn # Map over the chunks tasks = map(task_chunks) do chunk # Spawn a task to process this chunk - @spawn begin + StableTasks.@spawn begin # Where we map `f` over the chunk indices mapreduce(f, op, chunk; init) end @@ -144,6 +174,7 @@ import Base.Threads: nthreads, @threads, @spawn # Finally we join the results into a new vector return mapreduce(fetch, op, tasks; init) end -Base.@assume_effects :foldable function _mapreducetasks(f::F, op, taskrange, threaded::False; init) where F + +function _mapreducetasks(f::F, op, taskrange, threaded::False; init) where F mapreduce(f, op, taskrange; init) end diff --git a/GeometryOpsCore/src/constants.jl b/GeometryOpsCore/src/constants.jl new file mode 100644 index 000000000..e6c89318b --- /dev/null +++ b/GeometryOpsCore/src/constants.jl @@ -0,0 +1,8 @@ +"The semi-major axis of the WGS84 ellipsoid" +const WGS84_EARTH_SEMI_MAJOR_RADIUS = 6378137.0 + +"The inverse flattening of the WGS84 ellipsoid" +const WGS84_EARTH_INV_FLATTENING = 298.257223563 + +"The mean radius of the WGS84 ellipsoid, used for spherical manifold default" +const WGS84_EARTH_MEAN_RADIUS = 6371008.8 \ No newline at end of file diff --git a/GeometryOpsCore/src/types.jl b/GeometryOpsCore/src/types.jl deleted file mode 100644 index c923b12dc..000000000 --- a/GeometryOpsCore/src/types.jl +++ /dev/null @@ -1,166 +0,0 @@ -#= -# Types - -This defines core types that the GeometryOps ecosystem uses, -and that are usable in more than just GeometryOps. - -=# - -#= -## `Manifold` - -A manifold is mathematically defined as a topological space that resembles Euclidean space locally. - -In GeometryOps (and geodesy more generally), there are three manifolds we care about: -- [`Planar`](@ref): the 2d plane, a completely Euclidean manifold -- [`Spherical`](@ref): the unit sphere, but one where areas are multiplied by the radius of the Earth. This is not Euclidean globally, but all map projections attempt to represent the sphere on the Euclidean 2D plane to varying degrees of success. -- [`Geodesic`](@ref): the ellipsoid, the closest we can come to representing the Earth by a simple geometric shape. Parametrized by `semimajor_axis` and `inv_flattening`. - -Generally, we aim to have `Linear` and `Spherical` be operable everywhere, whereas `Geodesic` will only apply in specific circumstances. -Currently, those circumstances are `area` and `segmentize`, but this could be extended with time and https://github.com/JuliaGeo/SphericalGeodesics.jl. -=# - -export Planar, Spherical, Geodesic -export TraitTarget -export BoolsAsTypes, True, False, booltype - -""" - abstract type Manifold - -A manifold is mathematically defined as a topological space that resembles Euclidean space locally. - -We use the manifold definition to define the space in which an operation should be performed, or where a geometry lies. - -Currently we have [`Planar`](@ref), [`Spherical`](@ref), and [`Geodesic`](@ref) manifolds. -""" -abstract type Manifold end - -""" - Planar() - -A planar manifold refers to the 2D Euclidean plane. - -Z coordinates may be accepted but will not influence geometry calculations, which -are done purely on 2D geometry. This is the standard "2.5D" model used by e.g. GEOS. -""" -struct Planar <: Manifold -end - -""" - Spherical(; radius) - -A spherical manifold means that the geometry is on the 3-sphere (but is represented by 2-D longitude and latitude). - -By default, the radius is defined as the mean radius of the Earth, `6371008.8 m`. - -## Extended help - -!!! note - The traditional definition of spherical coordinates in physics and mathematics, - ``r, \\theta, \\phi``, uses the _colatitude_, that measures angular displacement from the `z`-axis. - - Here, we use the geographic definition of longitude and latitude, meaning - that `lon` is longitude between -180 and 180, and `lat` is latitude between - `-90` (south pole) and `90` (north pole). -""" -Base.@kwdef struct Spherical{T} <: Manifold - radius::T = 6371008.8 -end - -""" - Geodesic(; semimajor_axis, inv_flattening) - -A geodesic manifold means that the geometry is on a 3-dimensional ellipsoid, parameterized by `semimajor_axis` (``a`` in mathematical parlance) -and `inv_flattening` (``1/f``). - -Usually, this is only relevant for area and segmentization calculations. It becomes more relevant as one grows closer to the poles (or equator). -""" -Base.@kwdef struct Geodesic{T} <: Manifold - semimajor_axis::T = 6378137.0 - inv_flattening::T = 298.257223563 -end - -#= -## `TraitTarget` - -This struct holds a trait parameter or a union of trait parameters. -It's essentially a way to construct unions. -=# - -""" - TraitTarget{T} - -This struct holds a trait parameter or a union of trait parameters. - -It is primarily used for dispatch into methods which select trait levels, -like `apply`, or as a parameter to `target`. - -## Constructors -```julia -TraitTarget(GI.PointTrait()) -TraitTarget(GI.LineStringTrait(), GI.LinearRingTrait()) # and other traits as you may like -TraitTarget(TraitTarget(...)) -# There are also type based constructors available, but that's not advised. -TraitTarget(GI.PointTrait) -TraitTarget(Union{GI.LineStringTrait, GI.LinearRingTrait}) -# etc. -``` - -""" -struct TraitTarget{T} end -TraitTarget(::Type{T}) where T = TraitTarget{T}() -TraitTarget(::T) where T<:GI.AbstractTrait = TraitTarget{T}() -TraitTarget(::TraitTarget{T}) where T = TraitTarget{T}() -TraitTarget(::Type{<:TraitTarget{T}}) where T = TraitTarget{T}() -TraitTarget(traits::GI.AbstractTrait...) = TraitTarget{Union{map(typeof, traits)...}}() - - -Base.in(::Trait, ::TraitTarget{Target}) where {Trait <: GI.AbstractTrait, Target} = Trait <: Target - - - -#= -## `BoolsAsTypes` - -In `apply` and `applyreduce`, we pass `threading` and `calc_extent` as types, not simple boolean values. - -This is to help compilation - with a type to hold on to, it's easier for -the compiler to separate threaded and non-threaded code paths. - -Note that if we didn't include the parent abstract type, this would have been really -type unstable, since the compiler couldn't tell what would be returned! - -We had to add the type annotation on the `booltype(::Bool)` method for this reason as well. - -TODO: should we switch to `Static.jl`? -=# - -""" - abstract type BoolsAsTypes - -""" -abstract type BoolsAsTypes end - -""" - struct True <: BoolsAsTypes - -A struct that means `true`. -""" -struct True <: BoolsAsTypes end - -""" - struct False <: BoolsAsTypes - -A struct that means `false`. -""" -struct False <: BoolsAsTypes end - -""" - booltype(x) - -Returns a [`BoolsAsTypes`](@ref) from `x`, whether it's a boolean or a BoolsAsTypes. -""" -function booltype end - -@inline booltype(x::Bool)::BoolsAsTypes = x ? True() : False() -@inline booltype(x::BoolsAsTypes)::BoolsAsTypes = x diff --git a/GeometryOpsCore/src/types/algorithm.jl b/GeometryOpsCore/src/types/algorithm.jl new file mode 100644 index 000000000..5fb21c4a3 --- /dev/null +++ b/GeometryOpsCore/src/types/algorithm.jl @@ -0,0 +1,78 @@ +#= +# `Algorithm`s + +An `Algorithm` is a type that describes the algorithm used to perform some [`Operation`](@ref). + +An algorithm may be associated with one or many [`Manifold`](@ref)s. It may either have the manifold as a field, or have it as a static parameter (e.g. `struct GEOS <: Manifold{Planar}`). + + +Algorithms are: +* Ways to perform an operation +* For example: LHuilier, Bessel, Ericsson for spherical area +* May be manifold agnostic (like simplification) or restrictive (like GEOS only works on planar, PROJ algorithm for arclength and area only works on geodesic) +* May or may not carry manifolds around, but manifold should always be accessible from manifold(alg) - it's not necessary that fixed manifold args can skip carrying the manifold around, eg in the case of Proj{Geodesic}. + + +## Single manifold vs manifold independent algorithms + +Some algorithms only work on a single manifold (shoelace area only works on `Planar`) +and others work on any manifold (`simplify`, `segmentize`). They are allowed to dispatch +on the manifold type but store that manifold as a field, and the fundamental algorithm is the +same. For example the segmentize algorithm is the same (distance along line) but the implementation +varies slightly depending on the manifold (planar, spherical, geodesic, etc). + +Here's a simple example of two algorithms, one only on planar and one manifold independent: + +```julia + +struct MyExternalArbitraryPackageAlgorithm <: SingleManifoldAlgorithm{Planar} + kw1::Int + kw2::String +end # this already has the methods specified + +struct MyIndependentAlgorithm{M <: Manifold} <: ManifoldIndependentAlgorithm{M} + m::M + kw1::Int + kw2::String +end + +MyIndependentAlgorithm(m::Manifold; kw1 = 1, kw2 = "hello") = MyIndependentAlgorithm(m, kw1, kw2) +``` +=# + +export Algorithm, AutoAlgorithm, ManifoldIndependentAlgorithm, SingleManifoldAlgorithm, NoAlgorithm + +abstract type Algorithm{M <: Manifold} end + +struct AutoAlgorithm{T, M <: Manifold} <: Algorithm{M} + manifold::M + x::T +end + +AutoAlgorithm(m::Manifold; kwargs...) = AutoAlgorithm(m, kwargs) +AutoAlgorithm(; kwargs...) = AutoAlgorithm(AutoManifold(), kwargs) + + +abstract type ManifoldIndependentAlgorithm{M <: Manifold} <: Algorithm{M} end + +abstract type SingleManifoldAlgorithm{M <: Manifold} <: Algorithm{M} end + +struct NoAlgorithm{M <: Manifold} <: Algorithm{M} + m::M +end + +NoAlgorithm() = NoAlgorithm(Planar()) # TODO: add a NoManifold or AutoManifold type? +# Maybe AutoManifold +# and then we have DD.format like materialization + +function (Alg::Type{<: SingleManifoldAlgorithm{M}})(m::M; kwargs...) where {M} + # successful - the algorithm is designed for this manifold + # in this case, just return `Alg(; kwargs...)` + return Alg(; kwargs...) +end + +function (Alg::Type{<: SingleManifoldAlgorithm{M}})(m::Manifold; kwargs...) where {M} + # this catches the case where the algorithm doesn't match the manifold + # throw a WrongManifoldException and be done with it + throw(WrongManifoldException{typeof(m), M, Alg}()) +end diff --git a/GeometryOpsCore/src/types/booltypes.jl b/GeometryOpsCore/src/types/booltypes.jl new file mode 100644 index 000000000..6089f0bbb --- /dev/null +++ b/GeometryOpsCore/src/types/booltypes.jl @@ -0,0 +1,53 @@ +#= +# `BoolsAsTypes` + +In `apply` and `applyreduce`, we pass `threading` and `calc_extent` as types, not simple boolean values. + +This is to help compilation - with a type to hold on to, it's easier for +the compiler to separate threaded and non-threaded code paths. + +Note that if we didn't include the parent abstract type, this would have been really +type unstable, since the compiler couldn't tell what would be returned! + +We had to add the type annotation on the `booltype(::Bool)` method for this reason as well. + + +!!! note Static.jl + + Static.jl is a package that provides a way to store and manipulate static values. + But it creates a lot of invalidations since it breaks the assumption that operations + like `<`, `>` and `==` can only return booleans. So we don't use it here. + +=# + +export BoolsAsTypes, True, False, booltype + +""" + abstract type BoolsAsTypes + +""" +abstract type BoolsAsTypes end + +""" + struct True <: BoolsAsTypes + +A struct that means `true`. +""" +struct True <: BoolsAsTypes end + +""" + struct False <: BoolsAsTypes + +A struct that means `false`. +""" +struct False <: BoolsAsTypes end + +""" + booltype(x) + +Returns a [`BoolsAsTypes`](@ref) from `x`, whether it's a boolean or a BoolsAsTypes. +""" +function booltype end + +@inline booltype(x::Bool)::BoolsAsTypes = x ? True() : False() +@inline booltype(x::BoolsAsTypes)::BoolsAsTypes = x diff --git a/GeometryOpsCore/src/types/exceptions.jl b/GeometryOpsCore/src/types/exceptions.jl new file mode 100644 index 000000000..d9a3f8d82 --- /dev/null +++ b/GeometryOpsCore/src/types/exceptions.jl @@ -0,0 +1,45 @@ +#= +# Errors and exceptions + +We create a few custom exception types in this file, +that have nice show methods that we can use for certain errors. + +This makes it substantially easier to catch specific kinds of errors and show them. +For example, we can catch `WrongManifoldException` and show a nice error message, +and error hinters can be specialized to that as well. +=# + +export WrongManifoldException + +""" + WrongManifoldException{InputManifold, DesiredManifold, Algorithm} <: Exception + +This error is thrown when an `Algorithm` is called with a manifold that it was not designed for. + +It's mainly called for [`SingleManifoldAlgorithm`](@ref) types. +""" +struct WrongManifoldException{InputManifold, DesiredManifold, Algorithm} <: Base.Exception + description::String +end + +WrongManifoldException{I, D, A}() where {I, D, A} = WrongManifoldException{I, D, A}("") + +function Base.showerror(io::IO, e::WrongManifoldException{I,D,A}) where {I,D,A} + print(io, "Algorithm ") + printstyled(io, A; bold = true, color = :green) + print(io, " is only compatible with manifold ") + printstyled(io, D; bold = true, color = :blue) + print(io, ",\n but it was called with manifold ") + printstyled(io, I; bold = true, color = :red) + print(io, ".") + + println(io, """ + \n + To fix this issue, you can specify the manifold explicitly, + e.g. `$A($D(); kwargs...)`, when constructing the algorithm. + """) + if !isempty(e.description) + print(io, "\n\n") + print(io, e.description) + end +end \ No newline at end of file diff --git a/GeometryOpsCore/src/types/manifold.jl b/GeometryOpsCore/src/types/manifold.jl new file mode 100644 index 000000000..345036e9b --- /dev/null +++ b/GeometryOpsCore/src/types/manifold.jl @@ -0,0 +1,104 @@ + +#= +# `Manifold`s + +A manifold is mathematically defined as a topological space that resembles Euclidean space locally. + +In GeometryOps, this represents the domain or space on which your geometries live. +Manifolds can be accessible via the crs info of the geometry - OR can be specified explicitly. +For example you may pass planar geometry around using GeoJSON, but because the spec says GeoJSON is only geographic, +GeometryOps will interpret GeoJSON geometries as geographic on WGS84, unless told otherwise. + +In GeometryOps (and geodesy more generally), there are three manifolds we care about: +- [`Planar`](@ref): the 2d plane, a completely Euclidean manifold +- [`Spherical`](@ref): the unit sphere, but one where areas are multiplied by the radius of the Earth. This is not Euclidean globally, but all map projections attempt to represent the sphere on the Euclidean 2D plane to varying degrees of success. +- [`Geodesic`](@ref): the ellipsoid, the closest we can come to representing the Earth by a simple geometric shape. Parametrized by `semimajor_axis` and `inv_flattening`. +- [`AutoManifold`](@ref): a special manifold that automatically selects the best manifold for the operation when it's executed. Resolves to [`Planar`](@ref), [`Spherical`](@ref), or [`Geodesic`](@ref) depending on the input geometry. + +Generally, we aim to have `Linear` and `Spherical` be operable everywhere, whereas `Geodesic` will only apply in specific circumstances. +Currently, those circumstances are [`area`](@ref), [`arclength`](@ref), and [`segmentize`](@ref), but this could be extended with time and https://github.com/JuliaGeo/SphericalGeodesics.jl. +=# + +export Manifold, AutoManifold, Planar, Spherical, Geodesic + +""" + abstract type Manifold + +A manifold is mathematically defined as a topological space that resembles Euclidean space locally. + +We use the manifold definition to define the space in which an operation should be performed, or where a geometry lies. + +Currently we have [`Planar`](@ref), [`Spherical`](@ref), and [`Geodesic`](@ref) manifolds. +""" +abstract type Manifold end + +""" + AutoManifold() + +The `AutoManifold` is a special manifold that automatically selects the best manifold for the operation. +It does not carry any parameters, nor does it indicate anything about the nature of the space. + +This gets resolved to a specific manifold when an operation is applied, using the `format` method. +""" +struct AutoManifold <: Manifold end + +""" + Planar() + +A planar manifold refers to the 2D Euclidean plane. + +Z coordinates may be accepted but will not influence geometry calculations, which +are done purely on 2D geometry. This is the standard "2.5D" model used by e.g. GEOS. +""" +struct Planar <: Manifold +end + +""" + Spherical(; radius) + +A spherical manifold means that the geometry is on the 3-sphere (but is represented by 2-D longitude and latitude). + +## Extended help + +!!! note + The traditional definition of spherical coordinates in physics and mathematics, + ``r, \\theta, \\phi``, uses the _colatitude_, that measures angular displacement from the `z`-axis. + + Here, we use the geographic definition of longitude and latitude, meaning + that `lon` is longitude between -180 and 180, and `lat` is latitude between + `-90` (south pole) and `90` (north pole). +""" +Base.@kwdef struct Spherical{T} <: Manifold + radius::T = WGS84_EARTH_MEAN_RADIUS # this should be theWGS84 defined mean radius +end + +""" + Geodesic(; semimajor_axis, inv_flattening) + +A geodesic manifold means that the geometry is on a 3-dimensional ellipsoid, parameterized by `semimajor_axis` (``a`` in mathematical parlance) +and `inv_flattening` (``1/f``). + +Usually, this is only relevant for area and segmentization calculations. It becomes more relevant as one grows closer to the poles (or equator). +""" +Base.@kwdef struct Geodesic{T} <: Manifold + semimajor_axis::T = WGS84_EARTH_SEMI_MAJOR_RADIUS # WGS84 by default + inv_flattening::T = WGS84_EARTH_INV_FLATTENING # WGS84 by default +end + + +# specifically for manifolds +# not used now but will be used later +abstract type EllipsoidParametrization end + +struct SemimajorAxisInvFlattening{T} <: EllipsoidParametrization + semimajor_axis::T + inv_flattening::T +end + +# this should be the full Ellipsoid parametrization from +struct FullEllipsoidParametrization{T} <: EllipsoidParametrization + semimajor_axis::T + semiminor_axis::T + inv_flattening::T +end + diff --git a/GeometryOpsCore/src/types/operation.jl b/GeometryOpsCore/src/types/operation.jl new file mode 100644 index 000000000..a5841c0ca --- /dev/null +++ b/GeometryOpsCore/src/types/operation.jl @@ -0,0 +1,64 @@ +#= + +# Operations + +!!! warning + Operations are not yet implemented. If you want to implement them then you may do so at your own risk - or file a PR! + +Operations are callable structs, that contain the entire specification for what the algorithm will do. + +Sometimes they may be underspecified and only materialized fully when you see the geometry, so you can extract +the best manifold for those geometries. + +* Some conceptual thing that you do to a geometry +* Overloads on abstract type to decompose user input to have materialized algorithm, manifold, and geoms +* Run Operation{Alg{Manifold}}(trait, geom) at the lowest level +* Some indication on whether to use apply or applyreduce? Or are we going too far here + * if we do this, then we also need `operation_level` to return a geometry trait or traittarget + +Operations may look like: + +```julia +Arclength()(geoms) +Arclength(Geodesic())(geoms) +Arclength(Proj())(geoms) +Arclength(Proj(Geodesic(; ...)))(geoms) +Arclength(Ericsson())(geoms) # more precise, goes wonky if any points in a triangle are antipodal +Arclength(LHuilier())(geoms) # less precise, does not go wonky on antipodal points +``` + +Two argument operations, like polygon set operations, may look like: + +```julia +Union(intersection_alg(manifold); exact, target)(geom1, geom2) +``` + +Here intersection_alg can be Foster, which we already have in GeometryOps, or GEOS +but if we ever implement e.g. RelateNG in GeometryOps, we can add that in. +=# + +abstract type Operation{Alg <: Algorithm} end + +#= +Here's an example of how this might work: + +```julia +struct XPlusOneOperation{M <: Manifold} <: Operation{NoAlgorithm{M}} + m::M # the manifold always needs to be stored, since it's not a singleton + x::Int +end +``` + +```julia +struct Area{Alg<: Algorithm} <: Operation{Alg} + alg::Alg +end + +Area() = Area(NoAlgorithm(Planar())) + +function (op::Area{Alg})(data; threaded = _False(), init = 0.0) where {Alg <: Algorithm} + return GO.area(op.alg, data; threaded, init) +end +``` + +=# \ No newline at end of file diff --git a/GeometryOpsCore/src/types/traittarget.jl b/GeometryOpsCore/src/types/traittarget.jl new file mode 100644 index 000000000..9205dbd89 --- /dev/null +++ b/GeometryOpsCore/src/types/traittarget.jl @@ -0,0 +1,39 @@ + +#= +## `TraitTarget` + +This struct holds a trait parameter or a union of trait parameters. +It's essentially a way to construct unions. +=# + +export TraitTarget + +""" + TraitTarget{T} + +This struct holds a trait parameter or a union of trait parameters. + +It is primarily used for dispatch into methods which select trait levels, +like `apply`, or as a parameter to `target`. + +## Constructors +```julia +TraitTarget(GI.PointTrait()) +TraitTarget(GI.LineStringTrait(), GI.LinearRingTrait()) # and other traits as you may like +TraitTarget(TraitTarget(...)) +# There are also type based constructors available, but that's not advised. +TraitTarget(GI.PointTrait) +TraitTarget(Union{GI.LineStringTrait, GI.LinearRingTrait}) +# etc. +``` + +""" +struct TraitTarget{T} end +TraitTarget(::Type{T}) where T = TraitTarget{T}() +TraitTarget(::T) where T<:GI.AbstractTrait = TraitTarget{T}() +TraitTarget(::TraitTarget{T}) where T = TraitTarget{T}() +TraitTarget(::Type{<:TraitTarget{T}}) where T = TraitTarget{T}() +TraitTarget(traits::GI.AbstractTrait...) = TraitTarget{Union{map(typeof, traits)...}}() + + +Base.in(::Trait, ::TraitTarget{Target}) where {Trait <: GI.AbstractTrait, Target} = Trait <: Target diff --git a/Project.toml b/Project.toml index 333ac3c08..285ff9797 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.15" +version = "0.1.16" [deps] CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" @@ -34,7 +34,7 @@ ExactPredicates = "2.2.8" FlexiJoins = "0.1.30" GeoInterface = "1.2" GeometryBasics = "0.4.7, 0.5" -GeometryOpsCore = "=0.1.2" +GeometryOpsCore = "=0.1.3" LibGEOS = "0.9.2" LinearAlgebra = "1" Proj = "1" diff --git a/docs/Project.toml b/docs/Project.toml index fccffcc39..e6eb64742 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -27,6 +27,7 @@ GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6" GeoParquet = "e99870d8-ce00-4fdd-aeee-e09192881159" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab" +GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013" LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" @@ -41,3 +42,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9" + +[sources] +GeometryOps = { path = "../" } +GeometryOpsCore = { path = "../GeometryOpsCore" } \ No newline at end of file diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index c2d6c4aad..152899db9 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -5,7 +5,8 @@ module GeometryOps import GeometryOpsCore import GeometryOpsCore: TraitTarget, - Manifold, Planar, Spherical, Geodesic, + Manifold, Planar, Spherical, Geodesic, AutoManifold, WrongManifoldException, + Algorithm, AutoAlgorithm, ManifoldIndependentAlgorithm, SingleManifoldAlgorithm, NoAlgorithm, BoolsAsTypes, True, False, booltype, apply, applyreduce, flatten, reconstruct, rebuild, unwrap, _linearring, diff --git a/test/core/algorithm.jl b/test/core/algorithm.jl new file mode 100644 index 000000000..505a03909 --- /dev/null +++ b/test/core/algorithm.jl @@ -0,0 +1,29 @@ +using GeometryOpsCore + +@testset "Constructing NoAlgorithm" begin + @test NoAlgorithm() isa NoAlgorithm +end + +@testset "Constructing AutoAlgorithm" begin + @test AutoAlgorithm() isa AutoAlgorithm + @test AutoAlgorithm(; x = 1) == AutoAlgorithm(AutoManifold(), pairs((; x = 1))) +end + +@testset "SingleManifoldAlgorithm" begin + struct TestSMAlgorithm <: SingleManifoldAlgorithm{Planar} + end + + @test TestSMAlgorithm() isa TestSMAlgorithm + + @test_throws GeometryOpsCore.WrongManifoldException TestSMAlgorithm(Geodesic()) +end + +@testset "ManifoldIndependentAlgorithm" begin + struct TestMIDAlgorithm{M} <: ManifoldIndependentAlgorithm{M} + m::M + end + + @test TestMIDAlgorithm(Planar()) isa TestMIDAlgorithm + @test TestMIDAlgorithm(Spherical()) isa TestMIDAlgorithm + @test TestMIDAlgorithm(Geodesic()) isa TestMIDAlgorithm +end diff --git a/test/core/manifold.jl b/test/core/manifold.jl new file mode 100644 index 000000000..05ea182bf --- /dev/null +++ b/test/core/manifold.jl @@ -0,0 +1 @@ +# TODO: test manifold conversion from Geodesic to Spherical and vice versa \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 3fc144484..07b08fea2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,11 @@ using SafeTestsets include("helpers.jl") +@testset "Core" begin + @safetestset "Algorithm" begin include("core/algorithm.jl") end + @safetestset "Manifold" begin include("core/manifold.jl") end +end + @safetestset "Primitives" begin include("primitives.jl") end # Methods @safetestset "Angles" begin include("methods/angles.jl") end