GeometryOps.jl

Appearance

Full GeometryOps API documentation ​

Warning

This page is still very much WIP!

Documentation for GeometryOps's full API (only for reference!).

apply and associated functions ​

GeometryOpsCore.apply Function
julia
apply(f, target::Union{TraitTarget, GI.AbstractTrait}, obj; kw...)

Reconstruct a geometry, feature, feature collection, or nested vectors of either using the function f on the target trait.

f(target_geom) => x where x also has the target trait, or a trait that can be substituted. For example, swapping PolgonTrait to MultiPointTrait will fail if the outer object has MultiPolygonTrait, but should work if it has FeatureTrait.

Objects "shallower" than the target trait are always completely rebuilt, like a Vector of FeatureCollectionTrait of FeatureTrait when the target has PolygonTrait and is held in the features. These will always be GeoInterface geometries/feature/feature collections. But "deeper" objects may remain unchanged or be whatever GeoInterface compatible objects f returns.

The result is a functionally similar geometry with values depending on f.

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

  • calc_extent: true or false. Whether to calculate the extent. Defaults to false.

Example

Flipped point the order in any feature or geometry, or iterables of either:

julia
import GeoInterface as GI
import GeometryOps as GO
geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]),
                   GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])

flipped_geom = GO.apply(GI.PointTrait, geom) do p
    (GI.y(p), GI.x(p))
end
source
GeometryOpsCore.applyreduce Function
julia
applyreduce(f, op, target::Union{TraitTarget, GI.AbstractTrait}, obj; threaded, init, kw...)

Apply function f to all objects with the target trait, and reduce the result with an op like +.

The order and grouping of application of op is not guaranteed.

If threaded==true threads will be used over arrays and iterables, feature collections and nested geometries.

init specifies the initial value for the reduction. If not provided, the reduction uses the first result as the starting point (like reduce without init). For operations like vcat, you typically don't need to provide init. For numeric reductions like +, you may want to provide init=zero(T) to ensure type stability.

source
GeometryOps.reproject Function
julia
reproject(geometry; source_crs, target_crs, transform, always_xy, time)
reproject(geometry, source_crs, target_crs; always_xy, time)
reproject(geometry, transform; always_xy, time)

Reproject any GeoInterface.jl compatible geometry from source_crs to target_crs.

The returned object will be constructed from GeoInterface wrapper geometries, wrapping Vector{NTuple{D, Float64}}, where D is the dimension.

Tip

The Proj.jl package must be loaded for this method to work, since it is implemented in a package extension.

Arguments

  • geometry: Any GeoInterface.jl compatible geometries.

  • source_crs: the source coordinate reference system, as a GeoFormatTypes.jl object or a string.

  • target_crs: the target coordinate reference system, as a GeoFormatTypes.jl object or a string.

If these a passed as keywords, transform will take priority. Without it target_crs is always needed, and source_crs is needed if it is not retrievable from the geometry with GeoInterface.crs(geometry).

Keywords

  • always_xy: force x, y coordinate order, true by default. false will expect and return points in the crs coordinate order.

  • time: the time for the coordinates. Inf by default.

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

  • calc_extent: true or false. Whether to calculate the extent. Defaults to false.

source
GeometryOps.transform Function
julia
transform(f, obj)

Apply a function f to all the points in obj.

Points will be passed to f as an SVector to allow using CoordinateTransformations.jl and Rotations.jl without hassle.

SVector is also a valid GeoInterface.jl point, so will work in all GeoInterface.jl methods.

Example

julia
julia> import GeoInterface as GI

julia> import GeometryOps as GO

julia> geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]);

julia> f = CoordinateTransformations.Translation(3.5, 1.5)
Translation(3.5, 1.5)

julia> GO.transform(f, geom)
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Linea
rRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}(StaticArraysCo
re.SVector{2, Float64}[[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}(StaticA
rraysCore.SVector{2, Float64}[[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]], nothing, nothing)], nothing, nothing)

With Rotations.jl you need to actually multiply the Rotation by the SVector point, which is easy using an anonymous function.

julia
julia> using Rotations

julia> GO.transform(p -> one(RotMatrix{2}) * p, geom)
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearR
ing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}(StaticArraysCore.SVe
ctor{2, Int64}[[2, 1], [4, 3], [6, 5], [2, 1]], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}(StaticArraysCore.SVector{2, Int64
}[[4, 3], [6, 5], [7, 6], [4, 3]], nothing, nothing)], nothing, nothing)
source

General geometry methods ​

OGC methods ​

GeometryOps.contains Function
julia
contains(g1::AbstractGeometry, g2::AbstractGeometry)::Bool

Return true if the second geometry is completely contained by the first geometry. The interiors of both geometries must intersect and the interior and boundary of the secondary (g2) must not intersect the exterior of the first (g1).

contains returns the exact opposite result of within.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
point = GI.Point((1, 2))

GO.contains(line, point)
# output
true
source
julia
contains(g1)

Return a function that checks if its input contains g1. This is equivalent to x -> contains(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.coveredby Function
julia
coveredby([manifold::Manifold], g1, g2)::Bool

Return true if the first geometry is completely covered by the second geometry. The interior and boundary of the primary geometry (g1) must not intersect the exterior of the secondary geometry (g2).

Furthermore, coveredby returns the exact opposite result of covers. They are equivalent with the order of the arguments swapped.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
p1 = GI.Point(0.0, 0.0)
p2 = GI.Point(1.0, 1.0)
l1 = GI.Line([p1, p2])

GO.coveredby(p1, l1)
# output
true
source
julia
coveredby(g1)

Return a function that checks if its input is covered by g1. This is equivalent to x -> coveredby(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.covers Function
julia
covers(g1::AbstractGeometry, g2::AbstractGeometry)::Bool

Return true if the first geometry is completely covers the second geometry, The exterior and boundary of the second geometry must not be outside of the interior and boundary of the first geometry. However, the interiors need not intersect.

covers returns the exact opposite result of coveredby.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
l1 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)])
l2 = GI.LineString([(1.0, 1.0), (1.0, 2.0)])

GO.covers(l1, l2)
# output
true
source
julia
covers(g1)

Return a function that checks if its input covers g1. This is equivalent to x -> covers(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.crosses Function
julia
 crosses(geom1, geom2)::Bool

Return true if the intersection results in a geometry whose dimension is one less than the maximum dimension of the two source geometries and the intersection set is interior to both source geometries.

TODO: broken

Examples

julia
import GeoInterface as GI, GeometryOps as GO
# TODO: Add working example
source
julia
crosses(g1)

Return a function that checks if its input crosses g1. This is equivalent to x -> crosses(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.disjoint Function
julia
disjoint([manifold::Manifold], geom1, geom2)::Bool

Return true if the first geometry is disjoint from the second geometry.

Return true if the first geometry is disjoint from the second geometry. The interiors and boundaries of both geometries must not intersect.

Examples

julia
import GeometryOps as GO, GeoInterface as GI

line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
point = (2, 2)
GO.disjoint(point, line)

# output
true
source
julia
disjoint(g1)

Return a function that checks if its input is disjoint from g1. This is equivalent to x -> disjoint(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.intersects Function
julia
intersects(geom1, geom2)::Bool

Return true if the interiors or boundaries of the two geometries interact.

intersects returns the exact opposite result of disjoint.

Example

julia
import GeoInterface as GI, GeometryOps as GO

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
GO.intersects(line1, line2)

# output
true
source
julia
intersects(g1)

Return a function that checks if its input intersects g1. This is equivalent to x -> intersects(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.overlaps Function
julia
overlaps(geom1, geom2)::Bool

Compare two Geometries of the same dimension and return true if their intersection set results in a geometry different from both but of the same dimension. This means one geometry cannot be within or contain the other and they cannot be equal

Examples

julia
import GeometryOps as GO, GeoInterface as GI
poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
poly2 = GI.Polygon([[(1,1), (1,6), (6,6), (6,1), (1,1)]])

GO.overlaps(poly1, poly2)
# output
true
source
julia
overlaps(g1)

Return a function that checks if its input overlaps g1. This is equivalent to x -> overlaps(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.touches Function
julia
touches([manifold::Manifold], geom1, geom2)::Bool

Return true if the first geometry touches the second geometry. In other words, the two interiors cannot interact, but one of the geometries must have a boundary point that interacts with either the other geometry's interior or boundary.

Examples

julia
import GeometryOps as GO, GeoInterface as GI

l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)])
l2 = GI.Line([(1.0, 1.0), (1.0, -1.0)])

GO.touches(l1, l2)
# output
true
source
julia
touches(g1)

Return a function that checks if its input touches g1. This is equivalent to x -> touches(x, g1).

source

This functionality is experimental and may change at any time.

source
GeometryOps.within Function
julia
within(g1)

Return a function that checks if its input is within g1. This is equivalent to x -> within(x, g1).

source

This functionality is experimental and may change at any time.

source

Other general methods ​

GeometryOps.equals Function
julia
equals(geom1, geom2)::Bool

Compare two Geometries return true if they are the same geometry.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
poly2 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])

GO.equals(poly1, poly2)
# output
true
source
julia
equals(::T, geom_a, ::T, geom_b)::Bool

Two geometries of the same type, which don't have a equals function to dispatch off of should throw an error.

source
julia
equals(trait_a, geom_a, trait_b, geom_b)

Two geometries which are not of the same type cannot be equal so they always return false.

source
julia
equals(::GI.PointTrait, p1, ::GI.PointTrait, p2)::Bool

Two points are the same if they have the same x and y (and z if 3D) coordinates.

source
julia
equals(::GI.PointTrait, p1, ::GI.MultiPointTrait, mp2)::Bool

A point and a multipoint are equal if the multipoint is composed of a single point that is equivalent to the given point.

source
julia
equals(::GI.MultiPointTrait, mp1, ::GI.PointTrait, p2)::Bool

A point and a multipoint are equal if the multipoint is composed of a single point that is equivalent to the given point.

source
julia
equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)::Bool

Two multipoints are equal if they share the same set of points.

source
julia
equals(
    ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
    ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
)::Bool

Two lines/linestrings are equal if they share the same set of points going along the curve. Note that lines/linestrings aren't closed by definition.

source
julia
equals(
    ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
    ::GI.LinearRingTrait, l2,
)::Bool

A line/linestring and a linear ring are equal if they share the same set of points going along the curve. Note that lines aren't closed by definition, but rings are, so the line must have a repeated last point to be equal

source
julia
equals(
    ::GI.LinearRingTrait, l1,
    ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
)::Bool

A linear ring and a line/linestring are equal if they share the same set of points going along the curve. Note that lines aren't closed by definition, but rings are, so the line must have a repeated last point to be equal

source
julia
equals(
    ::GI.LinearRingTrait, l1,
    ::GI.LinearRingTrait, l2,
)::Bool

Two linear rings are equal if they share the same set of points going along the curve. Note that rings are closed by definition, so they can have, but don't need, a repeated last point to be equal.

source
julia
equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool

Two polygons are equal if they share the same exterior edge and holes.

source
julia
equals(::GI.PolygonTrait, geom_a, ::GI.MultiPolygonTrait, geom_b)::Bool

A polygon and a multipolygon are equal if the multipolygon is composed of a single polygon that is equivalent to the given polygon.

source
julia
equals(::GI.MultiPolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool

A polygon and a multipolygon are equal if the multipolygon is composed of a single polygon that is equivalent to the given polygon.

source
julia
equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool

Two multipolygons are equal if they share the same set of polygons.

source

This functionality is experimental and may change at any time.

source
GeometryOps.centroid Function
julia
centroid(geom, [T=Float64])::Tuple{T, T}

Returns the centroid of a given line segment, linear ring, polygon, or mutlipolygon.

source
GeometryOps.distance Function
julia
distance(point, geom, ::Type{T} = Float64)::T

Calculates the distance from the geometry g1 to the point. The distance will always be positive or zero.

The method will differ based on the type of the geometry provided: - The distance from a point to a point is just the Euclidean distance between the points. - The distance from a point to a line is the minimum distance from the point to the closest point on the given line. - The distance from a point to a linestring is the minimum distance from the point to the closest segment of the linestring. - The distance from a point to a linear ring is the minimum distance from the point to the closest segment of the linear ring. - The distance from a point to a polygon is zero if the point is within the polygon and otherwise is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - The distance from a point to a multigeometry or a geometry collection is the minimum distance between the point and any of the sub-geometries.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.signed_distance Function
julia
signed_distance(point, geom, ::Type{T} = Float64)::T

Calculates the signed distance from the geometry geom to the given point. Points within geom have a negative signed distance, and points outside of geom have a positive signed distance. - The signed distance from a point to a point, line, linestring, or linear ring is equal to the distance between the two. - The signed distance from a point to a polygon is negative if the point is within the polygon and is positive otherwise. The value of the distance is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - The signed distance from a point to a multigeometry or a geometry collection is the minimum signed distance between the point and any of the sub-geometries.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.area Function
julia
area(geom, [T = Float64])::T
area(manifold::Manifold, geom, [T = Float64])::T
area(algorithm::Algorithm, geom, [T = Float64])::T

Returns the area of a geometry or collection of geometries. This is computed slightly differently for different geometries:

julia
- The area of a point/multipoint is always zero.
- The area of a curve/multicurve is always zero.
- The area of a polygon is the absolute value of the signed area.
- The area multi-polygon is the sum of the areas of all of the sub-polygons.
- The area of a geometry collection, feature collection of array/iterable
    is the sum of the areas of all of the sub-geometries.

Manifold support

  • AutoManifold() (default): When the Proj extension is loaded, recognized geographic CRSs use a geodesic calculation on the CRS ellipsoid and recognized projected CRSs use native-unit Planar() calculations. Without Proj, geographic geometries use degree-based Spherical() calculations, while projected and unknown geometries use native-unit Planar() calculations. The manifold is selected once from the top-level input's CRS and applies to all geometries contained in that input.

  • Planar(): Uses the shoelace formula in native coordinate units squared, regardless of CRS.

  • Spherical(): Uses Girard's theorem for spherical polygons. Coordinates are interpreted as (longitude, latitude) in degrees. Returns area in square units of the sphere's radius (default: Earth's mean radius in meters).

  • Geodesic(): Uses geodesic calculations (requires Proj extension).

Projected map area is a planar grid measurement and can differ from surface area, particularly for projections with distortion.

Examples

julia
import GeometryOps as GO
import GeoInterface as GI

# CRS-free planar area (the AutoManifold default)
rect = GI.Polygon([[(0,0), (1,0), (1,1), (0,1), (0,0)]])
GO.area(rect)  # 1.0

# Spherical area (1/8 of Earth's surface)
octant = GI.Polygon([[(0.0, 0.0), (90.0, 0.0), (0.0, 90.0), (0.0, 0.0)]])
GO.area(GO.Spherical(), octant)  # ≈ 6.38e13 m²

# Spherical area with custom radius (unit sphere)
GO.area(GO.Spherical(radius=1.0), octant)  # ≈ π/2

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.signed_area Function
julia
signed_area(geom, [T = Float64])::T

Returns the signed area of a single geometry, based on winding order. This is computed slightly differently for different geometries:

julia
- The signed area of a point is always zero.
- The signed area of a curve is always zero.
- The signed area of a polygon is computed with the shoelace formula and is
positive if the polygon coordinates wind clockwise and negative if
counterclockwise.
- You cannot compute the signed area of a multipolygon as it doesn't have a
meaning as each sub-polygon could have a different winding order.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.angles Function
julia
angles(geom, ::Type{T} = Float64)

Returns the angles of a geometry or collection of geometries. This is computed differently for different geometries:

julia
- The angles of a point is an empty vector.
- The angles of a single line segment is an empty vector.
- The angles of a linestring or linearring is a vector of angles formed by the curve.
- The angles of a polygon is a vector of vectors of angles formed by each ring.
- The angles of a multi-geometry collection is a vector of the angles of each of the
    sub-geometries as defined above.

Result will be a Vector, or nested set of vectors, of type T where an optional argument with a default value of Float64.

source
GeometryOps.embed_extent Function
julia
embed_extent(obj)

Recursively wrap the object with a GeoInterface.jl geometry, calculating and adding an Extents.Extent to all objects.

This can improve performance when extents need to be checked multiple times, such when needing to check if many points are in geometries, and using their extents as a quick filter for obviously exterior points.

Keywords

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

source

Barycentric coordinates ​

GeometryOps.barycentric_coordinates Function
julia
barycentric_coordinates(method, geom, point; normalize = true)

Return the barycentric coordinates of point with respect to the vertices of geom, as a newly allocated Vector holding one weight per vertex.

method is an AbstractBarycentricCoordinateMethod, currently always MeanValue. geom may be any GeoInterface-compatible curve (a linear ring or linestring) with at least three points, and point any GeoInterface-compatible point.

If normalize is true (the default) the returned weights are scaled to sum to 1.

See also barycentric_coordinates!, which writes into a preallocated vector, and barycentric_interpolate, which uses these weights to interpolate values.

source
GeometryOps.barycentric_coordinates! Function
julia
barycentric_coordinates!(λs, method, geom, point; normalize = true)

Write the barycentric coordinates of point with respect to the vertices of geom into λs, and return λs.

This is the allocation-free form of barycentric_coordinates; λs must have one element per vertex of geom.

source
GeometryOps.barycentric_interpolate Function
julia
barycentric_interpolate(method, geom, values, point)

Interpolate values, given one per vertex of geom, to point using barycentric coordinates, and return the interpolated value.

geom may be a curve (linear ring or linestring) or a polygon; for polygons the values are taken around the exterior ring first and then around each hole, in the order GeoInterface.getring returns them. values may be numbers or, for the mean value method, colors.

This is equivalent to weighting values by barycentric_coordinates, but accumulates in place and so allocates less.

source

Other methods ​

GeometryOps.GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS Constant

The number of vertices past which we should use a STRtree for edge intersection checking.

source
GeometryOps._VecTypes Type

native Julia vector-like types with known size

source
GeometryOps.AbstractBarycentricCoordinateMethod Type
julia
abstract type AbstractBarycentricCoordinateMethod

Abstract supertype for barycentric coordinate methods. The subtypes may serve as dispatch types, or may cache some information about the target polygon.

API

The following methods must be implemented for all subtypes:

  • barycentric_coordinates!(λs::Vector{<: Real}, method::AbstractBarycentricCoordinateMethod, exterior::Vector{<: Point{2, T1}}, point::Point{2, T2})

  • barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, exterior::Vector{<: Point{2, T1}}, values::Vector{V}, point::Point{2, T2})::V

  • barycentric_interpolate(method::AbstractBarycentricCoordinateMethod, exterior::Vector{<: Point{2, T1}}, interiors::Vector{<: Vector{<: Point{2, T1}}} values::Vector{V}, point::Point{2, T2})::V

The rest of the methods will be implemented in terms of these, and have efficient dispatches for broadcasting.

source
GeometryOps.AdjacentEdgeLocator Type
julia
AdjacentEdgeLocator(m::Manifold, geom; exact)

Determines the location for a point which is known to lie on at least one edge of a set of polygons. This provides the union-semantics for determining point location in a GeometryCollection, which may have polygons with adjacent edges which are effectively in the interior of the geometry. Note that it is also possible to have adjacent edges which lie on the boundary of the geometry (e.g. a polygon contained within another polygon with adjacent edges).

The manifold m and the exact flag are stored in the struct (rather than threaded through every call) for consistency with how RelateGeometry holds them (Task 13); locate uses the stored values for all kernel queries.

The Java constructor signature is AdjacentEdgeLocator(Geometry geom); the manifold/exact parameters are the only additions.

source
GeometryOps.AntipodalEdgeSplit Type
julia
AntipodalEdgeSplit() <: GeometryCorrection

Split every edge whose endpoints map to exactly-antipodal unit vectors by inserting the lon/lat midpoint of the edge, so each edge has a well-defined great-circle arc. This is the remedy for the antipodal-edge ArgumentError thrown by relate on the Spherical manifold.

It can be called on any geometry as usual (AntipodalEdgeSplit()(geom)), or passed to GeometryOps.fix.

See also GeometryCorrection.

source
GeometryOps.ArcCrossingCounter Type
julia
ArcCrossingCounter(m::Spherical, exact, p, anchor, 0, false)

Counts ring edges crossing the reference meridian arc from the query point p to the anchor pole, in an incremental fashion — the spherical counterpart of RayCrossingCounter. As there, the location is only correct once all edges whose longitude interval contains p's longitude have been counted, and a query point found to lie on an edge is recorded in is_point_on_segment (final location LOC_BOUNDARY).

The final location is anchor's location, flipped once per crossing. Vertex grazing — an edge endpoint exactly on the reference arc — is resolved symbolically, S2-VertexCrossing style: the edge counts iff its off-arc endpoint lies strictly on the positive side of the meridian great circle, so a crossing pair of incident edges counts once, a same-side pair counts zero or twice (parity-equal), and edges collinear with the meridian count never (their chain's terminal edges decide). With exact = True() every branch below is decided by exact predicates.

source
GeometryOps.AutoAccelerator Type
julia
AutoAccelerator()

Let the algorithm choose the best accelerator based on the size of the input polygons.

Once we have prepared geometry, this will also consider the existing preparations on the geoms.

source
GeometryOps.BoundaryNodeRule Type
julia
BoundaryNodeRule

This functionality is experimental and may change at any time.

Abstract supertype for rules deciding which endpoints of a linear geometry are on its boundary, given the number of line ends meeting at the point (port of JTS BoundaryNodeRule). Concrete rules are Mod2Boundary (the OGC SFS default), EndpointBoundary, MultivalentEndpointBoundary, and MonovalentEndpointBoundary; each implements is_in_boundary(rule, boundary_count).

source
GeometryOps.CLibraryPlanarAlgorithm Type
julia
abstract type CLibraryPlanarAlgorithm <: GeometryOpsCore.SingleManifoldAlgorithm{Planar} end

This is a type which extends GeometryOpsCore.SingleManifoldAlgorithm{Planar}, and is used as an abstract supertype for some C library based algorithms.

The type requires that algorithm structs be arranged as:

julia
struct MyAlgorithm <: CLibraryPlanarAlgorithm
    manifold::Planar
    params::NamedTuple
end

Then you get a nice constructor for free, as well as the get(alg, key, value) and get(alg, key) do ... syntax. Plus the enforce method, which will check that given keyword arguments are present.

