Skip to content

Commit 30bf6b0

Browse files
asinghvi17asinghvi17
authored
Merge branch 'main' into as/clipping_crs
2 parents 881b221 + a7620d6 commit 30bf6b0

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

41 files changed

+3193
-196
lines changed

‎Project.toml‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "GeometryOps"
22
uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
33
authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
4-
version = "0.1.26"
4+
version = "0.1.30"
55

66
[deps]
77
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -15,6 +15,7 @@ GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
1515
GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
1616
GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013"
1717
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
18+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
1819
SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
1920
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
2021
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -51,7 +52,8 @@ GeometryOpsCore = "=0.1.7"
5152
LibGEOS = "0.9.2"
5253
LinearAlgebra = "1"
5354
Proj = "1"
54-
SortTileRecursiveTree = "0.1"
55+
Random = "1"
56+
SortTileRecursiveTree = "0.1.2"
5557
StaticArrays = "1"
5658
Statistics = "1"
5759
TGGeometry = "0.1"

‎docs/Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,6 @@ GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
2323
GeoDatasets = "ddc7317b-88db-5cb5-a849-8449e5df04f9"
2424
GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
2525
GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
26-
GeoInterfaceMakie = "0edc0954-3250-4c18-859d-ec71c1660c08"
2726
GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
2827
GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
2928
GeoParquet = "e99870d8-ce00-4fdd-aeee-e09192881159"
@@ -53,3 +52,4 @@ GeometryOpsCore = {path = "../GeometryOpsCore"}
5352

5453
[compat]
5554
DocumenterVitepress = ">=0.2.1"
55+
GeoInterface = "1.6"

‎docs/make.jl‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,7 @@ CairoMakie.activate!(px_per_unit = 2, type = "svg", inline = true) # TODO: make
77
# import packages that activate extensions
88
import FlexiJoins, LibGEOS, Proj, TGGeometry
99

10-
DocMeta.setdocmeta!(GeometryOps, :DocTestSetup, :(using GeometryOps; using GeometryBasics); recursive=true)
11-
12-
using GeoInterfaceMakie
10+
DocMeta.setdocmeta!(GeometryOps, :DocTestSetup, :(using GeometryOps; using GeometryBasics; using GeometryOps.GeometryOpsCore); recursive=true)
1311

1412
include(joinpath(dirname(@__DIR__), "benchmarks", "benchmark_plots.jl"))
1513

‎docs/src/experiments/spatial_index_building.jl‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
using CairoMakie
2+
import GeoInterface as GI, GeometryOps as GO
3+
using NaturalEarth
4+
5+
6+
using GeometryOps.SpatialTreeInterface
7+
import GeometryOps.SpatialTreeInterface as STI
8+
using SortTileRecursiveTree
9+
10+
function build_spatial_index_gif(geom, index_constructor, filename; plot_leaves = true, title = splitext(filename)[1], axis = (;), figure = (;), record = (;))
11+
fig = Figure(; figure...)
12+
ax = Axis(fig[1, 1]; title = title, axis...)
13+
14+
# Create a spatial index
15+
index = index_constructor(geom)
16+
ext = STI.node_extent(index)
17+
limits!(ax, ext.X[1], ext.X[2], ext.Y[1], ext.Y[2])
18+
19+
rects = Rect2f[Rect2f((NaN, NaN), (NaN, NaN))]
20+
colors = RGBAf[to_color(:transparent)]
21+
palette = Makie.wong_colors(0.7)
22+
23+
plt = poly!(ax, rects; color = colors)
24+
25+
to_rect2(extent) = Rect2f((extent.X[1], extent.Y[1]), (extent.X[2] - extent.X[1], extent.Y[2] - extent.Y[1]))
26+
27+
function dive_in(io, plt, node, level)
28+
if STI.isleaf(node) && plot_leaves
29+
push!(rects, to_rect2(STI.node_extent(node)))
30+
push!(colors, palette[level])
31+
else
32+
for child in STI.getchild(node)
33+
dive_in(io, plt, child, level + 1)
34+
end
35+
push!(rects, to_rect2(STI.node_extent(node)))
36+
push!(colors, palette[level])
37+
end
38+
update!(plt, rects; color = colors)
39+
recordframe!(io)
40+
return
41+
end
42+
43+
Makie.record(fig, filename; record...) do io
44+
empty!(rects)
45+
empty!(colors)
46+
dive_in(io, plt, index, 1)
47+
end
48+
49+
end

