GeometryOps.jl

https://juliageo.org/GeometryOps.jl/stable/
Admin
GeometryOps.jl is a pure Julia package for fast geometric calculations and operations on 2D and 2.5D...
Tokens:194,589
Snippets:1,583
Trust Score:7
Update:1 month ago
Show doc for...
Tokens:
Raw
### Install required packages

Source: https://juliageo.org/GeometryOps.jl/stable/tutorials/creating_geometry

Install the necessary geospatial and plotting packages for this tutorial.

```julia
using Pkg
Pkg.add(["GeoInterface", "GeometryOps", 
        "CoordinateTransformations",
        "CairoMakie", "GeoMakie"])
```

--------------------------------

### Equals Implementation and Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/equals

Documentation string and jldoctest example for the equals function.

```julia
"""
    equals(geom1, geom2)::Bool

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

# Examples
```jldoctest
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)
```
```

--------------------------------

### Voronoi Output Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/voronoi

Example output format for the voronoi function.

```julia
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)])])
```

--------------------------------

### Visualize geometry touching

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/touches

Example setup using Makie to visualize two lines that touch at their endpoints.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

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

f, a, p = lines(GI.getpoint(l1))
lines!(GI.getpoint(l2))
f
```

--------------------------------

### Setup Imports

Source: https://juliageo.org/GeometryOps.jl/stable/source/GeometryOpsMakieExt/spherical_cap

Required modules for spherical cap geometry and Makie visualization.

```julia
using Makie
import GeometryOps as GO
import GeometryOps.UnitSpherical
import LinearAlgebra
import Makie.GeometryBasics
import StaticArrays
```

--------------------------------

### Spatial Tree Setup and Search

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/SpatialTreeInterface/depth_first_search

Basic setup for creating an STRtree and performing a depth-first search to count intersecting extents.

```julia
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])
```

```julia
extents_tree = STRtree(extents_vec)
```

```julia
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
```

--------------------------------

### Implementation and examples of covers

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/covers

The core implementation of covers using coveredby, including a doctest example.

```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

```jldoctest
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
```
"""
covers(g1, g2)::Bool = GeometryOps.coveredby(g2, g1)

"""
    covers(g1)

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

--------------------------------

### Visualize disjoint geometries

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/disjoint

Example setup using Makie to visualize two disjoint LineStrings.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)])
l2 = GI.LineString([(2.0, 0.0), (2.75, 0.0)])
f, a, p = lines(GI.getpoint(l1), color = :blue)
scatter!(GI.getpoint(l1), color = :blue)
lines!(GI.getpoint(l2), color = :orange)
scatter!(GI.getpoint(l2), color = :orange)
f
```

--------------------------------

### Threaded Transformation Setup

Source: https://juliageo.org/GeometryOps.jl/stable/source/GeometryOpsProjExt/reproject

Clones transformations and contexts for multi-threaded execution.

```julia
        proj_transforms = [Proj.Transformation(Proj.proj_clone(transform.pj)) for _ in 1:ntasks]
```

```julia
        contexts = [Proj.proj_context_clone(context) for _ in 1:ntasks]
```

```julia
        foreach(Proj.proj_assign_context, getproperty.(proj_transforms, :pj), contexts)

        results = if _is3d(geom)
            functors = TaskFunctors(WithXYZ.(proj_transforms))
            apply(functors, GI.PointTrait(), geom; kw1...)
        else
            functors = TaskFunctors(WithXY.(proj_transforms))
            apply(functors, GI.PointTrait(), geom; kw1...)
        end
```

--------------------------------

### Benchmark Setup for Segmentize

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/segmentize

Initializes the benchmark suite and prepares a rectangle geometry for comparison against LibGEOS.

```julia
using BenchmarkTools: BenchmarkGroup
using Chairmarks: @be
using Main: plot_trials
using CairoMakie

import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG

segmentize_suite = BenchmarkGroup(["title:Segmentize", "subtitle:Segmentize a rectangle"])

rectangle = GI.Wrappers.Polygon([[(0.0, 50.0), (7.071, 57.07), (0.0, 64.14), (-7.07, 57.07), (0.0, 50.0)]])
lg_rectangle = GI.convert(LG, rectangle)
```

--------------------------------

### Initialize Ring

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Starts a new ring structure with the current and next nodes.

```julia
        currentnode = firstnode
        ring = [currentnode, nextnode]
        push!(rings, ring)