source
GeometryOps.Chaikin Type
julia
Chaikin(; iterations=1, manifold=Planar())

Smooths geometries using Chaikin's corner-cutting algorithm [1]. This algorithm "slices" off every corner of the geometry to smooth it out, equivalent to a sequence of quadratic Bezier curves.

Keywords

  • iterations: the number of times to apply the algorithm.

  • manifold: the Manifold to smooth the geometry on. Currently, Planar and Spherical are supported.

Extended help

The algorithm is very simple; for each corner of the line (a -> b -> c), insert two new points and remove b, such that a -> b -> c becomes a -> q -> r -> c, where q and r are the new points such that:

In practice the replacement happens on the level of each edge.

References

source
GeometryOps.ClosedRing Type
julia
ClosedRing() <: GeometryCorrection

This correction ensures that a polygon's exterior and interior rings are closed.

It can be called on any geometry correction as usual.

See also GeometryCorrection.

source
GeometryOps.ConvexConvexSutherlandHodgman Type
julia
ConvexConvexSutherlandHodgman{M <: Manifold} <: GeometryOpsCore.Algorithm{M}

Sutherland-Hodgman polygon clipping algorithm optimized for convex-convex intersection.

Both input polygons MUST be convex. If either polygon is non-convex, results are undefined.

This is simpler and faster than Foster-Hormann for small convex polygons, with O(n*m) complexity where n and m are vertex counts.

When intersecting many polygon pairs in a hot loop, pass a SutherlandHodgmanCache via the cache keyword argument to intersection to avoid intermediate allocations.

Spherical manifold

For Spherical() manifold, input polygons must have counter-clockwise winding when viewed from outside the sphere (i.e., the interior is on the left when traversing edges). Polygons with clockwise winding will produce incorrect results (typically a degenerate polygon). Use GO.fix(geom; corrections=[GO.ClosedRing(), GO.GeometryCorrection()]) or manually reverse the coordinates if your input has the wrong winding order.

Example

julia
import GeometryOps as GO, GeoInterface as GI

square1 = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
square2 = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]])

result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), square1, square2)
source
GeometryOps.CrossingEdgeSplit Type
julia
CrossingEdgeSplit() <: GeometryCorrection