‎ext/GeometryOpsProjExt/GeometryOpsProjExt.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,5 +7,6 @@ using GeometryOps, Proj
77

88
include("reproject.jl")
99
include("segmentize.jl")
10+
include("perimeter_area.jl")
1011

1112
end

‎ext/GeometryOpsProjExt/perimeter_area.jl‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
import GeometryOps: perimeter, area, applyreduce, TraitTarget, WithTrait
2+
import GeoInterface as GI
3+
4+
function perimeter(m::Geodesic, geom, ::Type{T} = Float64; init = zero(T), kwargs...) where T
5+
# Create a Proj geodesic object using the ellipsoid parameters from the Geodesic manifold
6+
proj_geodesic = Ref(Proj.geod_geodesic(m.semimajor_axis, 1/m.inv_flattening))
7+
proj_polygon = Ref(Proj._null(Proj.geod_polygon))
8+
9+
function _perimeter_geodesic_inner(trait, geom)
10+
@assert GI.npoint(geom) >= 2 "Geodesic perimeter requires at least 2 points"
11+
12+
# Initialize the polygon
13+
proj_polygon[] = Proj._null(Proj.geod_polygon)
14+
Proj.geod_polygon_init(proj_polygon, 1)
15+
16+
# Add all points to the polygon
17+
for point in GI.getpoint(trait, geom)
18+
lat, lon = GI.y(point), GI.x(point) # Proj expects lat, lon order
19+
Proj.geod_polygon_addpoint(proj_geodesic, proj_polygon, lat, lon)
20+
end
21+
22+
# Compute the polygon properties
23+
# geod_polygon_compute returns (num_vertices, perimeter, area)
24+
area_result, perimeter_result = Proj.geod_polygon_compute(proj_geodesic[], proj_polygon[], false, true)
25+
26+
return T(perimeter_result)
27+
end
28+
29+
return applyreduce(
30+
WithTrait(_perimeter_geodesic_inner),
31+
+,
32+
TraitTarget(GI.AbstractCurveTrait),
33+
geom; init, kwargs...
34+
)
35+
end
36+
37+
function GeometryOps.area(m::Geodesic, geom, ::Type{T} = Float64; threaded = False(), init = zero(T), kwargs...) where T
38+
# Create a Proj geodesic object using the ellipsoid parameters from the Geodesic manifold
39+
proj_geodesic = Ref(Proj.geod_geodesic(m.semimajor_axis, 1/m.inv_flattening))
40+
proj_polygon = Ref(Proj._null(Proj.geod_polygon))
41+
42+
function _area_geodesic_inner(trait, geom)
43+
@assert GI.npoint(geom) >= 2 "Geodesic perimeter requires at least 2 points"
44+
45+
# Initialize the polygon
46+
proj_polygon[] = Proj._null(Proj.geod_polygon)
47+
Proj.geod_polygon_init(proj_polygon, 0)
48+
49+
# Add all points to the polygon
50+
for point in GI.getpoint(trait, geom)
51+
lat, lon = GI.y(point), GI.x(point) # Proj expects lat, lon order
52+
Proj.geod_polygon_addpoint(proj_geodesic, proj_polygon, lat, lon)
53+
end
54+
55+
# Compute the polygon properties
56+
# geod_polygon_compute returns (num_vertices, perimeter, area)
57+
area_result, perimeter_result = Proj.geod_polygon_compute(proj_geodesic[], proj_polygon[], true, false)
58+
59+
return T(area_result)
60+
end
61+
62+
return applyreduce(
63+
WithTrait(_area_geodesic_inner),
64+
+,
65+
TraitTarget(GI.AbstractCurveTrait),
66+
geom; init, kwargs...
67+
)
68+
end