```

--------------------------------

### Visualizing geometry coverage

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/covers

Setup using Makie to visualize a line and a point to demonstrate coverage.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

p1 = (0.0, 0.0)
p2 = (1.0, 1.0)
l1 = GI.Line([p1, p2])

f, a, p = lines(GI.getpoint(l1))
scatter!(p1, color = :red)
f
```

--------------------------------

### Intersection Points Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/intersection

Demonstrates calculating intersection points between two lines.

```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_points(line1, line2)
```

--------------------------------

### Module Setup and Constants

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct

Initializes the module dependencies and defines numerical error constants used for robust calculations.

```julia
module RobustCrossProduct

using ..UnitSpherical: UnitSphericalPoint
using StaticArrays
using LinearAlgebra

include("utils.jl")
```

```julia
const DBL_ERR = eps(Float64) / 2
```

```julia
const SQRT3 = sqrt(3.0)
```

```julia
const ROBUST_CROSS_PROD_ERROR = 6 * DBL_ERR
```

```julia
const HAS_LONG_DOUBLE = precision(Float64) < precision(BigFloat)
```

```julia
const EXACT_CROSS_PROD_ERROR = DBL_ERR

isDoubleFloatsAvailable(args...) = false
```

--------------------------------

### Construct TraitTarget instances

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/types/traittarget

Examples of constructing TraitTarget using trait instances.

```julia
TraitTarget(GI.PointTrait())
TraitTarget(GI.LineStringTrait(), GI.LinearRingTrait()) # and other traits as you may like
TraitTarget(TraitTarget(...))
```

--------------------------------

### Constructing SphericalCap instances

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/UnitSpherical/cap

Examples demonstrating how to create a SphericalCap from a point and radius, geographic coordinates, or three points.

```julia
using GeometryOps
using GeoInterface

# Create a spherical cap from a point and radius
point = UnitSphericalPoint(1.0, 0.0, 0.0)  # Point on the unit sphere
cap = SphericalCap(point, 0.5)  # Cap with radius 0.5 radians
```

```julia
# Create a spherical cap from geographic coordinates
lat, lon = 40.0, -74.0  # New York City (approximate)
point = GeoInterface.Point(lon, lat)
cap = SphericalCap(point, 0.1)  # Cap with radius ~0.1 radians
```

```julia
# Create a spherical cap from three points (circumcircle)
p1 = UnitSphericalPoint(1.0, 0.0, 0.0)
p2 = UnitSphericalPoint(0.0, 1.0, 0.0)
p3 = UnitSphericalPoint(0.0, 0.0, 1.0)
cap = SphericalCap(p1, p2, p3)
```

--------------------------------

### Basic Geometry Buffering

Source: https://juliageo.org/GeometryOps.jl/stable/source/GeometryOpsLibGEOSExt/buffer

Example of creating a 1-unit buffer around a point geometry.

```julia
point = GI.Point(0, 0)
buffered = GO.buffer(GO.GEOS(), point, 1.0)  # 1 unit buffer
```

--------------------------------

### Module Imports and Setup

Source: https://juliageo.org/GeometryOps.jl/stable/source/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt

Initializes the module by importing necessary packages and filtering out internal Julia functions to prevent precompilation warnings.

```julia
module GeometryOpsLibGEOSExt

import GeometryOps as GO, LibGEOS as LG
import GeoInterface as GI

import GeometryOps: GEOS, enforce, True, False, booltype, booltype, WithTrait, TraitTarget

using GeometryOps
```

```julia
for name in filter(!in((:var"#eval", :eval, :var"#include", :include)), names(GeometryOps))
    @eval import GeometryOps: $name
end
```

--------------------------------

### Slerp usage example

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/UnitSpherical/slerp

Example of interpolating between two points on a unit sphere at a 0.5 proportion.

```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
```

--------------------------------

### Visualize Depth First Search

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/SpatialTreeInterface/depth_first_search

An example demonstrating how to visualize the depth-first search process using Makie and record frames.