Split every polygon ring at the points where two of its non-adjacent edges cross properly as great-circle arcs, reassembling the resulting loops as separate rings (even-odd semantics: both lobes of a figure-eight are kept, matching S2Builder's undirected split_crossing_edges repair). A polygon whose shell splits becomes a MultiPolygon; when it carries holes, each (likewise repaired) hole loop is assigned to the shell loop that contains it. This is the remedy for the ring-crossing ArgumentError thrown by prepare on the Spherical manifold.

Crossing points are constructed in Float64 (lon/lat of the exact crossing direction) — corrections construct geometry; they don't decide predicates — the same standard as AntipodalEdgeSplit's midpoint insertion.

Scope

The correction handles isolated pairwise crossings — the needle and bowtie class, where no edge participates in more than one crossing and no two crossings interleave around the ring. Rings with tangled multi-crossing topology beyond that throw an ArgumentError rather than emitting wrong geometry. Vertex touches (rings meeting at a point) are valid and are not split.

It can be called on any polygonal geometry as usual (CrossingEdgeSplit()(geom)), or passed to GeometryOps.fix. Because a split changes the geometry type (Polygon → MultiPolygon), apply it directly to MultiPolygon inputs rather than through fix's per-polygon traversal.

See also GeometryCorrection, AntipodalEdgeSplit.

source
GeometryOps.DE9IM Type
julia
DE9IM

This functionality is experimental and may change at any time.

An immutable DE-9IM intersection matrix. Entries are dimension codes (DIM_FALSE, DIM_P, DIM_L, DIM_A) stored row-major over (Interior, Boundary, Exterior) of A × B, matching the standard string form "212101212". Construct from a 9-character string or empty (all-F) via DE9IM(). Index with im[locA, locB], where the indices are the JTS location codes (0 = Interior, 1 = Boundary, 2 = Exterior — the internal LOC_* constants), not 1-based array positions: im[0, 0] is the Interior/Interior entry.

source
GeometryOps.DiffIntersectingPolygons Type
julia
DiffIntersectingPolygons() <: GeometryCorrection

This correction ensures that the polygons included in a multipolygon aren't intersecting. If any polygon's are intersecting, they will be made nonintersecting through the difference operation to create a unique set of disjoint (other than potentially connections by a single point) polygons covering the same area. See also GeometryCorrection, UnionIntersectingPolygons.

source
GeometryOps.DouglasPeucker Type
julia
DouglasPeucker <: SimplifyAlg

DouglasPeucker(; number, ratio, tol)

Simplifies geometries by removing points below tol distance from the line between its neighboring points.

Keywords

  • ratio: the fraction of points that should remain after simplify. Useful as it will generalise for large collections of objects.

  • number: the number of points that should remain after simplify. Less useful for large collections of mixed size objects.

  • tol: the minimum distance a point will be from the line joining its neighboring points.

Note: user input tol is squared to avoid unnecessary computation in algorithm.

source
GeometryOps.EdgeSourceInfo Type
julia
EdgeSourceInfo

Source topology of one input linework (port of JTS EdgeSourceInfo): index (0 = A, 1 = B), dim (DIM_A for a ring, DIM_L for a line), is_hole (ring role), and depth_delta (the signed side-labelling delta of an area ring, 0 for a line). Consumed by the edge merger to build OverlayLabels.

source
GeometryOps.EndpointBoundary Type
julia
EndpointBoundary()

This functionality is experimental and may change at any time.

BoundaryNodeRule under which any endpoint of a line is on the boundary, regardless of how many line ends meet there (JTS EndPointBoundaryNodeRule).

source
GeometryOps.FosterHormannClipping Type
julia
FosterHormannClipping{M <: Manifold, A <: Union{Nothing, Accelerator}} <: GeometryOpsCore.Algorithm{M}

Applies the Foster-Hormann clipping algorithm.

Arguments

  • manifold::M: The manifold on which the algorithm operates.

  • accelerator::A: The accelerator to use for the algorithm. Can be nothing for automatic choice, or a custom accelerator.

source
GeometryOps.GEOS Type
julia
GEOS(; params...)

A struct which instructs the method it's passed to as an algorithm to use the appropriate GEOS function via LibGEOS.jl for the operation.

Dispatch is generally carried out using the names of the keyword arguments. For example, segmentize will only accept a GEOS struct with only a max_distance keyword, and no other.

It's generally somewhat slower than the native Julia implementations, since it must convert to the LibGEOS implementation and back - so be warned!

Extended help

This uses the LibGEOS.jl package, which is a Julia wrapper around the C library GEOS (https://trac.osgeo.org/geos).

source
GeometryOps.GeodesicSegments Type
julia
GeodesicSegments(; max_distance::Real, equatorial_radius::Real=6378137, flattening::Real=1/298.257223563)

Warning

This is deprecated - call segmentize(Geodesic(; semimajor_axis, inv_flattening), geom; max_distance) instead.

A method for segmentizing geometries by adding extra vertices to the geometry so that no segment is longer than a given distance. This method calculates the distance between points on the geodesic, and assumes input in lat/long coordinates.

Warning

Any input geometries must be in lon/lat coordinates! If not, the method may fail or error.

Arguments

  • max_distance::Real: The maximum distance, in meters, between vertices in the geometry.

  • equatorial_radius::Real=6378137: The equatorial radius of the Earth, in meters. Passed to Proj.geod_geodesic.

  • flattening::Real=1/298.257223563: The flattening of the Earth, which is the ratio of the difference between the equatorial and polar radii to the equatorial radius. Passed to Proj.geod_geodesic.

One can also omit the equatorial_radius and flattening keyword arguments, and pass a geodesic object directly to the eponymous keyword.

This method uses the Proj/GeographicLib API for geodesic calculations.

source
GeometryOps.GeometryCorrection Type
julia
abstract type GeometryCorrection

This abstract type represents a geometry correction.

Interface

Any GeometryCorrection must implement two functions: * application_level(::GeometryCorrection)::AbstractGeometryTrait: This function should return the GeoInterface trait that the correction is intended to be applied to, like PointTrait or LineStringTrait or PolygonTrait. * (::GeometryCorrection)(::AbstractGeometryTrait, geometry)::(some_geometry): This function should apply the correction to the given geometry, and return a new geometry.

source
GeometryOps.IndexedPointInAreaLocator Type
julia
IndexedPointInAreaLocator(m::Manifold, geom; exact)

Determines the location (LOC_* code) of points relative to an areal geometry, using indexing for efficiency. This algorithm is suitable for use in cases where many points will be tested against a given area. The location is computed precisely: points located on the geometry boundary or segments return LOC_BOUNDARY.

Port of JTS IndexedPointInAreaLocator together with its private IntervalIndexedGeometry (the y-interval segment index; its isEmpty flag is index === nothing here, since a recursively empty polygonal geometry contributes no rings, hence no segments). JTS lazy-loads the index on the first locate; here the index is built in the constructor, since RelatePointLocator already creates the locator itself lazily on the first use per polygonal element (_get_poly_locator, the port of RelatePointLocator.getLocator).

On Spherical, the locator caches the element's rings in kernel space (polys — SphericalKernelRings; Layer 1 of the 2026-07-14 spherical-indexed-locator design), so queries never reconvert vertices, and — with indexed = true, the default — builds the longitude-interval edge index plus the parity anchor (Layer 2): the location of a reference pole, computed once by the exact scan, from which each query is a 1-D stab at its longitude and a crossing-parity count along its meridian arc to the pole. With indexed = false (the unprepared arm, where a one-shot query cannot amortize the build) it locates by the exact ring scan over the cached rings. The polys vector is empty on Planar, where the corresponding roles are played by index and the implicit EXTERIOR at the ray's far end.

source
GeometryOps.IntersectionAccelerator Type
julia
abstract type IntersectionAccelerator

The abstract supertype for all intersection accelerator types.

The idea is that these speed up the edge-edge intersection checking process, perhaps at the cost of memory.

The naive case is NestedLoop, which is just a nested loop, running in O(n*m) time.

Then we have SingleSTRtree, which is a single STRtree, running in O(n*log(m)) time.

Then we have DoubleSTRtree, which is a simultaneous double-tree traversal of two STRtrees.

Finally, we have AutoAccelerator, which chooses the best accelerator based on the size of the input polygons. This gets materialized in build_a_list for now. AutoAccelerator should also try to respect existing spatial indices, if they exist.

source
GeometryOps.LineOrientation Type
julia
Enum LineOrientation

Enum for the orientation of a line with respect to a curve. A line can be line_cross (crossing over the curve), line_hinge (crossing the endpoint of the curve), line_over (collinear with the curve), or line_out (not interacting with the curve).

source
GeometryOps.LinearBoundary Type
julia
LinearBoundary(lines, rule::BoundaryNodeRule)

Determines the boundary points of a linear geometry, using a BoundaryNodeRule. lines is an iterable of linestrings (any GeoInterface linestring-like geometries); the endpoint degree of every line endpoint is counted and the rule decides which degrees are boundary points.

Coordinate keys are normalized via _node_point (kernel.jl): exact (Float64, Float64) tuples with signed zeros normalized (-0.0 → +0.0), so lookups here agree with the NodeKey vertex-node identity from the kernel (Task 7) under Dict bit-pattern hashing.

Faithful to Java: only empty lines are skipped. Closed lines are NOT special-cased — a closed line contributes degree 2 to its closure vertex (both endpoints coincide), which is never a boundary under the Mod-2 or monovalent rules but would be under e.g. the endpoint rule.

source
GeometryOps.LinearSegments Type
julia
LinearSegments(; max_distance::Real)

Warning

This is deprecated - call segmentize(Planar(), geom; max_distance) instead.

A method for segmentizing geometries by adding extra vertices to the geometry so that no segment is longer than a given distance.

Here, max_distance is a purely nondimensional quantity and will apply in the input space. This is to say, that if the polygon is provided in lat/lon coordinates then the max_distance will be in degrees of arc. If the polygon is provided in meters, then the max_distance will be in meters.

source
GeometryOps.MeanValue Type
julia
MeanValue() <: AbstractBarycentricCoordinateMethod

This method calculates barycentric coordinates using the mean value method.

References

source
GeometryOps.Mod2Boundary Type
julia
Mod2Boundary()

This functionality is experimental and may change at any time.

The OGC SFS standard BoundaryNodeRule (and the RelateNG default): an endpoint is on the boundary iff an odd number of line ends meet it (the "Mod-2 rule"; JTS Mod2BoundaryNodeRule).

source
GeometryOps.MonotoneChainMethod Type
julia
MonotoneChainMethod()

This is an algorithm for the convex_hull function.

Uses DelaunayTriangulation.jl to compute the convex hull. This is a pure Julia algorithm which provides an optimal Delaunay triangulation.

See also convex_hull

source
GeometryOps.MonovalentEndpointBoundary Type
julia
MonovalentEndpointBoundary()

This functionality is experimental and may change at any time.

BoundaryNodeRule under which an endpoint is on the boundary iff exactly one line end meets it (JTS MonoValentEndPointBoundaryNodeRule).

source
GeometryOps.MultivalentEndpointBoundary Type
julia
MultivalentEndpointBoundary()

This functionality is experimental and may change at any time.

BoundaryNodeRule under which an endpoint is on the boundary iff more than one line end meets it (JTS MultiValentEndPointBoundaryNodeRule).

source
GeometryOps.NodeKey Type
julia
NodeKey{P}

Symbolic identity of a node (design D2). Vertex nodes key exactly by their coordinate (is_crossing == false, all point fields equal to the vertex); proper-crossing nodes key by their canonicalized defining segment pair (is_crossing == true, fields (pt, a1) and (b0, b1) are the two segments). No intersection coordinate is ever computed for the key. Construct via vertex_node and crossing_node.

source
GeometryOps.NodeSection Type
julia
NodeSection{P}

Represents a computed node along with the incident edges on either side of it (if they exist). This captures the information about a node in a geometry component required to determine the component's contribution to the node topology. A node in an area geometry always has edges on both sides of the node. A node in a linear geometry may have one or other incident edge missing, if the node occurs at an endpoint of the line.

The edges of an area node are assumed to be provided with CW-shell orientation (as per JTS norm). This must be enforced by the caller.

Port of JTS NodeSection, with one symbolic twist (design D2): the Java class stores the node as a Coordinate; here node is a NodeKey, so proper-crossing nodes never need a constructed coordinate. v0/v1 are coordinate tuples of type P (or nothing for a missing incident edge at a line endpoint); polygonal is the parent polygonal geometry of an area section, or nothing if the section is not on a polygon boundary — an opaque payload (abstract field type, as in Java): it is only compared by identity, so parameterizing on it would re-specialize the node machinery per input geometry type for no runtime gain.

The field order matches the Java constructor argument order (isA, dimension, id, ringId, poly, isNodeAtVertex, v0, nodePt, v1). (The struct is declared before the comparator helpers below because their signatures reference it; the Java file declares EdgeAngleComparator and isAreaArea first.)

source
GeometryOps.NodeSections Type
julia
NodeSections(node::NodeKey)

Collects the NodeSections of all geometry components incident on one node, and assembles them into the node's edge topology (create_node).

Port of JTS NodeSections; the Java class is keyed by the node Coordinate, here by the symbolic NodeKey (design D2).

source
GeometryOps.NodedArrangement Type
julia
NodedArrangement{P,T}

The exactly-noded arrangement of two input geometries (design §2.1). P is the manifold's kernel point type — exactly two instantiations, Tuple{Float64,Float64} (planar) and UnitSphericalPoint{Float64} (spherical) — so the engine is type-erased over the input geometry types. T is the OUTPUT point type node_point realizes into (point_type in the OverlayNG API); it defaults to P, which is what makes an emitted vertex bit-identical to the ingested one.

Fields:

  • segstrings: the ingested inputs as RelateSegmentStrings (A side first, then B side); NodedEdge.string_idx indexes here.

  • nodes: the symbolic node table (NodeTable).

  • seg_nodes: per-parent-segment ordered interior node-id lists, keyed by (string_idx, seg_idx); absent for unsplit segments.

  • edges: every noded sub-segment of every parent segment.

  • truncated: node ids whose incident-edge set was reduced by clip pruning (see clip_a/clip_b below). Empty (BitVector()) whenever no pruning ran, which is every construction that does not pass a clip envelope.

Construct with NodedArrangement(m, a, b) (raw geometries) or NodedArrangement(m, ssa, ssb) (pre-ingested segment strings); both take point_type to choose T.

source
GeometryOps.OverlayGraph Type
julia
OverlayGraph{P,T}

The topology graph of an overlay operation (port of JTS OverlayGraph). Holds the arrangement it was built from, the vector of half-edges (both orientations of every merged edge), and one representative outgoing half-edge index per node id (node_edges[nid], 0 if the node has no incident edges). P is the manifold kernel point type, so the graph is type-erased over input geometry types; T is the arrangement's output point type, carried so the builders reading node_point off it stay concrete.

source
GeometryOps.OverlayGraph Method
julia
OverlayGraph(m, arr::NodedArrangement, sources) -> OverlayGraph

Build the overlay graph from a noded arrangement and its EdgeSourceInfo table. Coincident noded edges are merged (JTS Edge.merge semantics), each merged edge becomes a symmetric OverlayEdge pair sharing one label, and every node's star is ordered CCW about its symbolic apex via the exact kernel comparator.

source
GeometryOps.OverlayLabel Type
julia
OverlayLabel

The topological label of one edge of the overlay graph (port of JTS OverlayLabel). Mutable plain fields, one instance shared between a symmetric OverlayEdge pair; the is_forward-parameterized accessors swap Left/Right for the reverse half-edge. Index 0 selects input A, index 1 input B.

source
GeometryOps.OverlayNG Type
julia
OverlayNG(; manifold = Planar(), exact = True(), point_type = ...)
OverlayNG(manifold::Manifold; kwargs...)

This functionality is experimental and may change at any time.

The exact-arrangement overlay algorithm, a port of the JTS OverlayNG engine by Martin Davis, extended to the sphere.

OverlayNG computes intersection, union, difference and symdifference of two geometries of any dimension. It is opt-in: the algorithm is the first argument.

julia
GO.intersection(GO.OverlayNG(), a, b)
GO.union(GO.OverlayNG(GO.Spherical()), a, b)

Foster–Hormann clipping remains the default engine for intersection, union and difference when no algorithm is given; those defaults are unchanged by the presence of this algorithm. symdifference is the one exception — see its docstring.

Keyword arguments

  • manifold: Planar() (default) or Spherical(). Spherical(; oriented) selects how a ring denotes its region, and overlay follows that choice throughout — nothing extra is needed here.

  • exact: True() (default) to decide every uncertain filter with an exact predicate, False() to stay in Float64. Leave it at the default unless you are measuring the cost of exactness.

  • point_type: the type of the coordinates in the result. Defaults to the manifold's own working point type — Tuple{Float64,Float64} on Planar(), UnitSphericalPoint{Float64} (3D unit-sphere xyz) on Spherical(). On the sphere Tuple{Float64,Float64} is also accepted and gives (lon, lat) degrees; see "Output coordinates" below for what that costs.

Output coordinates

The arrangement is exact and symbolic: rounding to Float64 happens exactly once, when a node's coordinate is realized for output. point_type chooses the chart that rounding lands in, and on the sphere the choice is not neutral.

Spherical() works in unit-sphere xyz throughout — an input vertex is converted to a UnitSphericalPoint at ingest and every predicate reads it there — so emitting UnitSphericalPoint{Float64} is a pass-through for every result vertex that is an input vertex: the coordinate that comes out is bit-for-bit the one that went in, and GO.area of a clipped cell agrees with the uncut original exactly rather than in all but the last ULP.

Emitting (lon, lat) instead sends that vertex back out through atan/asin, which is a rounding of a value that had an exact image in the format it was already in. Measured over 200 000 uniformly random directions, the round trip displaces a point by up to 3.2 ULPs of the unit sphere below 60° latitude, 7.9 at 75°, and 126 above 89.5°; swept along single parallels the worst displacement is 7.8 ULPs at 89°, 1 364 at 89.99° and 2 915 (6.5e-13 rad) at 89.999°. The growth is the chart's, not the arithmetic's — a degree of longitude is cos φ of an arc — and it is why a polar grid cell survives an overlay in xyz and does not in lon/lat.

For a crossing node there is no exact Float64 answer in either chart: the position is a Rational{BigInt} direction with no finite decimal form. What xyz buys there is rounding once (normalize the direction) rather than twice (normalize, then trigonometry).

Tuple{Float64,Float64} therefore exists for callers that need lon/lat coordinates back and are willing to pay for them, not as a lossless alternative.

Result shape

The result is a single GeoInterface geometry, the most specific one that fits:

  • one component of one dimension → Point / LineString / Polygon;

  • several components of one dimension → the corresponding Multi geometry;

  • components of several dimensions → a GeometryCollection, ordered areas, then lines, then points;

  • nothing → an empty MultiPoint/MultiLineString/MultiPolygon, whose dimension follows the OGC rule for the operation (min of the inputs for intersection, max for union and symmetric difference, the left input's for difference).

Lower-dimension components are included: the intersection of two polygons that share a boundary segment and also overlap is a GeometryCollection of the overlap polygon and the shared line. This is JTS's original (non-strict) overlay semantics.

target: asking for one dimension

Each operation takes a target keyword that narrows the result to a single dimension, chosen up front rather than by the data. It is the target of the Foster–Hormann entry points, with the same return shapes — a singular trait gives the Vector of atomic components, a Multi trait gives the one multi-geometry:

targetreturns
nothing (default)as above — most specific over every dimension
GI.PolygonTrait()Vector{<:Polygon}
GI.MultiPolygonTrait()MultiPolygon
GI.LineStringTrait()Vector{<:LineString}
GI.MultiLineStringTrait()MultiLineString
GI.PointTrait()Vector{<:Point}
GI.MultiPointTrait()MultiPoint
julia
#-- always a MultiPolygon, whatever `a` and `b` turn out to share
GO.intersection(GO.OverlayNG(), a, b; target = GI.MultiPolygonTrait())

An empty targeted result is the empty Vector or empty multi-geometry of that same concrete type, so the return type no longer depends on the inputs — which is the point: without target, code that only wants areas has to handle a GeometryCollection that appears only when the inputs happen to touch.

target also removes work. A target above the result's OGC dimension is answered from the input dimensions alone, with no noding at all (an areal target on a line ∩ area is empty for every possible input), and the builds the target cannot want are skipped. The saving is one-directional, because the three builds form a dependency chain: an areal target skips both the line and point builds, a line target skips the point build, and a point target skips neither. Result polygons are built even for a line or point target — whether a line lies inside the result area is part of deciding it, and there is no cheaper equivalent test.

Inputs

Point, MultiPoint, LineString, LinearRing, MultiLineString, Polygon and MultiPolygon, in any A × B combination. GeometryCollection inputs raise an ArgumentError.

Inputs must be valid: rings may not self-cross, and the components of a multi-geometry may not overlap. The engine nodes A against B but never A against itself, so an invalid input yields an undefined result rather than an error. Fix inputs first (e.g. with CrossingEdgeSplit / AntipodalEdgeSplit on the sphere) if you are not sure.

Robustness

All topological decisions are made by exact predicates on the input coordinates and on symbolic crossing keys. Float64 values are used only as filters with certified error bounds, and only appear in the output. There is therefore no snapping, no tolerance and no precision model to configure, and no robustness failure mode to retry around.

Limitation: the full sphere

On Spherical(), an overlay whose result covers the entire sphere and has no boundary at all cannot be returned, and raises an ArgumentError. A polygon denotes the region bounded by its rings, and a polygon with no rings already means the empty geometry, so GeometryOps has no spelling for the full sphere. Reformulate such an operation as a difference from the covering region.

source
GeometryOps.PROJ Type
julia
PROJ(; params...)

A struct which instructs the method it's passed to as an algorithm to use the appropriate PROJ function via Proj.jl for the operation.

Extended help

This is the default algorithm for reproject, and will also be the default algorithm for operations on geodesics like area and arclength.

source
GeometryOps.PlanarCircle Type
julia
PlanarCircle{T}

A circle in 2D Euclidean space, represented by its center and squared radius.

Use radius(circle) to get the actual radius (computes sqrt).

Experimental

This type is not part of the public API and is subject to change without a breaking version. It implements GeoInterface's PolygonTrait, so code using GeoInterface methods will remain compatible even if the concrete type changes.

Fields

  • center::Tuple{T, T}: The (x, y) coordinates of the center

  • radius_squared::T: The squared radius (avoids sqrt in distance comparisons)

source
GeometryOps.PlanarCircleRing Type
julia
PlanarCircleRing{T}

A lazy wrapper that presents a PlanarCircle as a LinearRing with interpolated points. Used internally for GeoInterface compatibility.

Internal

This is an internal type and not part of the public API.

source
GeometryOps.PointOrientation Type
julia
Enum PointOrientation

Enum for the orientation of a point with respect to a curve. A point can be point_in the curve, point_on the curve, or point_out of the curve.

source
GeometryOps.PreparedRelate Type
julia
PreparedRelate{ALG, G, RG, SS, T}

This functionality is experimental and may change at any time.

A prepared RelateNG instance for optimized repeated evaluation of topological relationships against a single geometry a (the "prepared mode" of JTS RelateNG.prepare). Holds:

  • alg: the RelateNG algorithm configuration,

  • input: the A geometry as passed to prepare, which re-prepares from it under a different alg (geom_a below caches alg's manifold's extents, so it cannot serve another),

  • geom_a: the A-side RelateGeometry, constructed with is_prepared = true and with its lazy locator/unique-points caches forced,

  • segs_a: the A segment strings, extracted once without an interaction-envelope filter so they serve any B geometry,

  • edge_tree: the prebuilt segment index over segs_a (_relate_edge_index, edge_intersector.jl — the stand-in for Java's cached MCIndexSegmentSetMutualIntersector), or nothing below the accelerator size threshold (where the nested loop wins).

Construct with prepare; evaluate with relate / relate_predicate.

Warning

Not safe for concurrent use: self-noding evaluations mutate the held RelateGeometry (edge re-extraction, element-id counter). Use one PreparedRelate per thread.

source
GeometryOps.RadialDistance Type
julia
RadialDistance <: SimplifyAlg

Simplifies geometries by removing points less than tol distance from the line between its neighboring points.

Keywords

  • ratio: the fraction of points that should remain after simplify. Useful as it will generalise for large collections of objects.

  • number: the number of points that should remain after simplify. Less useful for large collections of mixed size objects.

  • tol: the minimum distance between points.

Note: user input tol is squared to avoid unnecessary computation in algorithm.

source
GeometryOps.RayCrossingCounter Type
julia
RayCrossingCounter(m::Manifold, p; exact)

Counts the number of segments crossed by a horizontal ray extending to the right from a given point, in an incremental fashion. This can be used to determine whether a point lies in a polygonal geometry. The class determines the situation where the point lies exactly on a segment. This handles polygonal geometries with any number of shells and holes; ring orientation is unimportant. In order to compute a correct location for a given polygonal geometry, it is essential that all segments are counted which touch the ray or lie in any ring which may contain the point — which is what allows optimization by y-interval indexing, since segments whose y-extent misses the ray a priori cannot touch it.

The manifold m and the exact flag are stored in the struct (consistent with AdjacentEdgeLocator); the orientation test goes through rk_orient (JTS uses the extended-precision Orientation.index, matching exact = True()). The horizontal-ray sweep itself is coordinate-plane logic, exactly as in JTS.

source
GeometryOps.RelateEdge Type
julia
RelateEdge{P}

An edge of a RelateNode's wheel: the direction node → dir_pt, labeled per input geometry with the dimension of the geometry element the edge came from and the geometry's location on the left of, right of, and on the edge. Unknown dimensions are DIM_UNKNOWN_EDGE; unknown locations are LOC_NONE.

Port of JTS RelateEdge; the Java class stores its parent RelateNode to reach the node coordinate, here the symbolic NodeKey is stored directly (design D2).

source
GeometryOps.RelateGeometry Type
julia
RelateGeometry(m::Manifold, geom; exact, is_prepared = false,
               boundary_rule = Mod2Boundary())

The input-geometry facade of RelateNG: wraps one of the two operand geometries and caches its metadata — recursive emptiness, extent, dimension analysis (has_points/has_lines/has_areas), zero-length-line degeneracy — plus lazily created unique points and a RelatePointLocator.

The Java constructor signature is RelateGeometry(Geometry input, boolean isPrepared, BoundaryNodeRule bnRule); the manifold/exact parameters are the only additions (consistent with RelatePointLocator). Where the Java caches geomEnv = input.getEnvelopeInternal(), here extent is the union of the interaction bounds (rk_interaction_bounds) of the non-empty elements, or nothing if the geometry is empty.

source
GeometryOps.RelateNG Type
julia
RelateNG{M <: Manifold, A <: IntersectionAccelerator, E, BR <: BoundaryNodeRule}

This functionality is experimental and may change at any time.

The next-generation DE-9IM topological-relationship algorithm, a port of the JTS RelateNG algorithm by Martin Davis. Capabilities: 2. Efficient short-circuited evaluation of topological predicates (including matching custom DE-9IM matrix patterns).

  1. Robust evaluation: all answers are computed from exact predicates on the input coordinates (no constructed intersection points), so invalid topology does not cause failures.

  2. GeometryCollection inputs containing mixed types and overlapping polygons are supported, using union semantics.

  3. Zero-length LineStrings are treated as being topologically identical to Points.

  4. Support for BoundaryNodeRules.

All coordinates are evaluated as Float64: input coordinates are converted on extraction, and the exact-predicate machinery (adaptive orientation predicates, rational-arithmetic node coincidence) assumes Float64 inputs. Non-Float64 geometries are accepted but evaluated at Float64 precision.

Keyword arguments (all optional): manifold (default Planar()), accelerator (default AutoAccelerator), exact (default True()), boundary_rule (default Mod2Boundary, the OGC SFS rule).

Unprepared performance

Every unprepared evaluation rebuilds an extent-annotated view of both inputs (the stand-in for the envelope cache JTS keeps on each Geometry), one coordinate pass per call. Inputs that already carry extents at every level skip that pass — stamp them once with

julia
geom = GO.tuples(geom; calc_extent = true)

and repeated unprepared calls read the stored extents instead. When one geometry is queried many times, prepare it instead: the prepared form also caches the point locators and edge index.

See relate and relate_predicate for the entry points.

source
GeometryOps.RelateNode Type
julia
RelateNode(m::Manifold, node::NodeKey; exact)

The topology at a node between the edges of two input geometries: a list of the RelateEdges around the node in CCW order, ordered by their CCW angle with the positive X-axis.

Port of JTS RelateNode; the Java class is keyed by the node Coordinate, here by the symbolic NodeKey (design D2). The manifold and the exact flag (absent in Java) are stored for the edge-angle comparisons in add_edge!.

source
GeometryOps.RelatePointLocator Type
julia
RelatePointLocator(m::Manifold, geom; exact, is_prepared = false,
                   boundary_rule = Mod2Boundary())

Locates a point on a geometry, including mixed-type collections. The dimension of the containing geometry element is also determined. GeometryCollections are handled with union semantics; i.e. the location of a point is that location of that point on the union of the elements of the collection.

Union semantics for GeometryCollections has the following behaviours: 2. For a mixed-dimension (heterogeneous) collection a point may lie on two geometry elements with different dimensions. In this case the location on the largest-dimension element is reported.

  1. For a collection with overlapping or adjacent polygons, points on polygon element boundaries may lie in the effective interior of the collection geometry.

Supports specifying the BoundaryNodeRule to use for line endpoints (RelateGeometry passes its rule down here; the default matches Java's BoundaryNodeRule.OGC_SFS_BOUNDARY_RULE, i.e. Mod-2).

The Java constructor signature is RelatePointLocator(geom, isPrepared, bnRule); the manifold/exact parameters are the only additions (consistent with AdjacentEdgeLocator). As in JTS, prepared mode swaps the per-polygon SimplePointInAreaLocator ring loop for a cached IndexedPointInAreaLocator (indexed_point_in_area.jl), created lazily on the first use per polygonal element (Task 22); unprepared planar mode scans the rings directly on every query, while the spherical path caches a per-element locator in both modes (its kernel-space ring cache is what makes queries conversion-free). Repeated point location against one geometry is what prepare is for.

source
GeometryOps.RelateSegmentString Type
julia
RelateSegmentString{P}

Models a linear edge of a RelateGeometry: the coordinate vector of one line or one polygon ring, tagged with which input geometry it came from (is_a), its dimension, the element/ring ids assigned during extraction, and (for rings) the parent polygonal geometry.

In JTS this extends BasicSegmentString; here the coordinates are stored directly in pts. Segment indices are 1-based: segment i runs from pts[i] to pts[i + 1] (the Java equivalents are 0-based).

The geometry references are deliberately opaque (abstract field types, as they are in Java): parent_polygonal is only ever compared by identity and carried into NodeSections, so parameterizing on the input geometry type would only re-specialize the whole edge machinery per geometry-type pair (a pure compile-time cost). Only pts — the per-segment hot path — stays concretely typed, on the manifold's kernel point type P.

source
GeometryOps.SegSegClass Type
julia
SegSegClass

Combinatorial classification of the intersection of closed segments (a0,a1) × (b0,b1). kind is SS_PROPER only for a crossing in both segments' interiors (the node is symbolic: no coordinate exists for it anywhere in the engine). All vertex incidences are reported via the *_on_* flags, whose coordinates are exact input vertices.

source
GeometryOps.SimplifyAlg Type
julia
abstract type SimplifyAlg

Abstract type for simplification algorithms.

API

For now, the algorithm must hold the number, ratio and tol properties.

Simplification algorithm types can hook into the interface by implementing the _simplify(trait, alg, geom) methods for whichever traits are necessary.

source
GeometryOps.SphericalKernelRing Type
julia
SphericalKernelRing(m::Spherical, ring; exact, is_hole = false)

The cached kernel-space form of one ring: the converted UnitSphericalPoint vertex vector (pts — the boundary edge walk), its deduped open form (ded/n — the parity walk; aliases pts when the ring has no repeated vertices), the ring's denoted-region bit (_ring_interior_on_left, from the ring's winding or — on an oriented manifold — its declared role; the same bit edge topology and interaction bounds use), and — in enclosed-region mode — the definitional-exterior parity anchor (spherical_exterior_anchor; nothing on an oriented manifold, which never consults it, or for a degenerate vertex mass, where queries fall back to the wedge bootstrap).

rk_point_in_ring re-derived all of this from lon/lat on every query — vertex conversion alone was ~60% of a prepared spherical point query. The point-in-area locators (indexed_point_in_area.jl) convert each ring once and query on this form (Layer 1 of the 2026-07-14 spherical-indexed-locator design).

Repeated consecutive vertices are dropped from the parity walk (real rings carry them — NE 110m North Korea's sliver is [A, A, B, A]; JTS removes them at ingest, but this path receives the raw ring): a retraced edge lies exactly under the anchor midpoint and breaks the parity count. After dedup a ring with fewer than 3 distinct vertices bounds no area.

source
GeometryOps.SphericalRingPoints Type
julia
SphericalRingPoints(ring)

A ring's vertices as UnitSphericalPoints, converted on indexing rather than up front, with any repeated closing vertex excluded from length.

Exists so that the spherical predicates can hand a ring to the UnitSpherical ring primitives — which index into a vector of unit points — without allocating that vector. Converting lazily costs one UnitSphereFromGeographic per access; a ring walked several times in one predicate call pays that more than once, which is the trade a prepared target would remove.

source
GeometryOps.SutherlandHodgmanCache Type
julia
SutherlandHodgmanCache{P}()
SutherlandHodgmanCache(manifold::Manifold, [T = Float64])
SutherlandHodgmanCache(alg::ConvexConvexSutherlandHodgman, [T = Float64])

Preallocated buffers for ConvexConvexSutherlandHodgman clipping.

Pass this as the cache keyword argument to intersection to avoid allocating intermediate vectors on every call - useful when intersecting many polygon pairs in a hot loop. The returned polygon never aliases the cache, so results remain valid after the cache is reused.

The point type P must match the algorithm's manifold and float type: Tuple{T,T} for Planar(), UnitSpherical.UnitSphericalPoint{T} for Spherical(). The manifold/algorithm constructors take care of this.

Thread safety

A cache must not be shared across concurrent tasks or threads. Create one cache per task. The default (cache = nothing) allocates fresh buffers on each call and is always safe.

Example

julia
import GeometryOps as GO

alg = GO.ConvexConvexSutherlandHodgman()
cache = GO.SutherlandHodgmanCache(alg)
for (a, b) in polygon_pairs
    result = GO.intersection(alg, a, b; cache)
end
source
GeometryOps.TG Type
julia
TG(; params...)

A struct which instructs the method it's passed to as an algorithm to use the appropriate TG function via TGGeometry.jl for the operation.

It's generally a lot faster than the native Julia implementations, but only supports planar manifolds / operations. Also, it only supports geometric predicates, specifically the ones which the underlying tg library supports. These are:

equals, intersects, disjoint, contains, within, covers, coveredby, and touches.

Extended help

This uses the TGGeometry.jl package, which is a Julia wrapper around the tg C library (https://github.com/tidwall/tg).

source
GeometryOps.TopologyComputer Type
julia
TopologyComputer(predicate, geom_a::RelateGeometry, geom_b::RelateGeometry)

The DE-9IM accumulation engine of RelateNG: translates topological events into dimension updates on predicate and collects edge-intersection NodeSections per node (keyed by symbolic NodeKey, design D2) for evaluate_nodes!.

Port of JTS TopologyComputer.

source
GeometryOps.TopologyPredicate Type
julia
TopologyPredicate

The abstract supertype for strategy types implementing spatial predicates based on the DE-9IM topology model (port of JTS TopologyPredicate). Concrete predicates implement predicate_name(p), update_dim!(p, locA, locB, dim), finish!(p), is_known(p), and predicate_value(p), and may override the requirement flags and the init_dims!/init_bounds! hooks. Evaluate one against a pair of geometries with relate_predicate.

source
GeometryOps.TracingError Type
julia
TracingError{T1, T2} <: Exception

An error that is thrown when the clipping tracing algorithm fails somehow. This is a bug in the algorithm, and should be reported.

The polygons are contained in the exception object, accessible by try-catch or as err in the REPL.

source
GeometryOps.UnionIntersectingPolygons Type
julia
UnionIntersectingPolygons() <: GeometryCorrection

This correction ensures that the polygon's included in a multipolygon aren't intersecting. If any polygon's are intersecting, they will be combined through the union operation to create a unique set of disjoint (other than potentially connections by a single point) polygons covering the same area.

See also GeometryCorrection.

source
GeometryOps.VisvalingamWhyatt Type
julia
VisvalingamWhyatt <: SimplifyAlg

VisvalingamWhyatt(; kw...)

Simplifies geometries by removing points below tol distance from the line between its neighboring points.

Keywords

  • ratio: the fraction of points that should remain after simplify. Useful as it will generalise for large collections of objects.

  • number: the number of points that should remain after simplify. Less useful for large collections of mixed size objects.

  • tol: the minimum area of a triangle made with a point and its neighboring points.

Note: user input tol is doubled to avoid unnecessary computation in algorithm.

source
GeometryOps.Welzl Type
julia
Welzl{M <: Manifold} <: ManifoldIndependentAlgorithm{M}

Welzl's algorithm for computing the minimum bounding circle.

This is a randomized algorithm with expected O(n) time complexity. Works on any manifold given an appropriate distance function.

Constructor

julia
Welzl(; manifold=Planar())

Example

julia
import GeometryOps as GO

points = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
circle = GO.minimum_bounding_circle(GO.Welzl(), points)
source
GeometryOps._OverlayTopologyError Type
julia
_OverlayTopologyError(msg)

A robustness/topology error raised by the overlay engine (port of JTS TopologyException). Signals an inconsistency the builder could not resolve (e.g. a side-location conflict during area propagation, or a ring that cannot be closed).

source
Extents.extent Method
julia
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, whose extent also covers any enclosed axis point (a pole, say). Which region a ring bounds follows the manifold's interior mode (see Spherical): by default the region it encloses, independent of winding; with Spherical(; oriented = true) the region on the left of the stored vertex order, so a clockwise ring denotes the complement (and gets a box covering essentially the whole sphere).

Example

julia
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
source
GeometryOps._det Method
julia
_det(s1::Point2{T1}, s2::Point2{T2}) where {T1 <: Real, T2 <: Real}

Returns the determinant of the matrix formed by hcat'ing two points s1 and s2.

Specifically, this is:

julia
s1[1] * s2[2] - s1[2] * s2[1]
source
GeometryOps._equals_curves Method
julia
_equals_curves(c1, c2, closed_type1, closed_type2)::Bool

Two curves are equal if they share the same set of point, representing the same geometry. Both curves must must be composed of the same set of points, however, they do not have to wind in the same direction, or start on the same point to be equivalent. Inputs: c1 first geometry c2 second geometry closed_type1::Bool true if c1 is closed by definition (polygon, linear ring) closed_type2::Bool true if c2 is closed by definition (polygon, linear ring)

source
GeometryOps._exact_crossing_point Method

Exact intersection point of two properly crossing segments, as rationals.

source
GeometryOps._is_canonical_incidence Method
julia
_is_canonical_incidence(ss::RelateSegmentString, seg_index::Integer, pt)

The once-only rule for vertex incidences: whether segment seg_index of ss is the canonical owner of the intersection point pt. Segments are half-closed — a segment owns its start vertex but not its end vertex, except for the final segment of a non-closed string, which also owns its endpoint; in a closed ring the closing vertex is owned by the first segment (wraparound). This attributes every vertex of a segment string to exactly one of its segments, so an endpoint intersection enumerated against both incident segments produces node sections once, not twice.

Encodes the canonicality semantics of the Java addIntersections/RelateSegmentString.isContainingSegment pairing; delegates to is_containing_segment (its direct port).

source
GeometryOps._overlaps Method
julia
_overlaps(::GI.AbstractTrait, geom1, ::GI.AbstractTrait, geom2)::Bool

For any non-specified pair, all have non-matching dimensions, return false.

source
GeometryOps._overlaps Method
julia
_overlaps(::GI.LineTrait, line1, ::GI.LineTrait, line2)::Bool

If the lines overlap, meaning that they are collinear but each have one endpoint outside of the other line, return true. Else false.

source
GeometryOps._overlaps Method
julia
_overlaps(
    ::GI.MultiPointTrait, points1,
    ::GI.MultiPointTrait, points2,
)::Bool

If the multipoints overlap, meaning some, but not all, of the points within the multipoints are shared, return true.

source
GeometryOps._overlaps Method
julia
_overlaps(
    ::GI.MultiPolygonTrait, polys1,
    ::GI.MultiPolygonTrait, polys2,
)::Bool

Return true if at least one pair of polygons from multipolygons overlap. Else false.

source
GeometryOps._overlaps Method
julia
_overlaps(
    ::GI.MultiPolygonTrait, polys1,
    ::GI.PolygonTrait, poly2,
)::Bool

Return true if polygon overlaps with at least one of the polygons within the multipolygon. Else false.

source
GeometryOps._overlaps Method
julia
_overlaps(
    ::GI.PolygonTrait, poly1,
    ::GI.MultiPolygonTrait, polys2,
)::Bool

Return true if polygon overlaps with at least one of the polygons within the multipolygon. Else false.

source
GeometryOps._overlaps Method
julia
_overlaps(
    trait_a::GI.PolygonTrait, poly_a,
    trait_b::GI.PolygonTrait, poly_b,
)::Bool

If the two polygons intersect with one another, but are not equal, return true. Else false.

source
GeometryOps._overlaps Method
julia
_overlaps(
    ::Union{GI.LineStringTrait, GI.LinearRing}, line1,
    ::Union{GI.LineStringTrait, GI.LinearRing}, line2,
)::Bool

If the curves overlap, meaning that at least one edge of each curve overlaps, return true. Else false.

source
GeometryOps._overlay_ng Method
julia
_overlay_ng(m, op::_OverlayOpCode, a, b; exact=True(), tree_a=nothing, tree_b=nothing, target=nothing, point_type=_kernel_point_type(m))

Compute the overlay of a and b under op on manifold m, returning a GeoInterface geometry. Internal engine entry point for OverlayNG — point, line and area inputs are supported in any A×B combination. tree_a/tree_b accept caller-prebuilt segment indices (threaded to the noding substrate). target narrows the result to one dimension — see "Result targeting" below. point_type is the output coordinate type (OverlayNG's keyword of the same name).

point_type becomes a POSITIONAL ::Type{T} one call in, for the reason NodedArrangement does the same: a keyword-bound static parameter does not specialize, and every result type in this file is a function of T.

source
GeometryOps.add_area_vertex! Method
julia
add_area_vertex!(tc, is_area_a, loc_area, loc_target, dim_target, pt)

Adds topology for an area vertex interaction with a target geometry element. Assumes the target geometry element has highest dimension (i.e. if the point lies on two elements of different dimension, the location on the higher dimension element is provided. This is the semantic provided by RelatePointLocator.)

Note that in a GeometryCollection containing overlapping or adjacent polygons, the area vertex location may be INTERIOR instead of BOUNDARY.

Port of TopologyComputer.addAreaVertex.

source
GeometryOps.add_edge! Method
julia
add_edge!(n::RelateNode, is_a, dir_pt, dim, is_forward)

Adds or merges an edge to the node, keeping the wheel sorted by CCW angle with the positive X-axis. dim is the dimension of the geometry element containing the edge, is_forward the direction of the edge. Returns the created or merged edge for this point, or nothing for a malformed (nothing or zero-length) input edge.

Port of RelateNode.addEdge.

source
GeometryOps.add_intersections! Method
julia
add_intersections!(tc::TopologyComputer, ssA, seg_index_a, ssB, seg_index_b; m, exact)

Classify the intersection of one segment pair via rk_classify_intersection and add a NodeSection pair (one on ssA, one on ssB) to tc for each distinct intersection point:

  • SS_DISJOINT: nothing to add.

  • SS_PROPER: one section pair at the symbolic crossing_node.

  • SS_TOUCH/SS_COLLINEAR: a section pair at the vertex_node of each distinct flagged vertex — but only when both segments contain the vertex canonically (_is_canonical_incidence), which ensures endpoint intersections are added once only across adjacent segments.

Port of EdgeSegmentIntersector.addIntersections (private).

source
GeometryOps.add_line_end_on_geometry! Method
julia
add_line_end_on_geometry!(tc, is_line_a, loc_line_end, loc_target, dim_target, pt)

Add topology for a line end. The line end point must be "significant"; i.e. not contained in an area if the source is a mixed-dimension GC. loc_line_end is the location of the line end (Interior or Boundary); loc_target the location on the target geometry; dim_target the dimension of the interacting target geometry element (if any), or the dimension of the target.

Port of TopologyComputer.addLineEndOnGeometry.

source
GeometryOps.angles Method
julia
angles(geom, ::Type{T} = Float64)

Returns the angles of a geometry or collection of geometries. This is computed differently for different geometries:

julia
- The angles of a point is an empty vector.
- The angles of a single line segment is an empty vector.
- The angles of a linestring or linearring is a vector of angles formed by the curve.
- The angles of a polygon is a vector of vectors of angles formed by each ring.
- The angles of a multi-geometry collection is a vector of the angles of each of the
    sub-geometries as defined above.

Result will be a Vector, or nested set of vectors, of type T where an optional argument with a default value of Float64.

source
GeometryOps.antimeridian_split Method
julia
antimeridian_split(geom; antimeridian = 180.0, north_pole = nothing, pole_spacing = 5.0)

Split a lon/lat Polygon or MultiPolygon at the antimeridian, returning a GI.MultiPolygon whose pieces each stay within one 360°-wide longitude branch (none crosses the seam). The spherical region is preserved exactly — this is an encoding repair, not a geometry change.

Keyword arguments

  • antimeridian = 180.0: the seam longitude λ. The cut runs along the meridian λ (normalised to λn ∈ (-180, 180]); the default reproduces the ±180° antimeridian. Emitted longitudes lie in the closed branch [λn - 360, λn]: the two seam lips are exactly λn (west-side pieces) and λn - 360 (east-side pieces), and every non-seam vertex is strictly interior, in (λn - 360, λn). At the default seam this is the usual [-180, 180] with lips at +180 / −180.

  • north_pole = nothing: when set to (λp, φp), cut along a seam through the rotated pole at geographic (λp, φp) (CF rotated_latitude_longitude / +proj=ob_tran, north_pole_grid_longitude = 0). The returned coordinates are in the rotated frame, not geographic — this is deliberate.

  • pole_spacing = 5.0: maximum longitude step (degrees) of the constant-latitude row emitted along a pole edge of a pole-enclosing piece. nothing emits only the two branch corners. The corners are never optional — they are the topological product of the face walk; only the infill between them is controlled here.

Only PolygonTrait and MultiPolygonTrait inputs are supported; other traits throw an ArgumentError (LineString support is future work).

source
GeometryOps.area Method
julia
area(geom, [T = Float64])::T
area(manifold::Manifold, geom, [T = Float64])::T
area(algorithm::Algorithm, geom, [T = Float64])::T

Returns the area of a geometry or collection of geometries. This is computed slightly differently for different geometries:

julia
- The area of a point/multipoint is always zero.
- The area of a curve/multicurve is always zero.
- The area of a polygon is the absolute value of the signed area.
- The area multi-polygon is the sum of the areas of all of the sub-polygons.
- The area of a geometry collection, feature collection of array/iterable
    is the sum of the areas of all of the sub-geometries.

Manifold support

  • AutoManifold() (default): When the Proj extension is loaded, recognized geographic CRSs use a geodesic calculation on the CRS ellipsoid and recognized projected CRSs use native-unit Planar() calculations. Without Proj, geographic geometries use degree-based Spherical() calculations, while projected and unknown geometries use native-unit Planar() calculations. The manifold is selected once from the top-level input's CRS and applies to all geometries contained in that input.

  • Planar(): Uses the shoelace formula in native coordinate units squared, regardless of CRS.

  • Spherical(): Uses Girard's theorem for spherical polygons. Coordinates are interpreted as (longitude, latitude) in degrees. Returns area in square units of the sphere's radius (default: Earth's mean radius in meters).

  • Geodesic(): Uses geodesic calculations (requires Proj extension).

Projected map area is a planar grid measurement and can differ from surface area, particularly for projections with distortion.

Examples

julia
import GeometryOps as GO
import GeoInterface as GI

# CRS-free planar area (the AutoManifold default)
rect = GI.Polygon([[(0,0), (1,0), (1,1), (0,1), (0,0)]])
GO.area(rect)  # 1.0

# Spherical area (1/8 of Earth's surface)
octant = GI.Polygon([[(0.0, 0.0), (90.0, 0.0), (0.0, 90.0), (0.0, 0.0)]])
GO.area(GO.Spherical(), octant)  # ≈ 6.38e13 m²

# Spherical area with custom radius (unit sphere)
GO.area(GO.Spherical(radius=1.0), octant)  # ≈ π/2

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.barycentric_coordinates! Method
julia
barycentric_coordinates!(λs, method, geom, point; normalize = true)

Write the barycentric coordinates of point with respect to the vertices of geom into λs, and return λs.

This is the allocation-free form of barycentric_coordinates; λs must have one element per vertex of geom.

source
GeometryOps.barycentric_coordinates Method
julia
barycentric_coordinates(method, geom, point; normalize = true)

Return the barycentric coordinates of point with respect to the vertices of geom, as a newly allocated Vector holding one weight per vertex.

method is an AbstractBarycentricCoordinateMethod, currently always MeanValue. geom may be any GeoInterface-compatible curve (a linear ring or linestring) with at least three points, and point any GeoInterface-compatible point.

If normalize is true (the default) the returned weights are scaled to sum to 1.

See also barycentric_coordinates!, which writes into a preallocated vector, and barycentric_interpolate, which uses these weights to interpolate values.

source
GeometryOps.barycentric_interpolate Method
julia
barycentric_interpolate(method, geom, values, point)

Interpolate values, given one per vertex of geom, to point using barycentric coordinates, and return the interpolated value.

geom may be a curve (linear ring or linestring) or a polygon; for polygons the values are taken around the exterior ring first and then around each hole, in the order GeoInterface.getring returns them. values may be numbers or, for the mean value method, colors.

This is equivalent to weighting values by barycentric_coordinates, but accumulates in place and so allocates less.

source
GeometryOps.centroid Method
julia
centroid(geom, [T=Float64])::Tuple{T, T}

Returns the centroid of a given line segment, linear ring, polygon, or mutlipolygon.

source
GeometryOps.centroid_and_area Method
julia
centroid_and_area(geom, [T=Float64])::(::Tuple{T, T}, ::Real)

Returns the centroid and area of a given geometry.

source
GeometryOps.centroid_and_length Method
julia
centroid_and_length(geom, [T=Float64])::(::Tuple{T, T}, ::Real)

Returns the centroid and length of a given line/ring. Note this is only valid for line strings and linear rings.

source
GeometryOps.compare_to Method
julia
compare_to(ns::NodeSection, other::NodeSection)

Compare node sections by parent geometry, dimension, element id and ring id, and edge vertices. Sections are assumed to be at the same node point. Returns a negative/zero/positive Int (Java compareTo contract).

Port of NodeSection.compareTo.

source
GeometryOps.contains Method
julia
contains(g1::AbstractGeometry, g2::AbstractGeometry)::Bool

Return true if the second geometry is completely contained by the first geometry. The interiors of both geometries must intersect and the interior and boundary of the secondary (g2) must not intersect the exterior of the first (g1).

contains returns the exact opposite result of within.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
point = GI.Point((1, 2))

GO.contains(line, point)
# output
true
source
GeometryOps.contains Method
julia
contains(g1)

Return a function that checks if its input contains g1. This is equivalent to x -> contains(x, g1).

source
GeometryOps.contains Method

This functionality is experimental and may change at any time.

source
GeometryOps.convex_hull Function
julia
convex_hull([method], geometries)

Compute the convex hull of the points in geometries. Returns a GI.Polygon representing the convex hull.

Note that the polygon returned is wound counterclockwise as in the Simple Features standard by default. If you choose GEOS, the winding order will be inverted.

Warning

This interface only computes the 2-dimensional convex hull!

For higher dimensional hulls, use the relevant package (Qhull.jl, Quickhull.jl, or similar).

source
GeometryOps.count_segment! Method
julia
count_segment!(rcc::RayCrossingCounter, p1, p2)

Counts a segment with endpoints p1, p2. Port of RayCrossingCounter.countSegment.

source
GeometryOps.coverage Method
julia
coverage(geom, xmin, xmax, ymin, ymax, [T = Float64])::T

Returns the area of intersection between given geometry and grid cell defined by its minimum and maximum x and y-values. This is computed differently for different geometries:

  • The signed area of a point is always zero.

  • The signed area of a curve is always zero.

  • The signed area of a polygon is calculated by tracing along its edges and switching to the cell edges if needed.

  • The coverage of a geometry collection, multi-geometry, feature collection of array/iterable is the sum of the coverages of all of the sub-geometries.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.coveredby Method
julia
coveredby([manifold::Manifold], g1, g2)::Bool

Return true if the first geometry is completely covered by the second geometry. The interior and boundary of the primary geometry (g1) must not intersect the exterior of the secondary geometry (g2).

Furthermore, coveredby returns the exact opposite result of covers. They are equivalent with the order of the arguments swapped.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
p1 = GI.Point(0.0, 0.0)
p2 = GI.Point(1.0, 1.0)
l1 = GI.Line([p1, p2])

GO.coveredby(p1, l1)
# output
true
source
GeometryOps.coveredby Method
julia
coveredby(g1)

Return a function that checks if its input is covered by g1. This is equivalent to x -> coveredby(x, g1).

source
GeometryOps.coveredby Method

This functionality is experimental and may change at any time.

source
GeometryOps.covers Method
julia
covers(g1::AbstractGeometry, g2::AbstractGeometry)::Bool

Return true if the first geometry is completely covers the second geometry, The exterior and boundary of the second geometry must not be outside of the interior and boundary of the first geometry. However, the interiors need not intersect.

covers returns the exact opposite result of coveredby.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
l1 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)])
l2 = GI.LineString([(1.0, 1.0), (1.0, 2.0)])

GO.covers(l1, l2)
# output
true
source
GeometryOps.covers Method
julia
covers(g1)

Return a function that checks if its input covers g1. This is equivalent to x -> covers(x, g1).

source
GeometryOps.covers Method

This functionality is experimental and may change at any time.

source
GeometryOps.create_node Method
julia
create_node(m::Manifold, nss::NodeSections; exact)

Creates the node topology: prepares the sections, builds a RelateNode at the node and feeds it the sections via add_edges!. Per-polygon section groups are first rewritten into maximal-ring structure by polygon_node_convert (the PolygonNodeConverter.convert port). Returns the assembled node.

Port of NodeSections.createNode. The manifold/exact parameters (absent in Java, where createNode() is nullary) are threaded through for the angle comparisons in the converter and the node's edge wheel. (RelateNode is defined in relate_node.jl, included after this file; the reference resolves at call time.)

source
GeometryOps.create_node_section Method
julia
create_node_section(ss::RelateSegmentString, seg_index::Integer, node::NodeKey)

The NodeSection of this segment string at the node node, known to lie on segment seg_index.

Port of RelateSegmentString.createNodeSection, with the symbolic twist (design D2): the Java method takes the intersection Coordinate; here the node is identified by its NodeKey. For a vertex node (SS_TOUCH or other vertex incidences) the key carries the exact coordinate and the incident vertices are found as in Java. A proper-crossing node (SS_PROPER) lies strictly inside the segment, so the incident vertices are the segment endpoints and the node is never at a vertex.

source
GeometryOps.crosses Method
julia
 crosses(geom1, geom2)::Bool

Return true if the intersection results in a geometry whose dimension is one less than the maximum dimension of the two source geometries and the intersection set is interior to both source geometries.

TODO: broken

Examples

julia
import GeoInterface as GI, GeometryOps as GO
# TODO: Add working example
source
GeometryOps.crosses Method
julia
crosses(g1)

Return a function that checks if its input crosses g1. This is equivalent to x -> crosses(x, g1).

source
GeometryOps.crosses Method

This functionality is experimental and may change at any time.

source
GeometryOps.crossing_node Method
julia
crossing_node(a0, a1, b0, b1)::NodeKey

Node key of the proper crossing of segments (a0, a1) and (b0, b1). Canonicalize: each segment ordered lexicographically by (x, y); segments ordered lexicographically by their endpoint tuples — so any order/orientation of the same pair produces an identical key.

Only construct crossing keys for properly crossing segments (SS_PROPER from rk_classify_intersection): the exact rational slow path in rk_nodes_coincide divides by the segments' direction cross product, which is nonzero precisely when the crossing is proper.

source
GeometryOps.cut Method
julia
cut(geom, line, [T::Type])

Return given geom cut by given line as a list of geometries of the same type as the input geom. Return the original geometry as only list element if none are found. Line must cut fully through given geometry or the original geometry will be returned.

Note: This currently doesn't work for degenerate cases there line crosses through vertices.

Example

julia
import GeoInterface as GI, GeometryOps as GO

poly = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]])
line = GI.Line([(5.0, -5.0), (5.0, 15.0)])
cut_polys = GO.cut(poly, line)
GI.coordinates.(cut_polys)

# output
2-element Vector{Vector{Vector{Vector{Float64}}}}:
 [[[0.0, 0.0], [5.0, 0.0], [5.0, 10.0], [0.0, 10.0], [0.0, 0.0]]]
 [[[5.0, 0.0], [10.0, 0.0], [10.0, 10.0], [5.0, 10.0], [5.0, 0.0]]]
source
GeometryOps.difference Method

This functionality is experimental and may change at any time.

source
GeometryOps.difference Method
julia
difference(geom_a, geom_b, [T::Type]; target::Type, fix_multipoly = UnionIntersectingPolygons())

Return the difference between two geometries as a list of geometries. Return an empty list if none are found. The type of the list will be constrained as much as possible given the input geometries. Furthermore, the user can provide a taget type as a keyword argument and a list of target geometries found in the difference will be returned. The user can also provide a float type that they would like the points of returned geometries to be. If the user is taking a intersection involving one or more multipolygons, and the multipolygon might be comprised of polygons that intersect, if fix_multipoly is set to an IntersectingPolygons correction (the default is UnionIntersectingPolygons()), then the needed multipolygons will be fixed to be valid before performing the intersection to ensure a correct answer. Only set fix_multipoly to false if you know that the multipolygons are valid, as it will avoid unneeded computation.

Example

julia
import GeoInterface as GI, GeometryOps as GO

poly1 = GI.Polygon([[[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]]])
poly2 = GI.Polygon([[[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]]])
diff_poly = GO.difference(poly1, poly2; target = GI.PolygonTrait())
GI.coordinates.(diff_poly)

# output
1-element Vector{Vector{Vector{Vector{Float64}}}}:
 [[[6.5, 3.5], [5.0, 5.0], [0.0, 0.0], [5.0, -5.0], [6.5, -3.5], [3.0, 0.0], [6.5, 3.5]]]
source
GeometryOps.disjoint Method
julia
disjoint([manifold::Manifold], geom1, geom2)::Bool

Return true if the first geometry is disjoint from the second geometry.

Return true if the first geometry is disjoint from the second geometry. The interiors and boundaries of both geometries must not intersect.

Examples

julia
import GeometryOps as GO, GeoInterface as GI

line = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
point = (2, 2)
GO.disjoint(point, line)

# output
true
source
GeometryOps.disjoint Method
julia
disjoint(g1)

Return a function that checks if its input is disjoint from g1. This is equivalent to x -> disjoint(x, g1).

source
GeometryOps.disjoint Method

This functionality is experimental and may change at any time.

source
GeometryOps.distance Method
julia
distance(point, geom, ::Type{T} = Float64)::T

Calculates the distance from the geometry g1 to the point. The distance will always be positive or zero.

The method will differ based on the type of the geometry provided: - The distance from a point to a point is just the Euclidean distance between the points. - The distance from a point to a line is the minimum distance from the point to the closest point on the given line. - The distance from a point to a linestring is the minimum distance from the point to the closest segment of the linestring. - The distance from a point to a linear ring is the minimum distance from the point to the closest segment of the linear ring. - The distance from a point to a polygon is zero if the point is within the polygon and otherwise is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - The distance from a point to a multigeometry or a geometry collection is the minimum distance between the point and any of the sub-geometries.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.eachedge Method
julia
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.

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 UnitSphericalPoints: geographic (longitude, latitude) input is converted, UnitSphericalPoints pass through.

source
GeometryOps.edge_angle_compare Method
julia
edge_angle_compare(m::Manifold, ns1::NodeSection, ns2::NodeSection; exact)

Compares sections by the angle the entering edge (get_vertex(ns, 0)) makes with the positive X axis at the node, angles increasing CCW.

Port of NodeSection.EdgeAngleComparator (a static nested Comparator class in Java, compareAngle(ns1.nodePt, ns1.getVertex(0), ns2.getVertex(0))); here a comparator function over rk_compare_edge_dir with the symbolic node of ns1 as apex, taking the manifold and exact flag the kernel comparison needs. Use as a sort predicate via lt = (a, b) -> edge_angle_compare(m, a, b; exact) < 0.

At a crossing node (a NodeKey with is_crossing) the sections' v0 are normally among the four endpoints of the node's defining segments (the sections are built from those segments themselves), where the comparison is derived from the original endpoints. Sections merged onto the node by the D3 coincidence pass (TopologyComputer) may carry foreign directions; rk_compare_edge_dir then compares around the exact rational apex.

source
GeometryOps.edge_extents Method
julia
edge_extents(geom, [::Type{T}])

Return a vector of the extents of the edges (line segments) of geom.

source
GeometryOps.embed_extent Method
julia
embed_extent(obj)

Recursively wrap the object with a GeoInterface.jl geometry, calculating and adding an Extents.Extent to all objects.

This can improve performance when extents need to be checked multiple times, such when needing to check if many points are in geometries, and using their extents as a quick filter for obviously exterior points.

Keywords

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

source
GeometryOps.enforce Method
julia
enforce(alg::CLibraryPlanarAlgorithm, kw::Symbol, f)

Enforce the presence of a keyword argument in a GEOS algorithm, and return alg.params[kw].

Throws an error if the key is not present, and mentions f in the error message (since there isn't a good way to get the name of the function that called this method).

This applies to all CLibraryPlanarAlgorithm types, like GEOS and TG.

source
GeometryOps.equals Method
julia
equals(trait_a, geom_a, trait_b, geom_b)

Two geometries which are not of the same type cannot be equal so they always return false.

source
GeometryOps.equals Method
julia
equals(geom1, geom2)::Bool

Compare two Geometries return true if they are the same geometry.

Examples

julia
import GeometryOps as GO, GeoInterface as GI
poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
poly2 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])

GO.equals(poly1, poly2)
# output
true
source
GeometryOps.equals Method
julia
equals(
    ::GI.LinearRingTrait, l1,
    ::GI.LinearRingTrait, l2,
)::Bool

Two linear rings are equal if they share the same set of points going along the curve. Note that rings are closed by definition, so they can have, but don't need, a repeated last point to be equal.

source
GeometryOps.equals Method
julia
equals(
    ::GI.LinearRingTrait, l1,
    ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
)::Bool

A linear ring and a line/linestring are equal if they share the same set of points going along the curve. Note that lines aren't closed by definition, but rings are, so the line must have a repeated last point to be equal

source
GeometryOps.equals Method
julia
equals(::GI.MultiPointTrait, mp1, ::GI.MultiPointTrait, mp2)::Bool

Two multipoints are equal if they share the same set of points.

source
GeometryOps.equals Method
julia
equals(::GI.MultiPointTrait, mp1, ::GI.PointTrait, p2)::Bool

A point and a multipoint are equal if the multipoint is composed of a single point that is equivalent to the given point.

source
GeometryOps.equals Method
julia
equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool

Two multipolygons are equal if they share the same set of polygons.

source
GeometryOps.equals Method
julia
equals(::GI.MultiPolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool

A polygon and a multipolygon are equal if the multipolygon is composed of a single polygon that is equivalent to the given polygon.

source
GeometryOps.equals Method
julia
equals(::GI.PointTrait, p1, ::GI.MultiPointTrait, mp2)::Bool

A point and a multipoint are equal if the multipoint is composed of a single point that is equivalent to the given point.

source
GeometryOps.equals Method
julia
equals(::GI.PointTrait, p1, ::GI.PointTrait, p2)::Bool

Two points are the same if they have the same x and y (and z if 3D) coordinates.

source
GeometryOps.equals Method
julia
equals(::GI.PolygonTrait, geom_a, ::GI.MultiPolygonTrait, geom_b)::Bool

A polygon and a multipolygon are equal if the multipolygon is composed of a single polygon that is equivalent to the given polygon.

source
GeometryOps.equals Method
julia
equals(::GI.PolygonTrait, geom_a, ::GI.PolygonTrait, geom_b)::Bool

Two polygons are equal if they share the same exterior edge and holes.

source
GeometryOps.equals Method

This functionality is experimental and may change at any time.

source
GeometryOps.equals Method
julia
equals(
    ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
    ::GI.LinearRingTrait, l2,
)::Bool

A line/linestring and a linear ring are equal if they share the same set of points going along the curve. Note that lines aren't closed by definition, but rings are, so the line must have a repeated last point to be equal

source
GeometryOps.equals Method
julia
equals(
    ::Union{GI.LineTrait, GI.LineStringTrait}, l1,
    ::Union{GI.LineTrait, GI.LineStringTrait}, l2,
)::Bool

Two lines/linestrings are equal if they share the same set of points going along the curve. Note that lines/linestrings aren't closed by definition.

source
GeometryOps.equals Method
julia
equals(::T, geom_a, ::T, geom_b)::Bool

Two geometries of the same type, which don't have a equals function to dispatch off of should throw an error.

source
GeometryOps.extent_to_polygon Method
julia
extent_to_polygon(ext::Extents.Extent)

Convert an extent to a polygon.

Examples

julia
import GeometryOps as GO, Extents
    
ext = Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0))
GO.extent_to_polygon(ext)
# output
GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(1.0, 1.0), … (3) … , (1.0, 1.0)])])
source
GeometryOps.extract_segment_strings Method
julia
extract_segment_strings(rg::RelateGeometry, is_a::Bool, ext_filter)