‎src/GeometryOps.jl‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,12 +21,16 @@ export TraitTarget, Manifold, Planar, Spherical, Geodesic, apply, applyreduce, f
2121
using GeoInterface
2222
using LinearAlgebra, Statistics
2323

24+
using StaticArrays
25+
2426
import Tables, DataAPI
2527
import StaticArrays
2628
import DelaunayTriangulation # for convex hull and triangulation
2729
import ExactPredicates
2830
import Base.@kwdef
2931
import GeoInterface.Extents: Extents
32+
import SortTileRecursiveTree
33+
import SortTileRecursiveTree: STRtree
3034

3135
const GI = GeoInterface
3236

@@ -40,9 +44,17 @@ include("not_implemented_yet.jl")
4044
include("utils/utils.jl")
4145
include("utils/LoopStateMachine/LoopStateMachine.jl")
4246
include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
47+
include("utils/UnitSpherical/UnitSpherical.jl")
48+
49+
# Load utility modules in
50+
using .LoopStateMachine, .SpatialTreeInterface, .UnitSpherical
51+
52+
include("utils/NaturalIndexing.jl")
53+
using .NaturalIndexing
4354

44-
using .LoopStateMachine, .SpatialTreeInterface
4555

56+
# Load utility modules in
57+
using .NaturalIndexing, .SpatialTreeInterface, .LoopStateMachine
4658

4759
include("methods/angles.jl")
4860
include("methods/area.jl")
@@ -52,6 +64,7 @@ include("methods/centroid.jl")
5264
include("methods/convex_hull.jl")
5365
include("methods/distance.jl")
5466
include("methods/equals.jl")
67+
include("methods/perimeter.jl")
5568
include("methods/clipping/predicates.jl")
5669
include("methods/clipping/clipping_processor.jl")
5770
include("methods/clipping/coverage.jl")
@@ -77,6 +90,7 @@ include("transformations/flip.jl")
7790
include("transformations/reproject.jl")
7891
include("transformations/segmentize.jl")
7992
include("transformations/simplify.jl")
93+
include("transformations/smooth.jl")
8094
include("transformations/tuples.jl")
8195
include("transformations/transform.jl")
8296
include("transformations/forcedims.jl")

‎src/methods/area.jl‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -66,8 +66,12 @@ This is computed slightly differently for different geometries:
6666
Result will be of type T, where T is an optional argument with a default value
6767
of Float64.
6868
"""
69-
function area(geom, ::Type{T} = Float64; threaded=false) where T <: AbstractFloat
70-
applyreduce(WithTrait((trait, g) -> _area(T, trait, g)), +, _AREA_TARGETS, geom; threaded, init=zero(T))
69+
function area(geom, ::Type{T} = Float64; threaded=false, kwargs...) where T <: AbstractFloat
70+
area(Planar(), geom, T; threaded, kwargs...)
71+
end
72+
73+
function area(::Planar, geom, ::Type{T} = Float64; threaded=false, kwargs...) where T <: AbstractFloat
74+
applyreduce(WithTrait((trait, g) -> _area(T, trait, g)), +, _AREA_TARGETS, geom; threaded, init=zero(T), kwargs...)
7175
end
7276

7377
"""
@@ -151,3 +155,6 @@ function _signed_area(::Type{T}, geom) where T
151155
area += _area_component(p1, p2)
152156
return T(area / 2)
153157
end
158+
159+
160+
# ## Spherical

‎src/methods/barycentric.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ export MeanValue
2727
#=
2828
```@example barycentric
2929
import GeometryOps as GO, GeoInterface as GI
30-
using CairoMakie, GeoInterfaceMakie # plotting
30+
using CairoMakie # plotting
3131
# Define a polygon
3232
polygon_points = [
3333
(0.03, 0.05, 0.00), (0.07, 0.04, 0.02), (0.10, 0.04, 0.04),

0 commit comments

Comments
 (0)