```julia
f, a, all_rects_plot = lines(rects_vec; linewidth = 0.5, axis = (; aspect = DataAspect()));
Legend(f[1, 2], [Makie.PolyElement(; color = :transparent, strokecolor = :blue, strokewidth = 0.5), Makie.PolyElement(; color = :green), Makie.PolyElement(; color = :red)], ["Grid", "Branch node", "Leaf node"])

target_extent = Extents.Extent(X = (5, 7), Y = (5, 7))
target_rect_plot = poly!(a, ext2rect(target_extent); color = (:blue, 0.5))

current_rect = Observable(ext2rect(GI.extent(rects_vec[1])))
current_color = Observable(:red)

current_rect_plot = poly!(a, current_rect; color = current_color, alpha = 0.3)


record(f, "./depth_first_search.mp4"; framerate = 1) do io

    function leaf_f(i)
        current_rect[] = ext2rect(GI.extent(rects_vec[i]))
        current_color[] = :red
        recordframe!(io)
        return nothing
    end

    function pred_f(ext, target)
        current_rect[] = ext2rect(ext)
        current_color[] = :green
        recordframe!(io)
        result = Extents.intersects(ext, target)
        if result
            current_color[] = :forestgreen
            recordframe!(io)
        end
        return result
    end

    SpatialTreeInterface.depth_first_search(leaf_f, Base.Fix2(pred_f, target_extent), extents_tree)
end
```

--------------------------------

### Example: Calculating and Plotting a Centroid

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/centroid

Demonstrates creating a concave polygon and calculating its centroid for visualization.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

cshape = GI.Polygon([[(0,0), (0,3), (3,3), (3,2), (1,2), (1,1), (3,1), (3,0), (0,0)]])
f, a, p = poly(collect(GI.getpoint(cshape)); axis = (; aspect = DataAspect()))
```

```julia
cent = GO.centroid(cshape)
scatter!(GI.x(cent), GI.y(cent), color = :red)
f
```

--------------------------------

### Linear Segmentization Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/segmentize

Demonstrates densifying a polygon using the linear method with a specified maximum distance.

```julia
import GeometryOps as GO, GeoInterface as GI
rectangle = GI.Wrappers.Polygon([[(0.0, 50.0), (7.071, 57.07), (0, 64.14), (-7.07, 57.07), (0.0, 50.0)]])
linear = GO.segmentize(rectangle; max_distance = 5)
collect(GI.getpoint(linear))
```

--------------------------------

### Arclength operation examples

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/types/operation

Demonstrates various ways to instantiate and call Arclength operations with different algorithms and manifolds.

```julia
Arclength()(geoms)
Arclength(Geodesic())(geoms)
Arclength(Proj())(geoms)
Arclength(Proj(Geodesic(; ...)))(geoms)
Arclength(Ericsson())(geoms) # more precise, goes wonky if any points in a triangle are antipodal
Arclength(LHuilier())(geoms) # less precise, does not go wonky on antipodal points
```

--------------------------------

### Define the GeometryOpsTGGeometryExt module

Source: https://juliageo.org/GeometryOps.jl/stable/source/GeometryOpsTGGeometryExt/GeometryOpsTGGeometryExt

Setup the module by importing necessary dependencies from GeometryOps and TGGeometry.

```julia
module GeometryOpsTGGeometryExt

using GeometryOps: TG
import GeometryOps as GO

using TGGeometry
```

--------------------------------

### Flip point order in geometry

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/apply

Full example demonstrating how to flip point order within a polygon geometry.

```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
```

--------------------------------

### Throw MissingKeywordInAlgorithmException

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/types/exceptions

Example of manually constructing and throwing the MissingKeywordInAlgorithmException.

```julia
throw(MissingKeywordInAlgorithmException(GEOS(; tokl = 1.0), my_function, :tol))
```

--------------------------------

### Constructing a spatial tree for depth-first search

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/SpatialTreeInterface/depth_first_search

Setup for a spatial tree using SortTileRecursiveTree and Makie for visualization.

```julia
using Makie
using Extents, GeometryOps.SpatialTreeInterface, SortTileRecursiveTree
import GeoInterface as GI