Extract RelateSegmentStrings from the geometry which intersect a given extent (one per line, one per polygon ring). If ext_filter is nothing all edges are extracted.

Warning

nothing here means no filter (Java's prepared-mode null), while get_extent(rg) returns nothing for an empty geometry, where JTS's null Envelope intersects nothing. Never forward an empty geometry's extent as the filter — callers (the engine's computeAtEdges port) must early-return on empty inputs before extraction.

source
GeometryOps.finish! Method
julia
finish!(n::RelateNode, is_area_interior_a::Bool, is_area_interior_b::Bool)

Computes the final topology for the edges around this node. Although nodes lie on the boundary of areas or the interior of lines, in a mixed GC they may also lie in the interior of an area. This changes the locations of the sides and line to Interior.

Port of RelateNode.finish.

source
GeometryOps.flip Method
julia
flip(obj)

Swap all of the x and y coordinates in obj, otherwise keeping the original structure (but not necessarily the original type).

Keywords

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

  • calc_extent: true or false. Whether to calculate the extent. Defaults to false.

source
GeometryOps.forcexy Method
julia
forcexy(geom)

Force the geometry to be 2D. Works on any geometry, vector of geometries, feature collection, or table!

source
GeometryOps.forcexyz Function
julia
forcexyz(geom, z = 0)

Force the geometry to be 3D. Works on any geometry, vector of geometries, feature collection, or table!

The z parameter is the default z value - if a point has no z value, it will be set to this value. If it does, then the z value will be kept.

source
GeometryOps.foreach_pair_of_maybe_intersecting_edges_in_order Method
julia
foreach_pair_of_maybe_intersecting_edges_in_order(
    manifold::M, accelerator::A,
    f_on_each_a::FA,
    f_after_each_a::FAAfter,
    f_on_each_maybe_intersect::FI,
    geom_a,
    geom_b,
    ::Type{T} = Float64
) where {FA, FAAfter, FI, T, M <: Manifold, A <: IntersectionAccelerator}

Decompose geom_a and geom_b into edge lists (unsorted), and then, logically, perform the following iteration:

julia
for (a_edge, i) in enumerate(eachedge(geom_a))
    f_on_each_a(a_edge, i)
    for (b_edge, j) in enumerate(eachedge(geom_b))
        if may_intersect(a_edge, b_edge)
            f_on_each_maybe_intersect(a_edge, b_edge)
        end
    end
    f_after_each_a(a_edge, i)
end

This may not be the exact acceleration that is performed - but it is the logical sequence of events. It also uses the accelerator, and can automatically choose the best one based on an internal heuristic if you pass in an AutoAccelerator.

For example, the SingleSTRtree accelerator is used along with extent thinning to avoid unnecessary edge intersection checks in the inner loop.

source
GeometryOps.get_dimension_real Method
julia
get_dimension_real(rg::RelateGeometry)

Gets the actual non-empty dimension of the geometry. Zero-length LineStrings are treated as Points.

source
GeometryOps.get_location Method
julia
get_location(label, index, position, is_forward) -> Int8

The location of position (POS_LEFT / POS_RIGHT / POS_ON) of input index, for a containing half-edge whose orientation is is_forward. When the edge is the reverse half-edge (is_forward == false) the Left/Right stored sides are swapped, so a single label serves both members of a symmetric pair (port of JTS OverlayLabel.getLocation(index, position, isForward)).

source
GeometryOps.get_polygonal Method
julia
get_polygonal(nss::NodeSections, is_a::Bool)

The parent polygonal geometry of the first section of input geometry is_a that has one, or nothing.

Port of NodeSections.getPolygonal(boolean isA).

source
GeometryOps.get_polygonal Method
julia
get_polygonal(ns::NodeSection)

Gets the polygon this section is part of. Will be nothing if section is not on a polygon boundary.

Port of NodeSection.getPolygonal.

source
GeometryOps.intersection Method

This functionality is experimental and may change at any time.

source
GeometryOps.intersection Method
julia
intersection(geom_a, geom_b, [T::Type]; target::Type, fix_multipoly = UnionIntersectingPolygons())

Return the intersection between two geometries as a list of geometries. Return an empty list if none are found. The type of the list will be constrained as much as possible given the input geometries. Furthermore, the user can provide a target type as a keyword argument and a list of target geometries found in the intersection will be returned. The user can also provide a float type that they would like the points of returned geometries to be. If the user is taking a intersection involving one or more multipolygons, and the multipolygon might be comprised of polygons that intersect, if fix_multipoly is set to an IntersectingPolygons correction (the default is UnionIntersectingPolygons()), then the needed multipolygons will be fixed to be valid before performing the intersection to ensure a correct answer. Only set fix_multipoly to nothing if you know that the multipolygons are valid, as it will avoid unneeded computation.

Example

julia
import GeoInterface as GI, GeometryOps as GO

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
inter_points = GO.intersection(line1, line2; target = GI.PointTrait())
GI.coordinates.(inter_points)

# output
1-element Vector{Vector{Float64}}:
 [125.58375366067548, -14.83572303404496]
source
GeometryOps.intersection_area Function
julia
intersection_area(alg::Algorithm, geom_a, geom_b, [T = Float64]; kwargs...)

Area of the intersection of geom_a and geom_b, computed without constructing the intersection geometry.

Equivalent to area(manifold(alg), intersection(alg, geom_a, geom_b)), but the result polygon is never built. On Spherical() the result is in square units of the manifold radius.

The algorithm is required, and the supported list is closed:

algorithmmanifoldsnotes
ConvexConvexSutherlandHodgmanPlanar(), Spherical()takes the same cache keyword as intersection; with one, the call allocates nothing
OverlayNGPlanar(), Spherical()T sets the accumulator and return type only — the arrangement always runs at the algorithm's point_type
FosterHormannClippingas intersectionaccumulates during the trace for hole-free polygon pairs; other inputs delegate to the polygon path

Example

julia
import GeometryOps as GO, GeoInterface as GI

a = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
b = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]])

GO.intersection_area(GO.OverlayNG(), a, b)  # 1.0

alg = GO.ConvexConvexSutherlandHodgman()
cache = GO.SutherlandHodgmanCache(alg)
GO.intersection_area(alg, a, b; cache)      # 1.0, allocation-free
source
GeometryOps.intersection_points Method
julia
intersection_points(geom_a, geom_b, [T::Type])

Return a list of intersection tuple points between two geometries. If no intersection points exist, returns an empty list.

Example

jldoctest

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)]) line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)]) inter_points = GO.intersection_points(line1, line2)

**output**

1-element Vector{Tuple{Float64, Float64}}:  (125.58375366067548, -14.83572303404496)


<Badge type="info" class="source-link" text="source"><a href="https://github.com/JuliaGeo/GeometryOps.jl/blob/a2a8b8283322586a2cb4913a213b973e9108e257/src/methods/clipping/intersection.jl#L202-L220" target="_blank" rel="noreferrer">source</a></Badge>

</details>

<details class='jldocstring custom-block' open>
<summary><a id='GeometryOps.intersects-Tuple{Any, Any}' href='#GeometryOps.intersects-Tuple{Any, Any}'><span class="jlbinding">GeometryOps.intersects</span></a> <Badge type="info" class="jlObjectType jlMethod" text="Method" /></summary>



```julia
intersects(geom1, geom2)::Bool

Return true if the interiors or boundaries of the two geometries interact.

intersects returns the exact opposite result of disjoint.

Example

julia
import GeoInterface as GI, GeometryOps as GO

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
GO.intersects(line1, line2)

# output
true
source
GeometryOps.intersects Method
julia
intersects(g1)

Return a function that checks if its input intersects g1. This is equivalent to x -> intersects(x, g1).

source
GeometryOps.intersects Method

This functionality is experimental and may change at any time.

source
GeometryOps.is_containing_segment Method
julia
is_containing_segment(ss::RelateSegmentString, seg_index::Integer, pt)

Tests if a segment intersection point has that segment as its canonical containing segment. Segments are half-closed, and contain their start point but not the endpoint, except for the final segment in a non-closed segment string, which contains its endpoint as well. This test ensures that vertices are assigned to a unique segment in a segment string. In particular, this avoids double-counting intersections which lie exactly at segment endpoints.

source
GeometryOps.is_polygonal Method
julia
is_polygonal(rg::RelateGeometry)

Tests whether the geometry has polygonal topology. This is not the case if it is a GeometryCollection containing more than one polygon (since they may overlap or be adjacent). The significance is that polygonal topology allows more assumptions about the location of boundary vertices.

source
GeometryOps.is_self_noding_required Method
julia
is_self_noding_required(rg::RelateGeometry)

Indicates whether the geometry requires self-noding for correct evaluation of specific spatial predicates. Self-noding is required for geometries which may self-cross — i.e. lines, and overlapping elements in GeometryCollections. Self-noding is not required for polygonal geometries, since they can only touch at vertices.

source
GeometryOps.is_self_noding_required Method
julia
is_self_noding_required(tc::TopologyComputer)

Indicates whether the input geometries require self-noding for correct evaluation of specific spatial predicates. Self-noding is required for geometries which may have self-crossing linework, or may have lines lying in the boundary of an area. This ensures that node locations match in situations where a self-crossing and mutual crossing occur at the same logical location (here via the D3 coincidence-merge pass, since node identities are symbolic).

Currently self-noding is required for:

  • A geoms which require self-noding (lines or GCs, except for single-polygon GCs)

  • B geoms which are mixed A/L GCs

Port of TopologyComputer.isSelfNodingRequired.

source
GeometryOps.isclockwise Method
julia
isclockwise(line::Union{LineString, Vector{Position}})::Bool

Take a ring and return true if the line goes clockwise, or false if the line goes counter-clockwise. "Going clockwise" means, mathematically,

Example

julia
julia> import GeoInterface as GI, GeometryOps as GO
julia> ring = GI.LinearRing([(0, 0), (1, 1), (1, 0), (0, 0)]);
julia> GO.isclockwise(ring)
# output
true
source
GeometryOps.isconcave Method
julia
isconcave(poly::Polygon)::Bool

Take a polygon and return true or false as to whether it is concave or not.

Examples

julia
import GeoInterface as GI, GeometryOps as GO

poly = GI.Polygon([[(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]])
GO.isconcave(poly)

# output
false
source
GeometryOps.lazy_edge_extents Method
julia
lazy_edge_extents(geom)

Return an iterator over the extents of the edges (line segments) of geom. This is lazy but nonallocating.

source
GeometryOps.lazy_edgelist Method
julia
lazy_edgelist(geom, [::Type{T}])
lazy_edgelist(m::Manifold, geom, [::Type{T}])

Return an iterator over GI.Line objects with attached extents.

source
GeometryOps.locate Method
julia
locate(ael::AdjacentEdgeLocator, p)

Location (LOC_INTERIOR or LOC_BOUNDARY) of point p, which must lie on at least one polygon edge of the locator's geometry, under union semantics.

source
GeometryOps.locate Method
julia
locate(loc::IndexedPointInAreaLocator, p)

The location (LOC_* code) of point p in the locator's areal geometry. Port of IndexedPointInAreaLocator.locate.

source
GeometryOps.locate Method
julia
locate(loc::RelatePointLocator, p)

The location (LOC_* code) of point p relative to the locator's geometry, under GC union semantics.

source
GeometryOps.locate_area_vertex Method
julia
locate_area_vertex(rg::RelateGeometry, pt)

Locates a vertex of a polygon. A vertex of a Polygon or MultiPolygon is on the boundary; but a vertex of an overlapped polygon in a GeometryCollection may be in the interior.

source
GeometryOps.locate_line_end_with_dim Method
julia
locate_line_end_with_dim(loc::RelatePointLocator, p)

Locates a line endpoint, as a DL_* dimension-location code. In a mixed-dim GC, the line end point may also lie in an area. In this case the area location is reported. Otherwise, the dimloc is either DL_LINE_BOUNDARY or DL_LINE_INTERIOR, depending on the endpoint valence and the BoundaryNodeRule in place.

source
GeometryOps.locate_node Method
julia
locate_node(loc::RelatePointLocator, p, parent_polygonal)

The location (LOC_* code) of a point p which is known to be a node of the geometry (i.e. a vertex or on an edge). parent_polygonal is the polygonal element the point is a node of (or nothing).

source
GeometryOps.locate_node_with_dim Method
julia
locate_node_with_dim(loc::RelatePointLocator, p, parent_polygonal)

The dimension-location (DL_* code) of a point p which is known to be a node of the geometry.

source
GeometryOps.locate_with_dim Method
julia
locate_with_dim(loc::RelatePointLocator, p)

Computes the topological location (DL_* dimension-location code) of a single point in a geometry, including the dimension of the geometry element the point is located in (if not in the exterior). It handles both single-element and multi-element geometries. The algorithm for multi-part geometries takes into account the SFS Boundary Determination Rule.

source
GeometryOps.matches_entry Method

Match a single matrix entry against a pattern code (JTS IntersectionMatrix.matches).

source
GeometryOps.minimum_bounding_circle Function
julia
minimum_bounding_circle([algorithm], geometry)

Compute the minimum bounding circle of geometry.

Returns a circle geometry containing all points of the input. For planar geometries, returns a PlanarCircle; for spherical geometries, returns a SphericalCap.

Return type subject to change

The concrete return type (currently PlanarCircle for planar manifold) may change in future versions without a breaking release. However, the return type will always implement GeoInterface, so code using GeoInterface methods (e.g., GI.getexterior, GI.getpoint) will remain compatible.

Arguments

  • algorithm: The algorithm to use. Defaults to Welzl() which uses Welzl's expected O(n) algorithm.

  • geometry: Any geometry compatible with GeoInterface, or a vector of point-like objects.

Example

julia
import GeometryOps as GO, GeoInterface as GI

# From points
points = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)]
circle = GO.minimum_bounding_circle(points)

# From any geometry
polygon = GI.Polygon([[(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]])
circle = GO.minimum_bounding_circle(polygon)

# Access via GeoInterface for forward compatibility
ring = GI.getexterior(circle)
source
GeometryOps.node_point Method
julia
node_point(arr::NodedArrangement{P,T}, id) -> T

The realized output coordinate of node id, in the arrangement's output point type T (see output_point_type), memoized in the node table (design §2.6). The only place a constructed coordinate enters the substrate.

source
GeometryOps.overlaps Method
julia
overlaps(geom1, geom2)::Bool

Compare two Geometries of the same dimension and return true if their intersection set results in a geometry different from both but of the same dimension. This means one geometry cannot be within or contain the other and they cannot be equal

Examples

julia
import GeometryOps as GO, GeoInterface as GI
poly1 = GI.Polygon([[(0,0), (0,5), (5,5), (5,0), (0,0)]])
poly2 = GI.Polygon([[(1,1), (1,6), (6,6), (6,1), (1,1)]])

GO.overlaps(poly1, poly2)
# output
true
source
GeometryOps.overlaps Method
julia
overlaps(g1)

Return a function that checks if its input overlaps g1. This is equivalent to x -> overlaps(x, g1).

source
GeometryOps.overlaps Method

This functionality is experimental and may change at any time.

source
GeometryOps.polygon_node_convert Method
julia
polygon_node_convert(m::Manifold, poly_sections::Vector{<:NodeSection}; exact)

Converts the node sections at a polygon node where a shell and one or more holes touch, or two or more holes touch. This converts the node topological structure from the OGC "touching-rings" (AKA "minimal-ring") model to the equivalent "self-touch" (AKA "inverted/exverted ring" or "maximal ring") model. In the "self-touch" model the converted NodeSection corners enclose areas which all lie inside the polygon (i.e. they do not enclose hole edges). This allows RelateNode (Task 17) to use simple area-additive semantics for adding edges and propagating edge locations.

The input node sections are assumed to have canonical orientation (CW shells and CCW holes). The arrangement of shells and holes must be topologically valid. Specifically, the node sections must not cross or be collinear.

This supports multiple shell-shell touches (including ones containing holes), and hole-hole touches. This generalizes the relate algorithm to support both the OGC model and the self-touch model.

Converts a list of sections of valid polygon rings to have "self-touching" structure. There are the same number of output sections as input ones. Sorts (and thereby mutates) poly_sections; returns the converted sections.

Port of PolygonNodeConverter.convert. The angle sort goes through edge_angle_compare (the NodeSection.EdgeAngleComparator port), which is why the manifold and the exact flag are threaded in (the Java method is geometry-context-free).

source
GeometryOps.polygon_to_line Method
julia
polygon_to_line(poly::Polygon)

Converts a Polygon to LineString or MultiLineString

Examples

julia
import GeometryOps as GO, GeoInterface as GI

poly = GI.Polygon([[(-2.275543, 53.464547), (-2.275543, 53.489271), (-2.215118, 53.489271), (-2.215118, 53.464547), (-2.275543, 53.464547)]])
GO.polygon_to_line(poly)
# output
GeoInterface.Wrappers.LineString{false, false}([(-2.275543, 53.464547), … (3) … , (-2.275543, 53.464547)])
source
GeometryOps.polygonize Method
julia
polygonize(A::AbstractMatrix{Bool}; kw...)
polygonize(f, A::AbstractMatrix; kw...)
polygonize(xs, ys, A::AbstractMatrix{Bool}; kw...)
polygonize(f, xs, ys, A::AbstractMatrix; kw...)

Polygonize an AbstractMatrix of values, currently to a single class of polygons.

Returns a MultiPolygon for Bool values and f return values, and a FeatureCollection of Features holding MultiPolygon for all other values.

Function f should return either true or false or a transformation of values into simpler groups, especially useful for floating point arrays.

If xs and ys are ranges, they are used as the pixel/cell center points. If they are Vector of Tuple they are used as the lower and upper bounds of each pixel/cell.

Keywords

  • minpoints: ignore polygons with less than minpoints points.

  • values: the values to turn into polygons. By default these are union(A), If function f is passed these refer to the return values of f, by default union(map(f, A). If values Bool, false is ignored and a single MultiPolygon is returned rather than a FeatureCollection.

Example

julia
using GeometryOps
A = rand(100, 100)
multipolygon = polygonize(>(0.5), A);
source
GeometryOps.prepare Method
julia
prepare(alg::RelateNG, a; validate = <manifold-dependent>)::PreparedRelate

This functionality is experimental and may change at any time.

prepare is the generic entry point for prepared-geometry optimizations in GeometryOps; RelateNG is currently the only algorithm implementing it.

Creates a prepared relate instance to optimize the repeated evaluation of relationships against the single geometry a.

Port of RelateNG.prepare(Geometry) (the algorithm's boundary_rule plays the role of the prepare(Geometry, BoundaryNodeRule) overload). The A-side RelateGeometry is constructed with is_prepared = true, and the lazy caches that the Java instance accumulates across evaluations are forced eagerly:

  • the RelatePointLocator — the prepared flag selects the per-polygonal-element IndexedPointInAreaLocator caches, as in Java (see point_locator.jl / indexed_point_in_area.jl);

  • the unique-point set, when a has effective dimension P (the only case the P/P fast path consults it);

  • the A segment strings, extracted ONCE without an interaction-envelope filter (Java's prepared-mode envExtract = null in computeEdgesMutual) so the cache serves any future B, plus the prebuilt segment tree over them.

Note

Predicates whose evaluation requires self-noding (is_self_noding_required) bypass the cached edges entirely: as in the Java prepared branch, computeEdgesAll re-extracts the A edges per evaluation, filtered by the A/B interaction envelope.

Validation

validate controls a ring self-crossing check over the prepared geometry: a self-join of the segment set (reusing the prepared edge index) that detects PROPER crossings — transversal, interior to both edges — between non-adjacent edges of the same polygonal element. Shared-endpoint adjacency is excluded and vertex touches are NOT flagged: the scope is exactly the crossing class that breaks the engine's containment parity, not full OGC validity. On the first crossing found an ArgumentError is thrown naming the ring and the edge pair ("edge i crosses edge j"); the documented remedy is the CrossingEdgeSplit correction.

The default is manifold-dependent, and deliberately so:

  • Spherical → validate = true. A planar-valid ring whose edges cross when reinterpreted as great-circle arcs is undetectable by standard planar tooling (planar validity checks pass it), and undetected it can invert containment globally — the figure-eight's lobes cancel the curvature the interior bootstrap reads, so every query lands on the wrong side (Natural Earth 110m Sudan is a real instance). The check is a small fraction of the ~100 ms spherical prepare build.

  • Planar → validate = false. Planar invalidity of the same class is visible to ordinary planar tools, and the planar engine degrades gracefully under even-odd ray-crossing parity instead of inverting; JTS/GEOS prepared geometries do not validate either. A planar prepare costs ~600 µs, which a validation join would dominate, destroying the build-cost amortization. Pass validate = true to opt in.

source
GeometryOps.prepare Method
julia
prepare(alg::RelateNG, p::PreparedRelate; validate = <manifold-dependent>)::PreparedRelate

This functionality is experimental and may change at any time.

Returns p under the algorithm it was prepared with, so a PreparedRelate can be passed as the A geometry to any RelateNG entry point. Under any other alg — every setting is baked into the prepared structures — it is prepared afresh from p.input, at the full cost of prepare on every call, so prepare under the algorithm you query with. validate applies only to that rebuild.

source
GeometryOps.prepare_sections! Method
julia
prepare_sections!(nss::NodeSections)

Sorts the sections (by compare_to, the Java natural ordering) so that:

  • lines are before areas

  • edges from the same polygon are contiguous

Port of NodeSections.prepareSections.

source
GeometryOps.process_edge_intersections! Method
julia
process_edge_intersections!(computer, ssa_list, ssb_list,
    accelerator = AutoAccelerator(); m, exact)

Enumerate all extent-interacting segment pairs between the A-side segment strings ssa_list and the B-side segment strings ssb_list, feeding each pair through process_intersections! so the intersections are recorded on computer.

accelerator selects the enumeration strategy:

  • NestedLoop: a plain double loop over string pairs and segment pairs, with a per-pair segment-extent disjointness skip (on Planar).

  • Any tree-backed accelerator (canonically DoubleNaturalTree): a spatial index (_relate_edge_index, a natural-order RTree carrying segment owners as data) is built over the per-segment extents of each side and traversed with SpatialTreeInterface.dual_depth_first_search under the Extents.intersects predicate.

  • AutoAccelerator: picks NestedLoop below the clipping size threshold (GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS) and on manifolds without a segment-extent kernel (neither Planar nor Spherical), and the tree path otherwise.

After each processed pair is_result_known(computer) is consulted and the enumeration stops early once the predicate value is determined (the port of the Java noder's isDone() hook used by EdgeSetIntersector.process).

Warning

This enumerates A×B pairs only. JTS's EdgeSetIntersector also feeds A×A and B×B pairs (self-noding) with an id-ordering guard so each unordered pair is processed once. Calling this with the same list on both sides would process every pair twice — the engine's computeAtEdges port uses process_self_intersections! for the self-pair path instead.

source
GeometryOps.process_intersections! Method
julia
process_intersections!(tc::TopologyComputer, ss0, seg_index0, ss1, seg_index1; m, exact)

Classify the intersection of segment seg_index0 of ss0 with segment seg_index1 of ss1 and record any intersections in tc. The strings are ordered so the A geometry's string is processed first (the computer's A/B matrix updates rely on it).

Port of EdgeSegmentIntersector.processIntersections.

source
GeometryOps.process_self_intersections! Method
julia
process_self_intersections!(computer, ss_list,
    accelerator = AutoAccelerator(); m, exact)

Enumerate all extent-interacting segment pairs within the segment-string list ss_list (each unordered pair once, never a segment against itself), feeding each through process_intersections! — the self-noding (A×A or B×B) counterpart of process_edge_intersections!.

After each processed pair is_result_known(computer) is consulted and the enumeration stops early once the predicate value is determined.

source
GeometryOps.radius Method
julia
radius(circle::PlanarCircle)

Return the radius of the circle. Computes sqrt(radius_squared).

source
GeometryOps.relate Method
julia
relate(p::PreparedRelate, b, im_pattern::AbstractString)::Bool

This functionality is experimental and may change at any time.

Tests whether the topological relationship of the prepared geometry to b matches the DE-9IM pattern. Port of RelateNG.evaluate(Geometry, String).

source
GeometryOps.relate Method
julia
relate(p::PreparedRelate, b)::DE9IM

This functionality is experimental and may change at any time.

Computes the DE-9IM matrix for the topological relationship of the prepared geometry to b. Port of the instance method RelateNG.evaluate(Geometry).

source
GeometryOps.relate Method
julia
relate([alg::RelateNG], a, b, im_pattern::AbstractString)::Bool

This functionality is experimental and may change at any time.

Tests whether the topological relationship between geometries a and b matches the DE-9IM matrix pattern im_pattern (9 characters over 012TF*).

Port of RelateNG.relate(Geometry, Geometry, String).

source
GeometryOps.relate Method
julia
relate([alg::RelateNG], a, b)::DE9IM

This functionality is experimental and may change at any time.

Computes the DE9IM matrix for the topological relationship between geometries a and b.

Port of RelateNG.relate(Geometry, Geometry).

source
GeometryOps.relate_predicate Method
julia
relate_predicate(p::PreparedRelate, predicate::TopologyPredicate, b)::Bool

This functionality is experimental and may change at any time.

Tests whether the topological relationship of the prepared geometry to b satisfies the predicate. Port of the instance method RelateNG.evaluate(Geometry, TopologyPredicate) in prepared mode.

source
GeometryOps.relate_predicate Method
julia
relate_predicate(alg::RelateNG, predicate::TopologyPredicate, a, b)::Bool

This functionality is experimental and may change at any time.

Tests whether the topological relationship between geometries a and b satisfies the given TopologyPredicate. This is the core evaluation entry point (the port of RelateNG.evaluate(Geometry, TopologyPredicate), via the static RelateNG.relate(a, b, pred)).

Note

Predicates are mutable accumulators — pass a freshly constructed one (e.g. pred_intersects()) per evaluation.

source
GeometryOps.rk_crossing_dirs_ccw Method

CCW cyclic order of the four half-edge directions incident to the proper crossing of (a0,a1) × (b0,b1), starting from a1. Since the crossing is proper, b0/b1 are strictly on opposite sides of line(a0,a1): if b1 is to the left, CCW order is (a1, b1, a0, b0), else (a1, b0, a0, b1).

source
GeometryOps.segmentize Method
julia
segmentize([method = Planar()], geom; max_distance::Real, threaded)

Segmentize a geometry by adding extra vertices to the geometry so that no segment is longer than a given distance. This is useful for plotting geometries with a limited number of vertices, or for ensuring that a geometry is not too "coarse" for a given application.

Arguments

  • method::Manifold = Planar(): The method to use for segmentizing the geometry. At the moment, Planar (assumes a flat plane), Spherical (assumes geometry on a sphere and interpolates along great circles) and Geodesic (assumes geometry on the ellipsoidal Earth and uses Vincenty's formulae) are available.

  • geom: The geometry to segmentize. Must be a LineString, LinearRing, Polygon, MultiPolygon, or GeometryCollection, or some vector or table of those.

  • max_distance::Real: The maximum distance between vertices in the geometry. For Planar manifolds, this is in the units of the geometry. For Spherical and Geodesic, this is in the units of the manifold (the units of the radius of the sphere/ellipsoid). By default this is meters.

Spherical and Geodesic both assume that the input geometry is in lon/lat coordinates, in degrees.

Returns a geometry of similar type to the input geometry, but resampled.

source
GeometryOps.signed_area Method
julia
signed_area(geom, [T = Float64])::T

Returns the signed area of a single geometry, based on winding order. This is computed slightly differently for different geometries:

julia
- The signed area of a point is always zero.
- The signed area of a curve is always zero.
- The signed area of a polygon is computed with the shoelace formula and is
positive if the polygon coordinates wind clockwise and negative if
counterclockwise.
- You cannot compute the signed area of a multipolygon as it doesn't have a
meaning as each sub-polygon could have a different winding order.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.signed_distance Method
julia
signed_distance(point, geom, ::Type{T} = Float64)::T

Calculates the signed distance from the geometry geom to the given point. Points within geom have a negative signed distance, and points outside of geom have a positive signed distance. - The signed distance from a point to a point, line, linestring, or linear ring is equal to the distance between the two. - The signed distance from a point to a polygon is negative if the point is within the polygon and is positive otherwise. The value of the distance is the minimum distance from the point to an edge of the polygon. This includes edges created by holes. - The signed distance from a point to a multigeometry or a geometry collection is the minimum signed distance between the point and any of the sub-geometries.

Result will be of type T, where T is an optional argument with a default value of Float64.

source
GeometryOps.simplify Method
julia
simplify(obj; kw...)
simplify(::SimplifyAlg, obj; kw...)

Simplify a geometry, feature, feature collection, or nested vectors or a table of these.

RadialDistance, DouglasPeucker, or VisvalingamWhyatt algorithms are available, listed in order of increasing quality but decreasing performance.

PoinTrait and MultiPointTrait are returned unchanged.

The default behaviour is simplify(DouglasPeucker(; kw...), obj). Pass in other SimplifyAlg to use other algorithms.

Keywords

  • prefilter_alg: SimplifyAlg algorithm used to pre-filter object before using primary filtering algorithm.

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

  • calc_extent: true or false. Whether to calculate the extent. Defaults to false.

Keywords for DouglasPeucker are allowed when no algorithm is specified:

Keywords

  • ratio: the fraction of points that should remain after simplify. Useful as it will generalise for large collections of objects.

  • number: the number of points that should remain after simplify. Less useful for large collections of mixed size objects.

  • tol: the minimum distance a point will be from the line joining its neighboring points.

Example

Simplify a polygon to have six points:

julia
import GeoInterface as GI
import GeometryOps as GO

poly = GI.Polygon([[
    [-70.603637, -33.399918],
    [-70.614624, -33.395332],
    [-70.639343, -33.392466],
    [-70.659942, -33.394759],
    [-70.683975, -33.404504],
    [-70.697021, -33.419406],
    [-70.701141, -33.434306],
    [-70.700454, -33.446339],
    [-70.694274, -33.458369],
    [-70.682601, -33.465816],
    [-70.668869, -33.472117],
    [-70.646209, -33.473835],
    [-70.624923, -33.472117],
    [-70.609817, -33.468107],
    [-70.595397, -33.458369],
    [-70.587158, -33.442901],
    [-70.587158, -33.426283],
    [-70.590591, -33.414248],
    [-70.594711, -33.406224],
    [-70.603637, -33.399918]]])

simple = GO.simplify(poly; number=6)
GI.npoint(simple)

# output
6
source
GeometryOps.smooth Method
julia
smooth(alg::Algorithm, geom)
smooth(geom; kw...)

Smooths a geometry using the provided algorithm.

The default algorithm is Chaikin(), which can be used on the spherical or planar manifolds.

source
GeometryOps.symdifference Method
julia
symdifference([alg::OverlayNG], geom_a, geom_b)
symdifference(manifold::Manifold, geom_a, geom_b)

This functionality is experimental and may change at any time.

The symmetric difference of two geometries: everything that lies in exactly one of them, i.e. union(difference(a, b), difference(b, a)).

julia
import GeoInterface as GI, GeometryOps as GO

a = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
b = GI.Polygon([[(1.0, 0.0), (3.0, 0.0), (3.0, 2.0), (1.0, 2.0), (1.0, 0.0)]])
GO.area(GO.symdifference(a, b))

# output
4.0

Engine

symdifference is OverlayNG-only: there is no Foster–Hormann symmetric difference, and none is planned — symmetric difference is the operation the exact arrangement gets for free and the ent/exit tracer does not.

That makes it the one member of the overlay family whose algorithm-free form does not run Foster–Hormann: symdifference(a, b) is symdifference(OverlayNG(Planar()), a, b), and symdifference(m, a, b) is symdifference(OverlayNG(m), a, b). The defaults of intersection, union and difference are unaffected.

See OverlayNG for the result shape, the input contract, and the spherical full-sphere limitation.

source
GeometryOps.t_value Method
julia
t_value(sᵢ, sᵢ₊₁, rᵢ, rᵢ₊₁)

Returns the "T-value" as described in Hormann's presentation [2] on how to calculate the mean-value coordinate.

Here, sᵢ is the vector from vertex vᵢ to the point, and rᵢ is the norm (length) of sᵢ. s must be Point and r must be real numbers.

ᵢᵢᵢ₊₁ᵢᵢ₊₁ᵢᵢ₊₁


<Badge type="info" class="source-link" text="source"><a href="https://github.com/JuliaGeo/GeometryOps.jl/blob/a2a8b8283322586a2cb4913a213b973e9108e257/src/methods/barycentric.jl#L211-L227" target="_blank" rel="noreferrer">source</a></Badge>

</details>

<details class='jldocstring custom-block' open>
<summary><a id='GeometryOps.to_edgelist-Union{Tuple{Any}, Tuple{T}, Tuple{Any, Type{T}}} where T' href='#GeometryOps.to_edgelist-Union{Tuple{Any}, Tuple{T}, Tuple{Any, Type{T}}} where T'><span class="jlbinding">GeometryOps.to_edgelist</span></a> <Badge type="info" class="jlObjectType jlMethod" text="Method" /></summary>



```julia
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 UnitSphericalPoints — each edge carries the 3D UnitSpherical.spherical_arc_extent of its great-circle arc.