# Construct a tree of extents (in this case, just a grid)
xs, ys = 1.0:0.1:10.0, 1.0:0.1:10.0
rects_vec = vec([GI.Polygon([[(x, y), (x+0.1, y), (x+0.1, y+0.1), (x, y+0.1), (x, y)]]) for x in xs, y in ys])
ext2rect(ext) = Rect2f((ext.X[1], ext.Y[1]), (ext.X[2] - ext.X[1], ext.Y[2] - ext.Y[1]))
```

--------------------------------

### Perform Voronoi Tessellation with Clipping

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/voronoi

Examples demonstrating how to clip Voronoi tessellations using either a tuple of tuples or a tuple of vectors.

```julia
clip_poly2 = (clip_points, clip_order) # tuples
GO.voronoi(points; clip_polygon = clip_poly2, rng = rng)
```

```julia
clip_poly3 = (collect(clip_points), collect(clip_order)) # vectors
GO.voronoi(points; clip_polygon = clip_poly3, rng = rng)
```

--------------------------------

### Barycentric Interpolation Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/barycentric

Demonstrates polygon rendering and barycentric interpolation using the MeanValue method with CairoMakie.

```julia
import GeometryOps as GO, GeoInterface as GI
using CairoMakie # plotting
# Define a polygon
polygon_points = [
(0.03, 0.05, 0.00), (0.07, 0.04, 0.02), (0.10, 0.04, 0.04),
(0.14, 0.04, 0.06), (0.17, 0.07, 0.08), (0.20, 0.09, 0.10),
(0.22, 0.11, 0.12), (0.25, 0.11, 0.14), (0.27, 0.10, 0.16),
(0.30, 0.07, 0.18), (0.31, 0.04, 0.20), (0.34, 0.03, 0.22),
(0.37, 0.02, 0.24), (0.40, 0.03, 0.26), (0.42, 0.04, 0.28),
(0.44, 0.07, 0.30), (0.45, 0.10, 0.32), (0.46, 0.13, 0.34),
(0.46, 0.19, 0.36), (0.47, 0.26, 0.38), (0.47, 0.31, 0.40),
(0.47, 0.35, 0.42), (0.45, 0.37, 0.44), (0.41, 0.38, 0.46),
(0.38, 0.37, 0.48), (0.35, 0.36, 0.50), (0.32, 0.35, 0.52),
(0.30, 0.37, 0.54), (0.28, 0.39, 0.56), (0.25, 0.40, 0.58),
(0.23, 0.39, 0.60), (0.21, 0.37, 0.62), (0.21, 0.34, 0.64),
(0.23, 0.32, 0.66), (0.24, 0.29, 0.68), (0.27, 0.24, 0.70),
(0.29, 0.21, 0.72), (0.29, 0.18, 0.74), (0.26, 0.16, 0.76),
(0.24, 0.17, 0.78), (0.23, 0.19, 0.80), (0.24, 0.22, 0.82),
(0.24, 0.25, 0.84), (0.21, 0.26, 0.86), (0.17, 0.26, 0.88),
(0.12, 0.24, 0.90), (0.07, 0.20, 0.92), (0.03, 0.15, 0.94),
(0.01, 0.10, 0.97), (0.02, 0.07, 1.00)]
# Plot it!
# First, we'll plot the polygon using Makie's rendering:
f, a1, p1 = poly(
    Point2.(GO.forcexy(polygon_points));
    color = last.(polygon_points),
    colormap = cgrad(:jet, 18; categorical = true),
    axis = (; 
       type = Axis, aspect = DataAspect(), title = "Makie mesh based polygon rendering", subtitle = "CairoMakie"
    ),
    figure = (; size = (800, 400),)
)
hidedecorations!(a1)

ext = GO.Extents.Extent(X = (0, 0.5), Y = (0, 0.42))

a2 = Axis(
        f[1, 2],
        aspect = DataAspect(),
        title = "Barycentric coordinate based polygon rendering", subtitle = "GeometryOps",
        limits = (ext.X, ext.Y)
    )
hidedecorations!(a2)

p2box = poly!( # Now, we plot a cropping rectangle around the axis so we only show the polygon
    a2,
    GI.Polygon( # This is a rectangle with an internal hole shaped like the polygon.
        [
            Point2f[(ext.X[1], ext.Y[1]), (ext.X[2], ext.Y[1]), (ext.X[2], ext.Y[2]), (ext.X[1], ext.Y[2]), (ext.X[1], ext.Y[1])], # exterior
            reverse(Point2f.(GO.forcexy(polygon_points))) # hole
        ]
    ); color = :white, xautolimits = false, yautolimits = false
)
cb = Colorbar(f[2, :], p1.plots[1]; vertical = false, flipaxis = true)
# Finally, we perform barycentric interpolation on a grid,
xrange = LinRange(ext.X..., 400)
yrange = LinRange(ext.Y..., 400)
@time mean_values = GO.barycentric_interpolate.(
    (GO.MeanValue(),), # The barycentric coordinate algorithm (MeanValue is the only one for now)
    (GI.Polygon(GI.LinearRing.([polygon_points])),), # The polygon
    (last.(polygon_points,),),   # The values per polygon point - can be anything which supports addition and division
    tuple.(xrange, yrange')    # The points at which to interpolate
)
# and render!
hm = heatmap!(a2, xrange, yrange, mean_values; colormap = p1.colormap, colorrange = p1.plots[1].colorrange[], xautolimits = false, yautolimits = false)
translate!(hm, 0, 0, -1) # translate the heatmap behind the cropping polygon!
f # finally, display the figure
```

--------------------------------

### RobustCrossProduct Usage Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct

Demonstrates calculating the robust cross product of two UnitSphericalPoints and verifying the result.

```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
```

--------------------------------

### Calculate Geometry Coverage Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/coverage

Demonstrates setting up a rectangle and a cell to calculate the intersection area using GeometryOps.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

rect = GI.Polygon([[(-1,0), (-1,1), (1,1), (1,0), (-1,0)]])
cell = GI.Polygon([[(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)]])
xmin, xmax, ymin, ymax = 0, 2, 0, 2
f, a, p = poly(collect(GI.getpoint(cell)); axis = (; aspect = DataAspect()))
poly!(collect(GI.getpoint(rect)))
f
```

```julia
GO.coverage(rect, xmin, xmax, ymin, ymax)
```

--------------------------------

### Define and Plot Geometries

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/equals

Setup and visualization of two LineString geometries using Makie and GeoInterface.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

l1 = GI.LineString([(0.0, 0.0), (0.0, 10.0)])
l2 = GI.LineString([(0.0, -10.0), (0.0, 3.0)])
f, a, p = lines(GI.getpoint(l1), color = :blue)
scatter!(GI.getpoint(l1), color = :blue)
lines!(GI.getpoint(l2), color = :orange)
scatter!(GI.getpoint(l2), color = :orange)
f
```

--------------------------------

### Visualizing geometry intersection

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/intersects

Setup using Makie and GeoInterface to visualize two lines and verify their intersection.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

line1 = GI.Line([(124.584961,-12.768946), (126.738281,-17.224758)])
line2 = GI.Line([(123.354492,-15.961329), (127.22168,-14.008696)])
f, a, p = lines(GI.getpoint(line1))
lines!(GI.getpoint(line2))
f
```

--------------------------------

### Polygonize Matrix Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Basic usage of the polygonize function to convert a random matrix into a MultiPolygon.

```julia
using GeometryOps
A = rand(100, 100)
multipolygon = polygonize(>(0.5), A);
```

--------------------------------

### Flip coordinates using apply

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/apply

Example of flipping x and y coordinates of a geometry using the apply framework.

```julia
flipped_geom = GO.apply(GI.PointTrait(), geom) do p
    (GI.y(p), GI.x(p))
end
```

--------------------------------

### Example: Two Crossing Arcs

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/UnitSpherical/arc_intersection

Demonstrates usage of the intersection function for two arcs that cross.

```julia
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)
```

--------------------------------

### Minimum Bounding Circle Usage

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/minimum_bounding_circle

Examples of computing the minimum bounding circle from points or geometries.

```julia
import GeometryOps as GO, GeoInterface as GI
```

```julia
points = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)]
circle = GO.minimum_bounding_circle(points)
```

```julia
polygon = GI.Polygon([[(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]])
circle = GO.minimum_bounding_circle(polygon)
```

```julia
ring = GI.getexterior(circle)
```

--------------------------------

### Two-argument operation example

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/types/operation

Shows the structure for operations requiring two geometry arguments, such as set operations.

```julia
Union(intersection_alg(manifold); exact, target)(geom1, geom2)
```

--------------------------------

### Planar Geometry Smoothing Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/smooth

Demonstrates smoothing a LineString on a planar manifold with varying iteration counts.

```julia
using CairoMakie
import GeoInterface as GI, GeometryOps as GO

line = GI.LineString([(0.0, 0.0), (1.0, 1.0), (2.0, 0.0)])
smoothed = GO.smooth(line)
smoothed_2 = GO.smooth(line; iterations=2)

f, a, p = lines(line; label = "line")
lines!(a, smoothed; label = "smooth(line; iterations=1)")
lines!(a, smoothed_2; label = "smooth(line; iterations=2)")
axislegend(a)
f
```

--------------------------------

### Setup Makie peaks data

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Prepare a heatmap dataset using Makie's peaks function for polygonization.

```julia
using Makie, CairoMakie
using GeometryOps
n = 49
xs, ys = LinRange(-3, 3, n), LinRange(-3, 3, n)
zs = Makie.peaks(n)
z_max_value = maximum(abs.(extrema(zs)))
f, a, p = heatmap(
    xs, ys, zs;
    axis = (; aspect = DataAspect(), title = "Exact function")
)
cb = Colorbar(f[1, 2], p; label = "Z-value")
f
```

--------------------------------

### Initialize Buffers

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Sets up buffers for storing edges and rings.

```julia
rings = Vector{T}[]

strait = true
turning = false
```

--------------------------------

### Handle start of overlapping chain

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/clipping_processor

Initializes the start of an overlapping chain by determining the side and winding state.

```julia
            elseif !a_prev_is_b_prev && !a_prev_is_b_next
                b_side = a_next_is_b_prev ? b_next_side : b_prev_side
                start_chain_edge = b_side
                start_chain_idx = i
                same_winding = a_next_is_b_next
            end
        end
        a_prev = curr_pt
        curr_pt = a_next
        i = next_idx
    end
```

--------------------------------

### Update start of chain tags

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/clipping_processor

Updates PolyNode tags for the start of a chain with crossing and endpoint information.

```julia
                    start_pt = a_list[start_chain_idx]
                    a_list[start_chain_idx] = PolyNode(start_pt;
                        crossing = crossing,
                        endpoint = start_chain,
                    )
                    b_list[start_pt.neighbor] = PolyNode(b_list[start_pt.neighbor];
                        crossing = crossing,
                        endpoint = same_winding ? start_chain : end_chain,
                    )
                end
```

--------------------------------

### Handle start point removal with maximum distance

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/simplify

Removes the start point and inserts the point with the maximum remaining distance.

```julia
            if endpt_dist < max_dist
                insert!(results, searchsortedfirst(sorted_results, max_idx), max_idx)
                results[i+1] = results[2]
                sorted_results = @view results[2:i+1]
            end
        end
    end
    return points[sorted_results]
end
```

--------------------------------

### Handle start point removal with tolerance

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/simplify

Removes the start point and replaces it with the second point if the distance is within tolerance.

```julia
            if endpt_dist < max_tol
                results[i] = results[2]
                sorted_results = @view results[2:i]
            end
        else
```

--------------------------------

### Include Interface and Implementation Files

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/SpatialTreeInterface/SpatialTreeInterface

Loads the core interface definitions and implementation logic.

```julia
include("interface.jl")
include("implementations.jl")
```

--------------------------------

### Initialize Ring Creation Loop

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Prepares variables for the ring creation process.

```julia
while nkeys > 0
    found = false
    local firstnode, nextnodes, nodestatus
```

--------------------------------

### Get pixel bounds

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Retrieves the bounds for the current pixel.

```julia
x1, x2 = xs[i]
```

--------------------------------

### Get or write edge

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Retrieves or initializes an edge in the dictionary.

```julia
existingnodes = get!(() -> newnodes, dict, key)
```

--------------------------------

### Get Pixel Edges

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Extracts edges from the input array.

```julia
edges = _pixel_edges(f, xs, ys, A)
```

--------------------------------

### Area operation implementation

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/types/operation

Example of an Area operation struct and its callable implementation.

```julia
struct Area{Alg<: Algorithm} <: Operation{Alg}
    alg::Alg
end

Area() = Area(NoAlgorithm(Planar()))

function (op::Area{Alg})(data; threaded = _False(), init = 0.0) where {Alg <: Algorithm}
    return GO.area(op.alg, data; threaded, init)
end
```

--------------------------------

### Calculate rectangle angles

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/angles

Example of calculating the angles of a simple rectangle geometry.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie, CairoMakie

rect = GI.Polygon([[(0.0, 0.0), (0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)]])
f, a, p = poly(collect(GI.getpoint(rect)); axis = (; aspect = DataAspect()))
```

```julia
GO.angles(rect)  # [90, 90, 90, 90]
```

--------------------------------

### Create a single Point

Source: https://juliageo.org/GeometryOps.jl/stable/tutorials/creating_geometry

Initialize a basic 2D point geometry.

```julia
point = GI.Point(0, 0)
```

--------------------------------

### Custom operation implementation

Source: https://juliageo.org/GeometryOps.jl/stable/source/src/types/operation

Example of defining a custom operation struct that inherits from Operation.

```julia
struct XPlusOneOperation{M <: Manifold} <: Operation{NoAlgorithm{M}}
    m::M # the manifold always needs to be stored, since it's not a singleton
    x::Int
end
```

--------------------------------

### Check LineString Parallelism

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/orientation

Example usage of the isparallel function to compare two LineStrings.

```julia
julia import GeoInterface as GI, GeometryOps as GO julia> line1 = GI.LineString([(9.170356, 45.477985), (9.164434, 45.482551), (9.166644, 45.484003)]) GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(9.170356, 45.477985), (9.164434, 45.482551), (9.166644, 45.484003)], nothing, nothing)
julia> line2 = GI.LineString([(9.169356, 45.477985), (9.163434, 45.482551), (9.165644, 45.484003)]) GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(9.169356, 45.477985), (9.163434, 45.482551), (9.165644, 45.484003)], nothing, nothing)
julia> GO.isparallel(line1, line2) true
```

--------------------------------

### Validate endpoint distance

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/simplify

Checks if the start or endpoint distance meets the criteria for simplification.

```julia
        pre_pt, post_pt = points[sorted_results[end - 1]], points[sorted_results[2]]
        endpt_dist = _squared_distance_line(Float64, points[1], pre_pt, post_pt)
        if !isnothing(alg.tol)