source
GeometryOps.to_edgelist Method
julia
to_edgelist(ext::E, geom, [::Type{T}])::(::Vector{GI.Line}, ::Vector{Int})

Filter the edges of geom for those that intersect ext, and return:

  • a vector of GI.Line objects with attached extents,

  • a vector of indices into the original geometry.

source
GeometryOps.to_edges Method
julia
to_edges()

Convert any geometry or collection of geometries into a flat vector of Tuple{Tuple{Float64,Float64},Tuple{Float64,Float64}} edges.

source
GeometryOps.touches Method
julia
touches([manifold::Manifold], geom1, geom2)::Bool

Return true if the first geometry touches the second geometry. In other words, the two interiors cannot interact, but one of the geometries must have a boundary point that interacts with either the other geometry's interior or boundary.

Examples

julia
import GeometryOps as GO, GeoInterface as GI

l1 = GI.Line([(0.0, 0.0), (1.0, 0.0)])
l2 = GI.Line([(1.0, 1.0), (1.0, -1.0)])

GO.touches(l1, l2)
# output
true
source
GeometryOps.touches Method
julia
touches(g1)

Return a function that checks if its input touches g1. This is equivalent to x -> touches(x, g1).

source
GeometryOps.touches Method

This functionality is experimental and may change at any time.

source
GeometryOps.transform Method
julia
transform(f, obj)

Apply a function f to all the points in obj.

Points will be passed to f as an SVector to allow using CoordinateTransformations.jl and Rotations.jl without hassle.

SVector is also a valid GeoInterface.jl point, so will work in all GeoInterface.jl methods.

Example

julia
julia> import GeoInterface as GI

julia> import GeometryOps as GO

julia> geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]);

julia> f = CoordinateTransformations.Translation(3.5, 1.5)
Translation(3.5, 1.5)

julia> GO.transform(f, geom)
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.Linea
rRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}(StaticArraysCo
re.SVector{2, Float64}[[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Float64}}, Nothing, Nothing}(StaticA
rraysCore.SVector{2, Float64}[[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]], nothing, nothing)], nothing, nothing)

With Rotations.jl you need to actually multiply the Rotation by the SVector point, which is easy using an anonymous function.

julia
julia> using Rotations

julia> GO.transform(p -> one(RotMatrix{2}) * p, geom)
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearR
ing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}(StaticArraysCore.SVe
ctor{2, Int64}[[2, 1], [4, 3], [6, 5], [2, 1]], nothing, nothing), GeoInterface.Wrappers.LinearRing{false, false, Vector{StaticArraysCore.SVector{2, Int64}}, Nothing, Nothing}(StaticArraysCore.SVector{2, Int64
}[[4, 3], [6, 5], [7, 6], [4, 3]], nothing, nothing)], nothing, nothing)
source
GeometryOps.tuples Method
julia
tuples(obj)

Convert all points in obj to Tuples, wherever the are nested.

Returns a similar object or collection of objects using GeoInterface.jl geometries wrapping Tuple points.

Keywords

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

  • calc_extent: true or false. Whether to calculate the extent. Defaults to false.

source
GeometryOps.union Method

This functionality is experimental and may change at any time.

source
GeometryOps.union Method
julia
union(geom_a, geom_b, [::Type{T}]; target::Type, fix_multipoly = UnionIntersectingPolygons())

Return the union between two geometries as a list of geometries. Return an empty list if none are found. The type of the list will be constrained as much as possible given the input geometries. Furthermore, the user can provide a taget type as a keyword argument and a list of target geometries found in the difference will be returned. The user can also provide a float type 'T' that they would like the points of returned geometries to be. If the user is taking a intersection involving one or more multipolygons, and the multipolygon might be comprised of polygons that intersect, if fix_multipoly is set to an IntersectingPolygons correction (the default is UnionIntersectingPolygons()), then the needed multipolygons will be fixed to be valid before performing the intersection to ensure a correct answer. Only set fix_multipoly to false if you know that the multipolygons are valid, as it will avoid unneeded computation.

Calculates the union between two polygons.

Example

julia
import GeoInterface as GI, GeometryOps as GO

p1 = GI.Polygon([[(0.0, 0.0), (5.0, 5.0), (10.0, 0.0), (5.0, -5.0), (0.0, 0.0)]])
p2 = GI.Polygon([[(3.0, 0.0), (8.0, 5.0), (13.0, 0.0), (8.0, -5.0), (3.0, 0.0)]])
union_poly = GO.union(p1, p2; target = GI.PolygonTrait())
GI.coordinates.(union_poly)

# output
1-element Vector{Vector{Vector{Vector{Float64}}}}:
 [[[6.5, 3.5], [5.0, 5.0], [0.0, 0.0], [5.0, -5.0], [6.5, -3.5], [8.0, -5.0], [13.0, 0.0], [8.0, 5.0], [6.5, 3.5]]]
source
GeometryOps.vertex_node Method

Node key of a vertex node: keyed exactly by its coordinate.

source
GeometryOps.voronoi Method
julia
voronoi(geometries, [T = Float64]; clip_polygon = nothing, kwargs...)

Compute the Voronoi tessellation of the points in geometries. Returns a vector of GI.Polygon objects representing the Voronoi cells, in the same order as the input points.

Arguments

  • geometries: Any GeoInterface-compatible geometry or collection of geometries that can be decomposed into points

  • T: Float-type for returned polygons points (default: Float64)

Keyword Arguments

  • clip_polygon: what bounding shape should the Voronoi cells be clipped to? (default: nothing -> clipped to the convex hull) clip_polygon can of several types: (1) a GeoInterface polygon, (2) a two-element tuple where the first element is a list of tuple points and the second element is a list of integer indices to indicate the order of the provided points, or (3) a a two-element tuple where the first element is a tuple of tuple points and the second element is a tuple of integer indices to indicate the order of the provided points

    • crs: The CRS to attach to geometries. Defaults to nothing.

  • rng: random number generator to generating the voronoi tesselation

Warning

This interface only computes the 2-dimensional Voronoi tessellation! Only clipped voronoi tesselations can be created! Only T = Float64 or Float32 are guaranteed good results by the underlying package DelaunayTriangulation.

Note

The polygons are returned in the same order as the input points after flattening. Each polygon corresponds to the Voronoi cell of the point at the same index.

Examples

An example with default clipping to the convex hull.

julia
import GeometryOps as GO
import GeoInterface as GI
using Random

rng = Xoshiro(0)
points = [(rand(rng), rand(rng)) .* 5 for i in range(1, 3)]
GO.voronoi(points; rng = rng)
# output
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(4.310704285977424, 0.42985432929210976), … (2) … , (4.310704285977424, 0.42985432929210976)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7949144210695653, 0.4101636087384888), … (4) … , (3.7949144210695653, 0.4101636087384888)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(2.685897788908803, 0.3678259474564151), … (2) … , (2.685897788908803, 0.3678259474564151)])])

An example with clipping to a GeoInterface polygon.

julia
clip_points = ((0.0,0.0), (5.0,0.0), (5.0,5.0), (0.0,5.0), (0.0,0.0))
clip_order = (1, 2, 3, 4, 1)
clip_poly1 = GI.Polygon([collect(clip_points)])
GO.voronoi(points; clip_polygon = clip_poly1, rng = rng)
# output
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(5.0, 0.0), … (3) … , (5.0, 0.0)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7328227614527916, 0.0), … (3) … , (3.7328227614527916, 0.0)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(0.0, 5.0), … (3) … , (0.0, 5.0)])])

An example with clipping to a tuple of tuples.

julia
clip_poly2 = (clip_points, clip_order) # tuples
GO.voronoi(points; clip_polygon = clip_poly2, rng = rng)
# output
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(5.0, 0.0), … (3) … , (5.0, 0.0)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7328227614527916, 0.0), … (3) … , (3.7328227614527916, 0.0)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(0.0, 5.0), … (3) … , (0.0, 5.0)])])

An example with clipping to a tuple of vectors.

julia
clip_poly3 = (collect(clip_points), collect(clip_order)) # vectors
GO.voronoi(points; clip_polygon = clip_poly3, rng = rng)
# output
3-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(5.0, 0.0), … (3) … , (5.0, 0.0)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(3.7328227614527916, 0.0), … (3) … , (3.7328227614527916, 0.0)])])
 GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(0.0, 5.0), … (3) … , (0.0, 5.0)])])
source
GeometryOps.weighted_mean Method
julia
weighted_mean(weight::Real, x1, x2)

Returns the weighted mean of x1 and x2, where weight is the weight of x1.

Specifically, calculates x1 * weight + x2 * (1 - weight).

Note

The idea for this method is that you can override this for custom types, like Color types, in extension modules.

source
GeometryOps.within Method
julia
within(g1)

Return a function that checks if its input is within g1. This is equivalent to x -> within(x, g1).

source
GeometryOps.within Method

This functionality is experimental and may change at any time.

source

Core types ​

GeometryOpsCore.WGS84_EARTH_INV_FLATTENING Constant

The inverse flattening of the WGS84 ellipsoid

source
GeometryOpsCore.WGS84_EARTH_MEAN_RADIUS Constant

The mean radius of the WGS84 ellipsoid, used for spherical manifold default

source
GeometryOpsCore.WGS84_EARTH_SEMI_MAJOR_RADIUS Constant

The semi-major axis of the WGS84 ellipsoid

source
GeometryOpsCore.Algorithm Type
julia
abstract type Algorithm{M <: Manifold}