```

--------------------------------

### Requirement Initialization

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/geom_geom_processors

Initializes requirement flags based on input parameters.

```julia
    in_req_met = !in_require
    on_req_met = !on_require
    out_req_met = !out_require
```

--------------------------------

### Close chain from middle

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/clipping_processor

Closes an overlapping chain if the process started in the middle of the chain.

```julia
    if unmatched_end_chain_edge != unknown
        crossing = unmatched_end_chain_edge != start_chain_edge
```

--------------------------------

### Import Geospatial Libraries

Source: https://juliageo.org/GeometryOps.jl/stable/tutorials/geospatial_geometry

Load the required modules for geometry handling, projections, plotting, and data I/O.

```julia
# Geospatial packages
import GeoInterface as GI
import GeometryOps as GO
import GeoFormatTypes as GFT
# Coordinate transformation and projection
import CoordinateTransformations
import Proj
# Plotting
using CairoMakie
using GeoMakie
# Data loading
using GeoDataFrames, DataFrames
import GeoParquet, GeoJSON, Shapefile # activate native IO packages
```

--------------------------------

### Get union of values

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/polygonize

Resolves unique return values from the function applied to the matrix.

```julia
values = map(identity, sort!(Base.union(Iterators.map(f, A))))
```

--------------------------------

### Simplify Polygon with Point Count

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/simplify

Example of simplifying a polygon by specifying the target number of points.

```julia
using CairoMakie
import GeoInterface as GI
import GeometryOps as GO