The abstract supertype for all GeometryOps algorithms. These define how to perform a particular Operation.

An algorithm may be associated with one or many Manifolds. It may either have the manifold as a field, or have it as a static parameter (e.g. struct GEOS <: Algorithm{Planar}).

Interface

All Algorithms must implement the following methods:

  • rebuild(alg, manifold::Manifold) Rebuild algorithm alg with a new manifold as passed in the second argument. This may error and throw a WrongManifoldException if the manifold is not compatible with that algorithm.

  • manifold(alg::Algorithm) Return the manifold associated with the algorithm.

  • best_manifold(alg::Algorithm, input): Return the best manifold for that algorithm, in the absence of any other context. WARNING: this may change in future and is not stable!

The actual implementation is left to the implementation of that particular Operation.

Notable subtypes

  • AutoAlgorithm: Tells the Operation receiving it to automatically select the best algorithm for its input data.

  • ManifoldIndependentAlgorithm: An abstract supertype for an algorithm that works on any manifold. The manifold must be stored in the algorithm for a ManifoldIndependentAlgorithm, and accessed via manifold(alg).

  • SingleManifoldAlgorithm: An abstract supertype for an algorithm that only works on a single manifold, specified in its type parameter. SingleManifoldAlgorithm{Planar} is a special case that does not have to store its manifold, since that doesn't contain any information. All other SingleManifoldAlgorithms must store their manifold, since they do contain information.

  • NoAlgorithm: A type that indicates no algorithm is to be used, essentially the equivalent of nothing.

source
GeometryOpsCore.Applicator Type
julia
abstract type Applicator{F,T}

An abstract type for applicators that apply a function to a target object.

The type parameter F is the type of the function to apply, and T is the type of the target object.

A common dispatch pattern is to dispatch on F which may also be e.g. a ThreadFunctor.

Interface

All applicators must be callable by an index integer, and define the following methods:

  • rebuild(a::Applicator, f) - swap out the function and return a new applicator.

The calling convention is my_applicator(i::Int), so applicators must define this method.

source
GeometryOpsCore.ApplyToArray Type
julia
ApplyToArray(f, target, arr, kw)

Create an Applicator that applies a function to all elements of arr.

source
GeometryOpsCore.ApplyToFeatures Type
julia
ApplyToFeatures(f, target, fc, kw)

Create an Applicator that applies a function to all features of fc.

source
GeometryOpsCore.ApplyToGeom Type
julia
ApplyToGeom(f, target, geom, kw)

Create an Applicator that applies a function to all sub-geometries of geom.

source
GeometryOpsCore.AutoAlgorithm Type
julia
AutoAlgorithm{T, M <: Manifold}(manifold::M, x::T)

Indicates that the Operation should automatically select the best algorithm for its input data, based on the passed in manifold (may be an AutoManifold) and data x.

The actual implementation is left to the implementation of that particular Operation.

source
GeometryOpsCore.AutoManifold Type
julia
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.

source
GeometryOpsCore.BoolsAsTypes Type
julia
abstract type BoolsAsTypes
source
GeometryOpsCore.False Type
julia
struct False <: BoolsAsTypes

A struct that means false.

source
GeometryOpsCore.Geodesic Type
julia
Geodesic(; semimajor_axis, inv_flattening)

A geodesic manifold means that the geometry is on a 3-dimensional ellipsoid, parameterized by semimajor_axis ( in mathematical parlance) and inv_flattening ().

Usually, this is only relevant for area and segmentization calculations. It becomes more relevant as one grows closer to the poles (or equator).

source
GeometryOpsCore.Manifold Type
julia
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, Spherical, and Geodesic manifolds.

source
GeometryOpsCore.ManifoldIndependentAlgorithm Type
julia
abstract type ManifoldIndependentAlgorithm{M <: Manifold} <: Algorithm{M}

The abstract supertype for a manifold-independent algorithm, i.e., one which may work on any manifold.

The manifold is stored in the algorithm for a ManifoldIndependentAlgorithm, and accessed via manifold(alg).

source
GeometryOpsCore.MissingKeywordInAlgorithmException Type
julia
MissingKeywordInAlgorithmException{Alg, F} <: Exception

An error type which is thrown when a keyword argument is missing from an algorithm.

The alg argument is the algorithm struct, and the keyword argument is the keyword that was missing.

This error message is used in the enforce method.

Usage

This is of course not how you would actually use this error type, but it is how you construct and throw it.

julia
throw(MissingKeywordInAlgorithmException(GEOS(; tokl = 1.0), my_function, :tol))

Real world usage will often look like this:

julia
function my_function(alg::CLibraryPlanarAlgorithm, args...)
    mykwarg = enforce(alg, :mykwarg, my_function) # this will throw an error if :mykwarg is not present in alg
end
source
GeometryOpsCore.NoAlgorithm Type
julia
NoAlgorithm(manifold)

A type that indicates no algorithm is to be used, essentially the equivalent of nothing.

Stores a manifold within itself.

source
GeometryOpsCore.Operation Type
julia
abstract type Operation{Alg <: Algorithm} end

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.

source
GeometryOpsCore.Planar Type
julia
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.

source
GeometryOpsCore.SingleManifoldAlgorithm Type
julia
abstract type SingleManifoldAlgorithm{M <: Manifold} <: Algorithm{M}

The abstract supertype for a single-manifold algorithm, i.e., one which is known to only work on a single manifold.

The manifold may be accessed via manifold(alg).

source
GeometryOpsCore.Spherical Type
julia
Spherical(; radius, oriented = false)

A spherical manifold means that the geometry is on the 3-sphere (but is represented by 2-D longitude and latitude).

oriented selects how the interior of a polygon ring is interpreted:

  • oriented = false (the default): a ring's interior is the region it encloses — the smaller of the two regions it bounds — independent of winding direction. This matches how most of the ecosystem treats unoriented data by default (R's s2/sf, spherely, BigQuery), and means shapefile-convention data (clockwise shells) is read the same way as counterclockwise data. No region larger than a hemisphere can be represented.

  • oriented = true: polygon ring directions are known to be correct — exterior rings counterclockwise, interior rings clockwise — so the interior of the polygon is the region on the left of each ring's stored vertex order (the convention of S2's S2Polygon::InitOriented). This makes regions larger than a hemisphere representable, e.g. "the sphere minus a small cap" as a clockwise ring.

Extended help

Note

The traditional definition of spherical coordinates in physics and mathematics, , 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).

Note

With oriented = true, a ring may denote a region covering most of the sphere. Operations remain correct on such regions, but extent-based pruning degenerates (the region's bounding box is essentially the whole sphere), so spatial predicates against them fall back to slower paths.

Validity is manifold-dependent

A ring that is valid in lon/lat can be invalid on the sphere: two non-adjacent edges may cross when reinterpreted as great-circle arcs (a planar needle a few meters wide is enough — Natural Earth 110m Sudan is a real instance), and no planar validity tool can detect it. Prepared spherical predicates (GeometryOps.prepare) therefore validate against this class by default and throw an "edge i crosses edge j" error; the remedy is the GeometryOps.CrossingEdgeSplit correction, which splits each ring at its crossing points into separate loops (even-odd semantics).

source
GeometryOpsCore.TaskFunctors Type
julia
TaskFunctors(functors, tasks_per_thread)

A struct to hold the functors and tasks_per_thread, for internal use with _maptasks where functions have state that cannot be shared accross threads, such as Proj.Transformation.

functors must be an array or tuple of functors, one per thread, and tasks_per_thread must be an integer. This also allows you to control the number of tasks per thread that _maptasks launches, useful for tuning performance if you like.

source
GeometryOpsCore.TraitTarget Type
julia
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.
source
GeometryOpsCore.True Type
julia
struct True <: BoolsAsTypes

A struct that means true.

source
GeometryOpsCore.WithTrait Type
julia
WithTrait(f)

WithTrait is a functor that applies a function to a trait and an object.

Specifically, the calling convention is for f is changed from f(geom) to f(trait, geom; kw...).

This is useful to keep the trait materialized through the call stack, which can improve inferrability and performance.

source
GeometryOpsCore.WrongManifoldException Type
julia
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 thrown when constructing SingleManifoldAlgorithm types.

source
GeometryOpsCore._InitialValue Type
julia
_InitialValue()

Sentinel value for "no init provided". This is the same way Base mapreduce does this. It's a singleton struct so e.g. similar to Nothing.

This is meant to be the default value for applyreduce's init keyword argument.

source
GeometryOpsCore.apply Method
julia
apply(f, target::Union{TraitTarget, GI.AbstractTrait}, obj; kw...)

Reconstruct a geometry, feature, feature collection, or nested vectors of either using the function f on the target trait.

f(target_geom) => x where x also has the target trait, or a trait that can be substituted. For example, swapping PolgonTrait to MultiPointTrait will fail if the outer object has MultiPolygonTrait, but should work if it has FeatureTrait.

Objects "shallower" than the target trait are always completely rebuilt, like a Vector of FeatureCollectionTrait of FeatureTrait when the target has PolygonTrait and is held in the features. These will always be GeoInterface geometries/feature/feature collections. But "deeper" objects may remain unchanged or be whatever GeoInterface compatible objects f returns.

The result is a functionally similar geometry with values depending on f.

  • threaded: true or false. Whether to use multithreading. Defaults to false.

  • crs: The CRS to attach to geometries. Defaults to nothing.

  • calc_extent: true or false. Whether to calculate the extent. Defaults to false.

Example

Flipped point the order in any feature or geometry, or iterables of either:

julia
import GeoInterface as GI
import GeometryOps as GO
geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]),
                   GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])

flipped_geom = GO.apply(GI.PointTrait, geom) do p
    (GI.y(p), GI.x(p))
end
source
GeometryOpsCore.applyreduce Method
julia
applyreduce(f, op, target::Union{TraitTarget, GI.AbstractTrait}, obj; threaded, init, kw...)

Apply function f to all objects with the target trait, and reduce the result with an op like +.

The order and grouping of application of op is not guaranteed.

If threaded==true threads will be used over arrays and iterables, feature collections and nested geometries.

init specifies the initial value for the reduction. If not provided, the reduction uses the first result as the starting point (like reduce without init). For operations like vcat, you typically don't need to provide init. For numeric reductions like +, you may want to provide init=zero(T) to ensure type stability.

source
GeometryOpsCore.best_manifold Function
julia
best_manifold(alg::Algorithm, input)::Manifold

Return the best Manifold for the algorithm alg based on the given input.

May be any subtype of Manifold.

source
GeometryOpsCore.booltype Function
julia
booltype(x)

Returns a BoolsAsTypes from x, whether it's a boolean or a BoolsAsTypes.

source
GeometryOpsCore.flatten Method
julia
flatten(target::Type{<:GI.AbstractTrait}, obj)
flatten(f, target::Type{<:GI.AbstractTrait}, obj)

Lazily flatten any AbstractArray, iterator, FeatureCollectionTrait, FeatureTrait or AbstractGeometryTrait object obj, so that objects with the target trait are returned by the iterator.

If f is passed in it will be applied to the target geometries.

source
GeometryOpsCore.manifold Function
julia
manifold(alg::Algorithm)::Manifold

Return the manifold associated with the algorithm.

May be any subtype of Manifold.

source
GeometryOpsCore.rebuild Method
julia
rebuild(geom, child_geoms)

Rebuild a geometry from child geometries.

By default geometries will be rebuilt as a GeoInterface.Wrappers geometry, but rebuild can have methods added to it to dispatch on geometries from other packages and specify how to rebuild them.

(Maybe it should go into GeoInterface.jl)

source
GeometryOpsCore.reconstruct Method
julia
reconstruct(geom, components)

Reconstruct geom from an iterable of component objects that match its structure.

All objects in components must have the same GeoInterface.trait.

Usually used in combination with flatten.

source
GeometryOpsCore.reconstruct_table Method
julia
reconstruct_table(input, geometry_column_names, geometry_columns, other_column_names, args...; kwargs...)

Reconstruct a table from the given input, geometry column names, geometry columns, and other column names.

Any function that defines reconstruct_table must also define used_reconstruct_table_kwargs.

The input must be a table.

The function should return a best-effort attempt at a table of the same type as the input, with the new geometry column(s) and other columns.

The fallback implementation invokes Tables.materializer. But if you want to be efficient and pass e.g. arbitrary kwargs to the materializer, or materialize in a different way, you can do so by overloading this function for your desired input type.

This is "semi-public" API and while it may add optional arguments, it will not add new required positional arguments. All implementations must allow arbitrary kwargs to pass through and harvest what they need.

source
GeometryOpsCore.unwrap Function
julia
unwrap(target::Type{<:AbstractTrait}, obj)
unwrap(f, target::Type{<:AbstractTrait}, obj)

Unwrap the object to vectors, down to the target trait.

If f is passed in it will be applied to the target geometries as they are found.

source
GeometryOpsCore.used_reconstruct_table_kwargs Method
julia
used_reconstruct_table_kwargs(input)

Return a tuple of the kwargs that should be passed to reconstruct_table for the given input.

This is "semi-public" API, and required for any input type that defines reconstruct_table.

source

Submodules ​

These are internal submodules of GeometryOps. They are not part of the public API, so anything here may change or disappear in a patch release, but they are documented for reference (and because GeometryOps' own source code links to them).

UnitSpherical ​

GeometryOps.UnitSpherical.ArcIntersectionResult Type
julia
ArcIntersectionResult{T}

Result of computing the intersection between two great circle arcs.

Fields:

  • type::ArcIntersectionType: The type of intersection

  • points::Vector{UnitSphericalPoint{T}}: The intersection point(s)

  • fracs::Vector{Tuple{T,T}}: For each intersection point, the fractional positions (α, β) where α is the position along the first arc and β is the position along the second arc. Both are in [0, 1], where 0 is the start point and 1 is the end point.

source
GeometryOps.UnitSpherical.ArcIntersectionType Type
julia
ArcIntersectionType

Enumeration of the types of arc intersections:

  • arc_cross: The arcs cross at a single point in their interiors

  • arc_hinge: The arcs share exactly one endpoint

  • arc_overlap: The arcs are collinear and overlap in a segment

  • arc_disjoint: The arcs do not intersect

source
GeometryOps.UnitSpherical.GeographicFromUnitSphere Type
julia
GeographicFromUnitSphere()

A transformation that converts a UnitSphericalPoint in ℝ³ to a 2-tuple geographic point (longitude, latitude), in degrees.

Accepts any 3-element vector, but the input is assumed to be on the unit sphere.

Examples

julia
julia> using GeometryOps.UnitSpherical

julia> GeographicFromUnitSphere()(UnitSphericalPoint(0.5, 0.5, 1/√(2)))
(45.0, 44.99999999999999)

(the inaccuracy is due to the precision of the atan function)

source
GeometryOps.UnitSpherical.SphericalCap Type
julia
SphericalCap{T}
SphericalCap(point::UnitSphericalPoint{T}, radius::T)

A spherical cap represents a section of a unit sphere about some point, bounded by a radius. It is defined by a center point on the unit sphere and a radius (in radians).

source
GeometryOps.UnitSpherical.SphericalCap Method

Get the circumcenter of the triangle (a, b, c) on the unit sphere. Returns a normalized 3-vector.

source
GeometryOps.UnitSpherical.UnitSphereFromGeographic Type
julia
UnitSphereFromGeographic()

A transformation that converts a geographic point (latitude, longitude) to a [UnitSphericalPoint] in ℝ³.

Accepts any GeoInterface-compatible point.

Examples

julia
julia> import GeoInterface as GI; using GeometryOps.UnitSpherical

julia> UnitSphereFromGeographic()(GI.Point(45, 45))
3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
 0.5000000000000001
 0.5000000000000001
 0.7071067811865476
julia
julia> using GeometryOps.UnitSpherical

julia> UnitSphereFromGeographic()((45, 45))
3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
 0.5000000000000001
 0.5000000000000001
 0.7071067811865476
source
GeometryOps.UnitSpherical.UnitSphericalPoint Type
julia
UnitSphericalPoint(v)

A unit spherical point, i.e., point living on the 2-sphere (𝕊²), represented as Cartesian coordinates in ℝ³.

This currently has no support for heights, only going from lat long to spherical and back again.

Examples

julia
julia> using GeometryOps.UnitSpherical

julia> UnitSphericalPoint(1, 0, 0)
3-element UnitSphericalPoint{Int64} with indices SOneTo(3):
 1
 0
 0
source
Base.convert Method
julia
convert(::Type{<:Extents.Extent}, cap::SphericalCap) -> Extents.Extent{(:X, :Y, :Z)}

Convert cap to an outward-rounded Cartesian bounding box containing every unit-sphere point in the cap. The result is intended for Cartesian spatial indexes over unit-spherical data.

Extent extraction conventionally means that an object has a Cartesian extent. A SphericalCap is instead a spherical query object, so conversion makes this boundary crossing explicit.

Valid caps have radii in [0, π]. A radius at least π, or a non-finite or negative radius, conservatively returns the whole unit-sphere box.

The abstract target Extents.Extent and compatible concrete targets are accepted. Incompatible concrete targets throw an ArgumentError.

source
Extents.contains Method
julia
Extents.contains(big::SphericalCap, small::SphericalCap; strict=false)

Whether the closed cap big contains all of small.

source
Extents.grow Method
julia
Extents.grow(cap::SphericalCap, factor::Real) -> SphericalCap

Grow the cap's angular diameter by factor on each side.

source
Extents.intersects Method
julia
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  , 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.

source
Extents.intersects Method
julia
Extents.intersects(x::SphericalCap, y::SphericalCap)

Whether the two closed caps intersect, including at tangency.

source
Extents.union Method
julia
Extents.union(x::SphericalCap, y::SphericalCap; strict=false) -> SphericalCap

Return a cap containing both x and y.

If either input already contains the other it is returned unchanged. Otherwise, the centre lies between the input caps' centres, and the radius is rounded outward to preserve coverage.

strict is accepted for consistency with the Extents API and has no effect because both inputs use the same domain.

source
Extents.within Method
julia
Extents.within(small::SphericalCap, big::SphericalCap; strict=false)

Whether the closed cap small is within big.

source
GeometryOps.UnitSpherical._arc_fraction Method
julia
_arc_fraction(p, a, b) -> T

Compute the fractional position of point p along the arc from a to b.

Returns a value in [0, 1] where 0 corresponds to a and 1 corresponds to b. Uses spherical distance to compute the fraction.

source
GeometryOps.UnitSpherical._find_collinear_arc_intersection Method
julia
_find_collinear_arc_intersection(a1, b1, a2, b2, T) -> ArcIntersectionResult{T}

Find the intersection of two collinear arcs (arcs on the same great circle).

Returns an ArcIntersectionResult with type arc_overlap if the arcs overlap, or arc_disjoint if they don't.

source
GeometryOps.UnitSpherical._make_hinge_result Method
julia
_make_hinge_result(T, point, α, β) -> ArcIntersectionResult{T}

Create an ArcIntersectionResult for a hinge intersection.

source
GeometryOps.UnitSpherical._points_equal Method
julia
_points_equal(a::UnitSphericalPoint, b::UnitSphericalPoint, tol) -> Bool

Check if two points are equal within a tolerance.

source
GeometryOps.UnitSpherical.circumcenter_on_unit_sphere Method
julia
circumcenter_on_unit_sphere(a, b, c)

Return the center of the circle passing through the three UnitSphericalPoints a, b and c, as a UnitSphericalPoint.

Of the two antipodal circumcenters, this returns the one on the same hemisphere as the input points, i.e. the center of the smaller of the two circles they bound.

source
GeometryOps.UnitSpherical.point_on_spherical_arc Method
julia
point_on_spherical_arc(p::UnitSphericalPoint, a::UnitSphericalPoint, b::UnitSphericalPoint) -> Bool

Check if point p lies on the great circle arc from a to b.

The arc is the shorter path along the great circle connecting a and b. Returns true if p is on the arc (including endpoints), false otherwise.

Examples

julia
using GeometryOps.UnitSpherical: UnitSphericalPoint, point_on_spherical_arc
a = UnitSphericalPoint(1.0, 0.0, 0.0)
b = UnitSphericalPoint(0.0, 1.0, 0.0)
mid = UnitSphericalPoint(1/√2, 1/√2, 0.0)
point_on_spherical_arc(mid, a, b)
# output
true
source
GeometryOps.UnitSpherical.slerp Method
julia
slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01::Number)

Interpolate between a and b, at a proportion i01 between 0 and 1 along the path from a to b.

Uses the tangent-vector form cos(r)·a + sin(r)·dir — where r = i01 · spherical_distance(a, b) and dir = normalize(robust_cross_product(a, b) × a) is the unit tangent at a pointing toward b. This avoids the 1/sin(Ω) divisor of the classic sin((1-t)Ω)/sin(Ω) · a + sin(tΩ)/sin(Ω) · b formulation, which collapses for near- and exactly-antipodal inputs. Adapted from Google's S2 geometry library (see S2::Interpolate and S2::GetPointOnLine).

For exactly antipodal a and b the great circle is mathematically ambiguous; robust_cross_product returns a deterministic perpendicular via its symbolic-perturbation branch, so the result is still a well-defined unit vector on some great circle through both points.

Examples

julia
julia> using GeometryOps.UnitSpherical

julia> slerp(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(0, 1, 0), 0.5)
3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
 0.7071067811865476
 0.7071067811865475
 0.0
source
GeometryOps.UnitSpherical.spherical_arc_extent Method
julia
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 UnitSphericalPoints, 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 is used.

Example

julia
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
source
GeometryOps.UnitSpherical.spherical_arc_intersection Method
julia
spherical_arc_intersection(a1, b1, a2, b2) -> ArcIntersectionResult

Compute the intersection between two great circle arcs on the unit sphere.

The first arc goes from a1 to b1, and the second arc goes from a2 to b2. All points should be UnitSphericalPoint instances.

Returns an ArcIntersectionResult containing:

  • The type of intersection (cross, hinge, overlap, or disjoint)

  • The intersection point(s) if they exist

  • The fractional positions along each arc for each intersection point

Examples

julia
# Two arcs that cross
a1 = UnitSphereFromGeographic()((-45.0, 0.0))
b1 = UnitSphereFromGeographic()((45.0, 0.0))
a2 = UnitSphereFromGeographic()((0.0, -45.0))
b2 = UnitSphereFromGeographic()((0.0, 45.0))
result = spherical_arc_intersection(a1, b1, a2, b2)
# result.type == arc_cross, with one intersection point at (0°, 0°)
source
GeometryOps.UnitSpherical.spherical_distance Method
julia
spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint)

Compute the great-circle distance (central angle, in radians) between two points on the unit sphere. Returns a Number whose type follows from the inputs.

Uses the atan2(‖x × y‖, x · y) form, which is numerically stable across the full [0, π] range — unlike acos(x · y), which loses precision for nearly-identical points. Adapted from Google's S2 geometry library (see Vector3::Angle).

Extended help

Doctests

julia
julia> using GeometryOps.UnitSpherical

julia> spherical_distance(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(0, 1, 0))
1.5707963267948966
julia
julia> using GeometryOps.UnitSpherical

julia> spherical_distance(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(1, 0, 0))
0.0
source
GeometryOps.UnitSpherical.spherical_exterior_anchor Method
julia
spherical_exterior_anchor(pts, n) -> Union{UnitSphericalPoint{Float64}, Nothing}

A reference point exterior BY DEFINITION of the enclosed-region semantics of the ring pts[1:n]: the antipode of the ring's normalized vertex mass (the sum of the unit vertex directions). For any ring whose enclosed region is meaningfully smaller than a hemisphere, the vertex mass points into the cap the vertices bound, so its antipode lies in the larger — exterior — region.

Returns nothing when the mass norm is tiny (below 1e-6 per vertex): near-hemisphere or vertex-symmetric rings, whose vertices spread over a near-great circle. There the enclosed/complement distinction is itself near-degenerate (the turning-angle winding tolerance already treats exact hemispheres permissively — see _ring_is_ccw), so callers fall back to the winding-consistent wedge bootstrap of spherical_ring_contains.

source
GeometryOps.UnitSpherical.spherical_orient Method
julia
spherical_orient(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint) -> Int

Determine the orientation of point c with respect to the great circle arc from a to b.

Returns:

  • 1 if c is to the left of the arc (counter-clockwise)

  • -1 if c is to the right of the arc (clockwise)

  • 0 if c is on the great circle (collinear)

Uses robust_cross_product for numerical stability with nearly identical or antipodal points.

Examples

julia
using GeometryOps.UnitSpherical: UnitSphericalPoint, spherical_orient
a = UnitSphericalPoint(1.0, 0.0, 0.0)
b = UnitSphericalPoint(0.0, 1.0, 0.0)
c = UnitSphericalPoint(0.0, 0.0, 1.0)
spherical_orient(a, b, c)
# output
1

Extended help

Why this does not simply call robust_cross_product

The common path uses the unnormalized cross(a - b, a + b): orientation needs only the sign of its dot product with c. The squared degeneracy test avoids normalizing the cross product or taking a square root. Rounding at the boundary may change 0 to a sign, but cannot flip +1 to -1.

For nearly equal or antipodal a and b, the cross-product direction becomes unstable. Those cases fall back to robust_cross_product.

source
GeometryOps.UnitSpherical.spherical_ring_contains Method
julia
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.

  • on_arc(q, a, b)::Bool: boundary membership; default point_on_spherical_arc. 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.

source
GeometryOps.UnitSpherical.spherical_ring_encloses Method
julia
spherical_ring_encloses(pts, n, q;
    anchor, orient, on_arc, proper_crossing) -> Union{Bool, Nothing}

Whether q lies in the region ENCLOSED by the ring pts[1:n] (the closing edge pts[n] → pts[1] is implied; boundary points count as enclosed): even-odd crossing parity of the arc from q to a reference point that is exterior by definition of the enclosed-region semantics — anchor, by default the antipode of the normalized vertex mass (spherical_exterior_anchor).

Winding-independent, like spherical_ring_contains composed with a winding test — but where that composition bootstraps the interior from a local wedge at one edge and a global turning-angle sum, both of which a ring that self-intersects on the sphere defeats (a figure-eight's lobes cancel the turning angle, and the wedge answer is anchored to whichever lobe hosts the edge — S2's forced-through behavior, globally inverted on real data), parity from a definitionally exterior point degrades to even-odd semantics: both lobes enclosed, the far side out.

Returns nothing — callers fall back conservatively — when:

  • anchor === nothing (degenerate vertex mass, see spherical_exterior_anchor);

  • q is (nearly) antipodal to the anchor (the test arc is ill-defined: q sits at the center of the vertex mass);

  • the anchor lies exactly ON a ring edge (the test arc ends on the ring); or

  • proper_crossing reports a crossing as too close to call (-1; never with exact injected predicates).

The orient/on_arc/proper_crossing predicates are injectable exactly as in spherical_ring_contains; on_test_arc(v, a, b) decides whether a point already known to lie on the great circle of (a, b) lies on the closed minor arc (the vertex-grazing resolution below — exact callers inject their span test).

source
GeometryOps.UnitSpherical.to_unit_spherical_points Method
julia
to_unit_spherical_points(ring) -> Vector{<:UnitSphericalPoint}

Convert a ring (linear ring or any GeoInterface point iterator) to a vector of UnitSphericalPoints, treating geographic input as (longitude, latitude). UnitSphericalPoints pass through unchanged.

source

RobustCrossProduct ​

GeometryOps.UnitSpherical.RobustCrossProduct.ensureNormalizable Method
julia
ensureNormalizable(p::AbstractVector)

Scales a 3-vector as necessary to ensure that the result can be normalized without loss of precision due to floating-point underflow.

This matches S2's EnsureNormalizable function.

source
GeometryOps.UnitSpherical.RobustCrossProduct.exact_cross_product Method
julia
exact_cross_product(a::AbstractVector, b::AbstractVector)

Compute the cross product using arbitrary precision arithmetic. This is used when standard floating-point arithmetic is not accurate enough.

This matches S2's ExactCrossProd function, first trying higher precision if available, then exact arithmetic, then symbolic perturbation.

source
GeometryOps.UnitSpherical.RobustCrossProduct.isNormalizable Method
julia
isNormalizable(v::AbstractVector)

Returns true if the given vector's magnitude is large enough such that the angle to another vector of the same magnitude can be measured using angle calculations without loss of precision due to floating-point underflow.

This matches S2's IsNormalizable function.

source
GeometryOps.UnitSpherical.RobustCrossProduct.isUnitLength Method
julia
isUnitLength(v::AbstractVector)

Check if a vector has unit length within a small tolerance.

Returns true if the vector's magnitude is approximately 1.0.

The tolerance adapts to the element type of the vector to handle both Float64 (high precision) and Float32 (lower precision, e.g., from GeoJSON) inputs. Uses 16 * eps(T) which provides adequate margin for unit vectors constructed from trigonometric functions while still being tight enough to catch errors.

source
GeometryOps.UnitSpherical.RobustCrossProduct.isless_vector Method
julia
isless_vector(a::AbstractVector, b::AbstractVector)

Lexicographic comparison of vectors. This is used to establish a consistent order for symbolic perturbations.

Returns true if a comes before b in lexicographic order.

source
GeometryOps.UnitSpherical.RobustCrossProduct.min_stable_norm Method
julia
min_stable_norm(::Type{T})

Return the smallest stable-cross-product norm that meets ROBUST_CROSS_PROD_ERROR in precision T.

This is S2's kMinNorm: s2edge_crossings.cc#L129-L131.

source
GeometryOps.UnitSpherical.RobustCrossProduct.normalizableFromExact Method
julia
normalizableFromExact(xf::Vector{BigFloat})

Converts a BigFloat vector to a double-precision vector, scaling the result as necessary to ensure that the result can be normalized without loss of precision due to floating-point underflow.

This matches S2's NormalizableFromExact function.

source
GeometryOps.UnitSpherical.RobustCrossProduct.normalization_needed Method
julia
normalization_needed(v::AbstractVector)

Determines if a vector's magnitude is too small for reliable normalization. Returns true if the vector needs special handling to avoid numerical issues.

This is essentially the opposite of isNormalizable.

source
GeometryOps.UnitSpherical.RobustCrossProduct.robust_cross_product Method
julia
robust_cross_product(a::AbstractVector, b::AbstractVector)

Compute a robust version of a × b (cross product) for unit vectors.

This method handles the case where a and b are very close together or antipodal by computing a stable perpendicular to both points.

The implementation follows Google's S2 Geometry Library to ensure numerical stability even in difficult cases.

Returns a unit-length vector that is perpendicular to both input vectors.

Examples

julia
using GeometryOps.UnitSpherical: UnitSphericalPoint, robust_cross_product
a = UnitSphericalPoint(1, 0, 0);
b = UnitSphericalPoint(0, 1, 0);
result = robust_cross_product(a, b)
isapprox(result, UnitSphericalPoint(0, 0, 1))
# output
true
source
GeometryOps.UnitSpherical.RobustCrossProduct.stable_cross_product Method
julia
getStableCrossProd(a::AbstractVector, b::AbstractVector)

Computes a numerically stable cross product between unit vectors.

This implements the algorithm from S2's GetStableCrossProd function, computing (a-b)×(a+b) which yields better numerical stability when the vectors are nearly identical.

Returns a tuple of (result, success) where:

  • result is the computed cross product vector (not normalized)

  • success is a boolean indicating if the computation was sufficiently accurate

source
GeometryOps.UnitSpherical.RobustCrossProduct.symbolic_cross_product Method
julia
symbolic_cross_product(a::AbstractVector, b::AbstractVector)

Compute a symbolic cross product when exact arithmetic yields zero. This implements the symbolic perturbation model used in S2 geometry.

Returns a vector that is the symbolic cross product.

source
GeometryOps.UnitSpherical.RobustCrossProduct.symbolic_cross_product_sorted Method
julia
symbolic_cross_product_sorted(a::AbstractVector, b::AbstractVector)

Helper function to compute the symbolic cross product when points are collinear. Assumes that a < b lexicographically.

This implements the symbolic perturbation model described in S2 geometry.

source

SpatialTreeInterface ​

GeometryOps.SpatialTreeInterface.FlatNoTree Type
julia
FlatNoTree(iterable_of_geoms_or_extents)

Represents a flat collection with no tree structure, i.e., a brute force search. This is cost free, so particularly useful when you don't want to build a tree!

source
GeometryOps.SpatialTreeInterface.child_indices_extents Method
julia
child_indices_extents(node)

Return an iterator over the indices and extents of the children of a node.

Each value of the iterator should take the form (i, extent).

This can only be invoked on leaf nodes!

source
GeometryOps.SpatialTreeInterface.depth_first_search Method
julia
depth_first_search(f, predicate, tree)

Call f(i) for each index i in the tree that satisfies predicate(extent(i)).

This is generic to anything that implements the SpatialTreeInterface, particularly the methods isleaf, getchild, and child_indices_extents.

Example

julia
using Extents, GeometryOps.SpatialTreeInterface, SortTileRecursiveTree
# Construct a tree of extents (in this case, just a grid)
xs, ys = 1:10, 1:10
extents_vec = vec([Extents.Extent(X = (x, x+1), Y = (y, y+1)) for x in xs, y in ys])

# Construct a tree of extents - this is an STRtree,
# but works on any SpatialTreeInterface-compatible tree.
extents_tree = STRtree(extents_vec)

# Count the number of extents that intersect the extent (1,1) x (2,2)
target_extent = Extents.Extent(X = (1,2), Y = (1,2))
count = 0
SpatialTreeInterface.depth_first_search(
    (i::Int -> global count += 1), # `f`
    Base.Fix1(Extents.intersects, target_extent), 
    extents_tree
)
count
# output
4
source
GeometryOps.SpatialTreeInterface.dual_depth_first_search Method
julia
dual_depth_first_search(f, predicate, tree1, tree2)

Executes a dual depth-first search over two trees, descending into the children of nodes i and j when predicate(node_extent(i), node_extent(j)) is true, and pruning that branch when predicate(node_extent(i), node_extent(j)) is false.

Finally, calls f(i1, i2) for each leaf-level index i1::Int in tree1 and i2::Int in tree2 that satisfies predicate(extent(i1), extent(i2)).

Here, f(i1::Int, i2::Int) may be any function that takes two integers as arguments.

It may optionally return an Action to alter the control flow of the Action(:full_return, true). Return Action(:full_return, true) from this function and break out of the recursion.

This is generic to anything that implements the SpatialTreeInterface, particularly the methods isleaf, getchild, node_extent and child_indices_extents.

Each visited node's extent is computed once and carried into the recursion. Trees that derive their extents rather than storing them should also define node_extent_is_expensive.

source
GeometryOps.SpatialTreeInterface.getchild Function
julia
getchild(node)
getchild(node, i)

Accessor function to get the children of a node.

If invoked as getchild(node), return an iterator over all the children of a node. This may be lazy, like a Base.Generator, or it may be materialized.

If invoked as getchild(node, i), return the i-th child of a node.

source
GeometryOps.SpatialTreeInterface.getchild Method
julia
getchild(node, i)

Return the i-th child of a node.

source
GeometryOps.SpatialTreeInterface.isleaf Method
julia
isleaf(node)

Return true if the node is a leaf node, i.e., there are no "children" below it. getchild should still work on leaf nodes, though, returning an iterator over the extents stored in the node - and similarly for getnodes.

source
GeometryOps.SpatialTreeInterface.isspatialtree Method
julia
isspatialtree(tree)::Bool

Return true if the object is a spatial tree, false otherwise.

Implementation notes

For type stability, if your spatial tree type is MyTree, you should define isspatialtree(::Type{MyTree}) = true, and isspatialtree(::MyTree) will forward to that method automatically.

source
GeometryOps.SpatialTreeInterface.nchild Method
julia
nchild(node)

Return the number of children of a node.

source
GeometryOps.SpatialTreeInterface.node_extent Method
julia
node_extent(node)

Return the extent like object of the node. Falls back to GI.extent by default, which falls back to Extents.extent.

Generally, defining Extents.extent(node) is sufficient here, and you won't need to define this

The reason we don't use that directly is to give users of this interface a way to define bounding boxes that are not extents, like spherical caps and other such things.

source
GeometryOps.SpatialTreeInterface.node_extent_is_expensive Method
julia
node_extent_is_expensive(node)::Bool

Return true if node_extent computes the node's extent instead of reading one the node already stores. Defaults to false.

When true, dual_depth_first_search caches a node's child extents rather than re-deriving them once per opposing child, at the cost of a small vector per visited node.

Implementation notes

Define this on the type - node_extent_is_expensive(::Type{MyNode}) = true - so that it is known at compile time; node_extent_is_expensive(::MyNode) forwards there automatically.

source
GeometryOps.SpatialTreeInterface.query Method
julia
query(tree, predicate)

Return a sorted list of indices of the tree that satisfy the predicate.

source
GeometryOps.SpatialTreeInterface.sanitize_predicate Method
julia
sanitize_predicate(pred)

Convert a predicate to a function that returns a Boolean.

If pred is an Extent, convert it to a function that returns a Boolean by intersecting with the extent. If pred is a geometry, convert it to an extent first, then wrap in Extents.intersects.

Otherwise, return the predicate unchanged.

Users and developers may overload this function to provide custom behaviour when something is passed in.

source
GeometryOps.SpatialTreeInterface.spatialtree Method
julia
spatialtree([manifold::Manifold], geometries)

Build a default packed STR RTree from an iterable of geometries, over their extents on manifold (Planar() if not given).

missing and nothing entries are skipped, as are empty geometries. Queries return indices into geometries regardless, so the skipped entries do not shift the answers. If no valid geometries remain, return nothing.

source

FlexibleRTrees ​

GeometryOps.FlexibleRTrees.BulkLoadAlgorithm Type
julia
BulkLoadAlgorithm

Supertype for the algorithms that decide the leaf order of an RTree. Packing is always "union consecutive runs of nodecapacity, bottom-up"; the algorithm only chooses the order, via a loadorder method.

source
GeometryOps.FlexibleRTrees.HPR Type
julia
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.

source
GeometryOps.FlexibleRTrees.RTree Type
julia
RTree(algorithm::BulkLoadAlgorithm, data; nodecapacity = 16, extents = nothing, indices = nothing)

A packed R-tree over the extents of data (anything GI.extent accepts — geometries, or Extents.Extents themselves), of any dimensionality, bulk loaded in the order chosen by algorithm.

Pass a vector as indices to index only data[indices], leaving the rest of data out of the tree — for a collection only part of which can be indexed, say one with missing entries. data stays whole, so queries still report positions in it.

Pass a vector as extents to index by precomputed extents instead of GI.extent — for payload elements that carry no extent of their own, or extents computed in another coordinate space. One per indexed element, in order: per element of data normally, per element of indices alongside indices. The tree takes ownership of the vector (Unsorted aliases it as the leaf level rather than copying).

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 data, which the tree keeps as tree.data so hits map straight back to elements wherever the tree travels.

source
GeometryOps.FlexibleRTrees.RTree Method
julia
RTree(m::Manifold, algorithm::BulkLoadAlgorithm, data; nodecapacity = 16, indices = nothing)

Build the tree over each element's extent on the manifold m, via Extents.extent(m, x). On Spherical() the leaves are the 3D Cartesian boxes of the elements as regions on the unit sphere, covering the arc bulge and enclosed poles that vertex extents miss.

source
GeometryOps.FlexibleRTrees.RTreeNode Type
julia
RTreeNode{T, E}

A cursor into one node of an RTree: 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.

source
GeometryOps.FlexibleRTrees.STR Type
julia
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.

source
GeometryOps.FlexibleRTrees.Unsorted Type
julia
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.

source
GeometryOps.FlexibleRTrees.hilbert_key Method
julia
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.

source
GeometryOps.FlexibleRTrees.loadorder Method
julia
loadorder(algorithm, extents::Vector{<:Extents.Extent}, nodecapacity)

The permutation (an AbstractVector{Int}) in which algorithm packs extents into leaves. Implement this for a new BulkLoadAlgorithm subtype to plug in another ordering. Return Base.OneTo for the identity — the constructor then skips the reorder copy and aliases extents as the leaf level.

source
GeometryOps.FlexibleRTrees.query Method
julia
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.

source

NaturalIndexing ​

GeometryOps.NaturalIndexing.NaturalIndex Type
julia
NaturalIndex{E <: Extents.Extent}

A natural tree index. Stored in a vector in NaturalIndex.

  • nodecapacity is the "spread", number of children per node

  • extent is the extent of the tree

  • levels is a vector of NaturalLevels

source
GeometryOps.NaturalIndexing.NaturalIndex Method
julia
NaturalIndex(m::Manifold, geoms; nodecapacity = 32)

Index each geometry's extent on the manifold m, via Extents.extent(m, geom). On Spherical() the leaf extents are the 3D Cartesian boxes of the geometries as regions on the unit sphere.

source
GeometryOps.NaturalIndexing.NaturalIndexNode Type
julia
NaturalIndexNode{E <: Extents.Extent}

A reference to a node in the natural tree. Kind of like a tree cursor.

  • parent_index is a pointer to the parent index

  • level is the level of the node in the tree

  • index is the index of the node in the level

  • extent is the extent of the node

source
GeometryOps.NaturalIndexing.NaturalLevel Type
julia
NaturalLevel{E <: Extents.Extent}

A level in the natural tree. Stored in a vector in NaturalIndex.

  • extents is a vector of extents of the children of the node
source
GeometryOps.NaturalIndexing.NaturallyIndexedRing Type
julia
NaturallyIndexedRing(points; nodecapacity = 32)

A linear ring that contains a natural index.

Warning

This will be removed in favour of prepared geometry - the idea here is just to test what interface works best to store things in.

source
GeometryOps.NaturalIndexing._number_of_levels Method
julia
_number_of_levels(nodecapacity::Int, ngeoms::Int)

Calculate the number of levels in a natural tree for a given number of geometries and node capacity.

How this works

The number of keys in a level is given by ngeoms / nodecapacity ^ level.

The number of levels is the smallest integer such that the number of keys in the last level is 1. So it goes - if that makes sense.

source

LoopStateMachine ​

GeometryOps.LoopStateMachine Module
julia
LoopStateMachine

Utilities for returning state from functions that run inside a loop.

This is used in e.g clipping, where we may need to break or transition states.

The main entry point is to return an Action from a function that is wrapped in a @controlflow f(...) macro in a loop. When a known Action (currently, :continue, :break, :return, or :full_return actions) is returned, it is processed by the @controlflow macro, which allows the function to break out of the loop early, continue to the next iteration, or return a value, basically a way to provoke syntactic behaviour from a function called from a inside a loop, where you do not have access to that loop.

Example

julia
source
GeometryOps.LoopStateMachine.Action Type
julia
Action(name::Symbol, [x])

Create an Action with the name name and optional contents x.

Actions are returned from functions wrapped in a @controlflow macro, which does something based on the return value of that function if it is an Action.

Available actions

  • :continue: continue to the next iteration of the loop. This is the continue keyword in Julia. The contents of the action are not used.

  • :break: break out of the loop. This is the break keyword in Julia. The contents of the action are not used.

  • :return: cause the function executing the loop to return with the wrapped value.

  • :full_return: cause the function executing the loop to return Action(:full_return, x). This is very useful to terminate recursive funtions, like tree queries terminating after you have found a single intersecting segment.

source
GeometryOps.LoopStateMachine.@controlflow Macro
julia
@controlflow f(...)

Process the result of f(...) and return the result if it's not an Action(@ref LoopStateMachine.Action).

If it is an Action, then process it according to the following rules, and throw an error if it's not recognized. :continue, :break, :return, or :full_return are valid actions.

  • :continue: continue to the next iteration of the loop. This is the continue keyword in Julia. The contents of the action are not used.

  • :break: break out of the loop. This is the break keyword in Julia. The contents of the action are not used.

  • :return: cause the function executing the loop to return with the wrapped value.

  • :full_return: cause the function executing the loop to return Action(:full_return, x). This is very useful to terminate recursive funtions, like tree queries terminating after you have found a single intersecting segment.

Warning

Only use this inside a loop, otherwise you'll get a syntax error, especially if you use :continue or :break.

Examples

source

  1. Chaikin, G. An algorithm for high speed curve generation. Computer Graphics and Image Processing 3 (1974), 346-349 ↩︎

  2. K. Hormann and N. Sukumar. Generalized Barycentric Coordinates in Computer Graphics and Computational Mechanics. Taylor & Fancis, CRC Press, 2017. ↩︎