original = 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(original; number=6)

f, a, p = poly(original; label = "Original")
poly!(simple; label = "Simplified")
axislegend(a)
f
```

--------------------------------

### Simplify Polygon Example

Source: https://juliageo.org/GeometryOps.jl/stable/source/transformations/simplify

Reduces a polygon to a specific number of points using the simplify function.

```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
```

--------------------------------

### Simplify geometry with GEOS

Source: https://juliageo.org/GeometryOps.jl/stable/source/GeometryOpsLibGEOSExt/simplify

Demonstrates basic usage of the GEOS simplification algorithm with a tolerance parameter.

```julia
import GeometryOps as GO, GeoInterface as GI
import LibGEOS # activate the GEOS algorithm
using GADM # get data
using CairoMakie # to plot

france = GI.geometry(GI.getfeature(GADM.get("FRA"), 1))

simplified = GO.simplify(GO.GEOS(; tol = 1), france)

f, a, p = lines(france; label = "Original")
lines!(a, simplified; label = "Simplified")
axislegend(a)
f
```

--------------------------------

### Include module dependencies

Source: https://juliageo.org/GeometryOps.jl/stable/source/utils/UnitSpherical/UnitSpherical

Loads internal files for point definitions and robust cross product calculations.

```julia
using TestItems # this is a thin package that allows TestItems.@testitem to be parsed.
include("point.jl")

include("robustcrossproduct/RobustCrossProduct.jl")
```

--------------------------------

### MonotoneChainMethod Implementation Steps

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/convex_hull

Internal implementation steps for the MonotoneChainMethod, including point extraction, triangulation, and polygon conversion.

```julia
convex_hull(geometries) = convex_hull(MonotoneChainMethod(), geometries)
```

```julia
function convex_hull(::MonotoneChainMethod, geometries)
```

```julia
    points = collect(flatten(tuples, GI.PointTrait, geometries))
```

```julia
    hull = DelaunayTriangulation.convex_hull(points)
```

```julia
    point_vec = DelaunayTriangulation.get_points(hull)[DelaunayTriangulation.get_vertices(hull)]
    return GI.Polygon([GI.LinearRing(point_vec)])
end
```

--------------------------------

### Visualize Geometry Containment

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/geom_relations/contains

Example using Makie to visualize two LineStrings where one contains the other.

```julia
import GeometryOps as GO
import GeoInterface as GI
using Makie
using CairoMakie

l1 = GI.LineString([(0.0, 0.0), (1.0, 0.0), (0.0, 0.1)])
l2 = GI.LineString([(0.25, 0.0), (0.75, 0.0)])
f, a, p = lines(GI.getpoint(l1), color = :blue)
scatter!(GI.getpoint(l1), color = :blue)
lines!(GI.getpoint(l2), color = :orange)
scatter!(GI.getpoint(l2), color = :orange)
f
```

--------------------------------

### Coverage Implementation Constants and Methods

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/coverage

Core implementation constants and the main coverage dispatch methods.

```julia
const UNKNOWN, NORTH, EAST, SOUTH, WEST = 0:4

"""
    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.
"""
function coverage(geom, xmin, xmax, ymin, ymax,::Type{T} = Float64; threaded=false) where T <: AbstractFloat
    applyreduce(+, _COVERAGE_TARGETS, geom; threaded, init=zero(T)) do g
        _coverage(T, GI.trait(g), g, T(xmin), T(xmax), T(ymin), T(ymax))
    end
end

function coverage(geom, cell_ext::Extents.Extent, ::Type{T} = Float64; threaded=false) where T <: AbstractFloat
    (xmin, xmax), (ymin, ymax) = values(cell_ext)
    return coverage(geom, xmin, xmax, ymin, ymax, T; threaded = threaded)
end
```

--------------------------------

### Update PolyNode Chain Tags

Source: https://juliageo.org/GeometryOps.jl/stable/source/methods/clipping/clipping_processor

Updates the start of a chain with endpoint and crossing or bouncing tags.

```julia
        start_pt = a_list[start_chain_idx]
        a_list[start_chain_idx] = PolyNode(start_pt;
            crossing = crossing,
            endpoint = start_chain,
        )
        b_list[start_pt.neighbor] = PolyNode(b_list[start_pt.neighbor];
            crossing = crossing,
            endpoint = same_winding ? start_chain : end_chain,
        )
    end
end
```
Tokens:
Raw