4  Geometric processing

In this chapter we give an overview of some of the advanced geometric processing features of the framework, highlighting the ones that will be most useful in future chapters.

From now on, we will assume that the main module of the framework is loaded in all code examples:

using GeoStats

We will also assume that the Makie.jl backend is loaded:

import CairoMakie as Mke

4.1 Geometries

We provide a vast list of geometries, which are organized into two main classes, Primitive and Polytope geometries. A geometry is a Primitive if it can be represented efficiently without discretization. For example, we can represent a Box with two corner points or a Ball with center and radius:

box = Box((0, 0, 0), (1, 1, 1))
ball = Ball((0, 0, 2), 0.5)

viz([box, ball], color = ["teal", "slategray3"])

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

Other examples include the Cylinder:

cyl = Cylinder(1.0)

viz(cyl)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

And the Torus:

torus = Torus((1, 1, 0), (-1, -1, 0), (1, -1, 0), 0.5)

viz(torus)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

The full list can be obtained with Julia’s subtypes function:

subtypes(Primitive)

20-element Vector{Any}:
 Ball
 BezierCurve
 Box
 Circle
 Cone
 ConeSurface
 Cylinder
 CylinderSurface
 Disk
 Ellipsoid
 Frustum
 FrustumSurface
 Line
 ParaboloidSurface
 ParametrizedCurve
 Plane
 Point
 Ray
 Sphere
 Torus

Note

Geometries have their own documentation. For example, to learn more about the Cylinder geometry, its parameters and defaults, we can type the following in the Julia REPL:

?Cylinder

A geometry is a Polytope if it can be represented as a combination of “flat sides”, which are also Polytope themselves. A 3D Polytope is called a Polyhedron, a 2D Polytope is called a Polygon and a 1D Polytope is called a polygonal Chain. All these geometries are represented internally with a list of vertices.

First, let’s take a look into the Polyhedron geometries:

subtypes(Polyhedron)

4-element Vector{Any}:
 Hexahedron
 Pyramid
 Tetrahedron
 Wedge

The Hexahedron is a generalization of a 3D Box in the sense that it doesn’t need to be aligned with the coordinate system:

hex = Hexahedron((0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0),
                 (0, 0, 1), (1, 0, 1), (1, 1, 1), (0, 1, 1))

viz(hex)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

In this case, we need to store all the 8 vertices instead of just the corner points. Other examples of Polyhedron include the Tetrahedron and the Pyramid.

Now, let’s move to the Polygon geometries:

subtypes(Polygon)

2-element Vector{Any}:
 Ngon
 PolyArea

We provide two types of Polygon that meet different application requirements.

The Ngon is a polygon without holes. Its vertices are stored in static memory, and they are mostly used for discretization of other geometries and geospatial domains. We provide type aliases to construct Ngon with a specific number N of vertices:

Triangle, Quadrangle, Pentagon, …, Decagon

The Quadrangle is a generalization of the 2D Box in the sense that it doesn’t need to be aligned with the coordinate system:

t = Triangle((0, 0), (1, 0), (1, 1))
q = Quadrangle((1, 1), (2, 1), (2, 2), (1, 2))

viz([t, q], color = ["teal", "slategray3"])

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

The PolyArea is a polygon with or without holes. Its vertices are stored in dynamic memory, and they are mostly used for representation of polygonal areas in GIS:

outer = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
hole1 = [(0.2, 0.2), (0.2, 0.4), (0.4, 0.4), (0.4, 0.2)]
hole2 = [(0.6, 0.2), (0.6, 0.4), (0.8, 0.4), (0.8, 0.2)]
poly  = PolyArea([outer, hole1, hole2])

viz(poly)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

In the example above, the first list of vertices represents the external boundary of the PolyArea, also known as the outer Ring. The other two lists represent the two internal boundaries, or inner rings. A single list of vertices can be used, in which case the PolyArea doesn’t have holes.

Note

The topology of a 2D Polygon is defined by the orientation of its rings. Most geometric processing algorithms expect counter clockwise (CCW) orientation for the outer ring, and clockwise (CW) orientation for the inner rings:

r = rings(poly)

orientation.(r)

3-element Vector{OrientationType}:
 CCW::OrientationType = 1
 CW::OrientationType = 0
 CW::OrientationType = 0

Tip for all users

The Repair transform covered in the next chapter can be used to repair the orientations of rings in polygons.

Finally, let’s take a look into the polygonal Chain:

subtypes(Chain)

3-element Vector{Any}:
 Ring
 Rope
 Segment

These are 1-dimensional polytopes connecting Points in sequence. We have seen the Rings in the PolyArea and Ngon geometries:

viz(r)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

A Ring is closed meaning that its first and last Points are connected with a Segment. A Rope is an open Ring without the closing Segment:

viz(open.(r))

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

We can obtain the list of segments of a Chain with the segments function:

collect(segments(first(r)))

4-element Vector{Segment{𝔼{2}, Cartesian2D{NoDatum, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Segment((x: 0.0 m, y: 0.0 m), (x: 1.0 m, y: 0.0 m))
 Segment((x: 1.0 m, y: 0.0 m), (x: 1.0 m, y: 1.0 m))
 Segment((x: 1.0 m, y: 1.0 m), (x: 0.0 m, y: 1.0 m))
 Segment((x: 0.0 m, y: 1.0 m), (x: 0.0 m, y: 0.0 m))

The Segment geometry is a Chain with just 2 vertices:

viz(Segment((0, 0), (1, 1)))

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

Finally, there is the Multi-geometry, which is a set of geometries seen as a single geometry. This is very common in GIS to represent disconnected areas on a geographic map that are related to each other (e.g., political areas):

Multi([Point(1, 2), Point(2, 3)])

MultiPoint
├─ Point(x: 1.0 m, y: 2.0 m)
└─ Point(x: 2.0 m, y: 3.0 m)

Multi(r)

MultiRing
├─ Ring((x: 0.0 m, y: 0.0 m), ..., (x: 0.0 m, y: 1.0 m))
├─ Ring((x: 0.2 m, y: 0.2 m), ..., (x: 0.4 m, y: 0.2 m))
└─ Ring((x: 0.6 m, y: 0.2 m), ..., (x: 0.8 m, y: 0.2 m))

Multi([t, q])

MultiNgon
├─ Triangle((x: 0.0 m, y: 0.0 m), (x: 1.0 m, y: 0.0 m), (x: 1.0 m, y: 1.0 m))
└─ Quadrangle((x: 1.0 m, y: 1.0 m), ..., (x: 1.0 m, y: 2.0 m))

4.2 Predicates

Julia provides support for unicode characters in variable and function names. We leverage this feature to define commonly used geometric predicates with intuitive mathematical notation:

p = Point(0.0, 0.0)
b = Ball((0.5, 0.5), 1.0)

viz([p, b], color = ["teal", "slategray3"])

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

p ∈ b

true

b1 = Box((0, 0), (1, 1))
b2 = Box((0.5, 0.5), (2, 2))

viz([b1, b2], color = ["teal", "slategray3"])

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

b1 ⊆ b2

false

Tip for all users

The symbol ∈ is obtained in Julia by typing \in and pressing the TAB key on the keyboard. We could have used the syntax p in b or in(p, b) as well. Similarly, the symbol ⊆ is obtained by typing \subseteq. We could have used the syntax issubseteq(b1, b2) as well.

If you don’t know the \(\LaTeX\) name of a symbol, you can copy/paste it in the Julia REPL in help mode:

?∈

Some predicates don’t have well-established mathematical notation. For example, a polygon issimple if it doesn’t have holes nor self-intersections:

q = Quadrangle((0, 0), (1, 0), (1, 1), (0.6, 0.4))

viz(q)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

issimple(q)

true

It isconvex if all line segments connecting two points of the polygon are inside the polygon:

isconvex(q)

false

A very useful predicate is intersects (with a “s” at the end):

outer = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
hole1 = [(0.2, 0.2), (0.2, 0.4), (0.4, 0.4), (0.4, 0.2)]
hole2 = [(0.6, 0.2), (0.6, 0.4), (0.8, 0.4), (0.8, 0.2)]
poly  = PolyArea([outer, hole1, hole2])
ball1 = Ball((0.5, 0.5), 0.05)
ball2 = Ball((0.3, 0.3), 0.05)

viz([poly, ball1, ball2], color = ["slategray3", "teal", "brown"])

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

intersects(poly, ball1)

true

intersects(poly, ball2)

false

It tells whether or not the geometries intersect, without actually computing the intersection. The intersection itself is considered a geometric operation as discussed in the next section.

Please consult the official documentation for the full list of predicates.

4.3 Operations

Geometric operations transform a geometry or a set of geometries into a new geometry or number. For example, the intersection of two segments can be a Point, a Segment or nothing:

s1 = Segment((0.0, 0.0), (1.0, 0.0))
s2 = Segment((0.5, 0.0), (2.0, 0.0))

s1 ∩ s2

Segment
├─ Point(x: 0.5 m, y: 0.0 m)
└─ Point(x: 1.0 m, y: 0.0 m)

Tip for advanced users

For performance-sensitive applications, it is wise to replace the ∩ operation by its the 3-argument version named intersection:

intersection(s1, s2) do I
  if I == Crossing
    return 1
  else
    return 0
  end
end

0

The example above uses Julia’s do-syntax to define a function in place. The function takes the intersection type I and creates branches that return the same type (Int in this case) for type stability. The more we reduce the number of branches and types, the more the Julia compiler will be able to infer the output type.

Likewise, the intersection of two 2D geometries can be obtained with:

poly1 = Ngon((0, 1), (1, 0), (2, 1), (3, 0), (4, 1), (2, 2))
poly2 = Ngon((1.0, 0.5), (3.5, 0.0), (3.5, 1.5), (1.0, 1.5))

poly  = poly1 ∩ poly2

viz(poly1)
viz!(poly2, color = "teal", alpha = 0.2)
viz!(boundary(poly), color = "red")
Mke.current_figure()

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

The previous example makes use of the boundary of a geometry, which is very useful to know:

boundary(poly)

Ring
├─ Point(x: 1.0 m, y: 0.5 m)
├─ Point(x: 1.4166666666666667 m, y: 0.4166666666666667 m)
├─ Point(x: 2.0 m, y: 1.0 m)
├─ Point(x: 2.875 m, y: 0.125 m)
├─ Point(x: 3.0833333333333335 m, y: 0.08333333333333333 m)
├─ Point(x: 3.5 m, y: 0.4999999999999999 m)
├─ Point(x: 3.5 m, y: 1.25 m)
├─ Point(x: 3.0 m, y: 1.5 m)
└─ Point(x: 1.0 m, y: 1.5 m)

Some operations like measure (length, area or volume) produce numbers instead of geometries. For example, the area of the polygon above is:

area(poly)

2.4479166666666665 m^2

The measure of the boundary is known as the perimeter of the geometry:

perimeter(poly)

7.598044863023599 m

All Polytope geometries have vertices:

vertices(poly)

9-element CircularVector(::Vector{Point{𝔼{2}, Cartesian2D{NoDatum, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}):
 Point(x: 1.0 m, y: 0.5 m)
 Point(x: 1.4166666666666667 m, y: 0.4166666666666667 m)
 Point(x: 2.0 m, y: 1.0 m)
 Point(x: 2.875 m, y: 0.125 m)
 Point(x: 3.0833333333333335 m, y: 0.08333333333333333 m)
 Point(x: 3.5 m, y: 0.4999999999999999 m)
 Point(x: 3.5 m, y: 1.25 m)
 Point(x: 3.0 m, y: 1.5 m)
 Point(x: 1.0 m, y: 1.5 m)

Please consult the official documentation for the full list of operations.

4.4 Algorithms

Any other function that is not a predicate nor an operation is called a geometric processing “algorithm” in the framework. We provide a list of advanced algorithms for discretization, simplification, refinement, convex hull, etc.

Below we illustrate some of these algorithms, which will be useful in future examples:

points = rand(Point, 100, crs=Cartesian2D)

100-element Vector{Point{𝔼{2}, Cartesian2D{NoDatum, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Point(x: 0.4041032663827714 m, y: 0.21106135658332847 m)
 Point(x: 0.1735116868658847 m, y: 0.9218222262047227 m)
 Point(x: 0.9797364421806157 m, y: 0.422814088993162 m)
 Point(x: 0.7390860525281144 m, y: 0.01961430765102179 m)
 Point(x: 0.3162711274531963 m, y: 0.24051326979982934 m)
 Point(x: 0.18563656353941305 m, y: 0.15760551003859402 m)
 Point(x: 0.44086393844800076 m, y: 0.7765008215821155 m)
 Point(x: 0.6315769049861976 m, y: 0.901629086873425 m)
 Point(x: 0.8452483046703105 m, y: 0.628306934661662 m)
 Point(x: 0.03553405934028142 m, y: 0.011075950306849824 m)
 Point(x: 0.28988561447025285 m, y: 0.2983872850203414 m)
 Point(x: 0.48712634692439016 m, y: 0.7351152256092213 m)
 Point(x: 0.6746648063891224 m, y: 0.7161685799246165 m)
 ⋮
 Point(x: 0.5968013398134062 m, y: 0.6805836367494245 m)
 Point(x: 0.31900298671619 m, y: 0.44974721394677886 m)
 Point(x: 0.6202387301808445 m, y: 0.22792559463924855 m)
 Point(x: 0.666437184604755 m, y: 0.24073407463221164 m)
 Point(x: 0.43133020978647574 m, y: 0.918094017984448 m)
 Point(x: 0.6633200676542864 m, y: 0.4668826538041201 m)
 Point(x: 0.11161870696420417 m, y: 0.6977709278333655 m)
 Point(x: 0.7814870232550928 m, y: 0.9382715827458015 m)
 Point(x: 0.48902494600650814 m, y: 0.6248620120538431 m)
 Point(x: 0.4150572459510692 m, y: 0.4930178053067755 m)
 Point(x: 0.39419736463609467 m, y: 0.5623553170314027 m)
 Point(x: 0.9227205477159082 m, y: 0.9491488530252664 m)

Note

The crs option specifies the Coordinate Reference System (CRS) of the points. It will be explored in the chapter Map projections. By default, the rand function generates geometries with Cartesian3D CRS:

rand(Triangle, 3)

3-element Vector{Triangle{𝔼{3}, Cartesian3D{NoDatum, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Triangle((x: 0.755918 m, y: 0.174028 m, z: 0.387394 m), (x: 0.646953 m, y: 0.225008 m, z: 0.026513 m), (x: 0.834323 m, y: 0.995491 m, z: 0.254951 m))
 Triangle((x: 0.057521 m, y: 0.991609 m, z: 0.639189 m), (x: 0.0777913 m, y: 0.98413 m, z: 0.0358643 m), (x: 0.751387 m, y: 0.984089 m, z: 0.286951 m))
 Triangle((x: 0.761334 m, y: 0.934571 m, z: 0.165763 m), (x: 0.74681 m, y: 0.409709 m, z: 0.200255 m), (x: 0.782853 m, y: 0.86964 m, z: 0.20033 m))

viz(boundingbox(points))
viz!(points, color = "black")
Mke.current_figure()

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

viz(convexhull(points))
viz!(points, color = "black")
Mke.current_figure()

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

Geometries can be discretized into geospatial domains (i.e., collections of geometries). We have seen some of these domains in previous chapters, including the GeometrySet and the CartesianGrid. Here, we will focus on the Mesh subtypes:

subtypes(Mesh)

5-element Vector{Any}:
 RectilinearGrid
 RegularGrid
 SimpleMesh
 StructuredGrid
 TransformedMesh

subtypes(Grid)

5-element Vector{Any}:
 RectilinearGrid
 RegularGrid
 SimpleMesh{M, C, V, GridTopology{Dim}} where {M<:Meshes.Manifold, C<:CRS, V<:AbstractArray{Point{M, C}, 1}, Dim}
 StructuredGrid
 TransformedMesh{M, C, GridTopology{Dim}} where {M<:Meshes.Manifold, C<:CRS, Dim}

There are two main functions to perform discretization: discretize and simplexify. The discretize function takes a geometry and a discretization algorithm as input, and produces a Mesh. For example, we can discretize a 2D ball with regular spacing along each parametric dimension (i.e., polar coordinates):

ball = Ball((0, 0), 1)

mesh = discretize(ball, RegularDiscretization(20, 50))

1050 SimpleMesh
  1051 vertices
  ├─ Point(x: 0.047619047619047616 m, y: 0.0 m)
  ├─ Point(x: 0.09523809523809523 m, y: 0.0 m)
  ├─ Point(x: 0.14285714285714285 m, y: 0.0 m)
  ├─ Point(x: 0.19047619047619047 m, y: 0.0 m)
  ├─ Point(x: 0.23809523809523808 m, y: 0.0 m)
  ⋮
  ├─ Point(x: 0.8503840296981238 m, y: -0.10742848591226112 m)
  ├─ Point(x: 0.8976275869035751 m, y: -0.11339673512960896 m)
  ├─ Point(x: 0.9448711441090264 m, y: -0.1193649843469568 m)
  ├─ Point(x: 0.9921147013144778 m, y: -0.12533323356430465 m)
  └─ Point(x: 0.0 m, y: 0.0 m)
  1050 elements
  ├─ Quadrangle(1, 2, 23, 22)
  ├─ Quadrangle(2, 3, 24, 23)
  ├─ Quadrangle(3, 4, 25, 24)
  ├─ Quadrangle(4, 5, 26, 25)
  ├─ Quadrangle(5, 6, 27, 26)
  ⋮
  ├─ Triangle(1051, 946, 967)
  ├─ Triangle(1051, 967, 988)
  ├─ Triangle(1051, 988, 1009)
  ├─ Triangle(1051, 1009, 1030)
  └─ Triangle(1051, 1030, 1)

We can visualize the elements of the mesh with the showsegments option:

viz(mesh, showsegments = true)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

The meshes produced by the discretize function can have mixed element types. In this case, the mesh is made of Quadrangle and Triangle geometries. The function simplexify is preferred if the application requires a mesh of simplices (i.e., triangles or tetrahedra):

tmesh = simplexify(ball)

5050 SimpleMesh
  2551 vertices
  ├─ Point(x: 0.0196078431372549 m, y: 0.0 m)
  ├─ Point(x: 0.0392156862745098 m, y: 0.0 m)
  ├─ Point(x: 0.058823529411764705 m, y: 0.0 m)
  ├─ Point(x: 0.0784313725490196 m, y: 0.0 m)
  ├─ Point(x: 0.09803921568627451 m, y: 0.0 m)
  ⋮
  ├─ Point(x: 0.9337550130018614 m, y: -0.1179606904134632 m)
  ├─ Point(x: 0.9532082424394003 m, y: -0.12041820479707702 m)
  ├─ Point(x: 0.972661471876939 m, y: -0.12287571918069083 m)
  ├─ Point(x: 0.9921147013144778 m, y: -0.12533323356430465 m)
  └─ Point(x: 0.0 m, y: 0.0 m)
  5050 elements
  ├─ Triangle(53, 52, 1)
  ├─ Triangle(1, 2, 53)
  ├─ Triangle(54, 53, 2)
  ├─ Triangle(2, 3, 54)
  ├─ Triangle(55, 54, 3)
  ⋮
  ├─ Triangle(2551, 2296, 2347)
  ├─ Triangle(2551, 2347, 2398)
  ├─ Triangle(2551, 2398, 2449)
  ├─ Triangle(2551, 2449, 2500)
  └─ Triangle(2551, 2500, 1)

viz(tmesh, showsegments = true)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

Various triangulation algorithms are provided such as DehnTriangulation (Devadoss and O’Rourke 2011), DelaunayTriangulation (Cheng, Dey, and Shewchuk 2012) and HeldTriangulation (Held 2001). After the mesh is created, it is possible to refine the elements:

mesh = discretize(Box((0, 0), (1, 1)))

2 SimpleMesh
  4 vertices
  ├─ Point(x: 0.0 m, y: 0.0 m)
  ├─ Point(x: 1.0 m, y: 0.0 m)
  ├─ Point(x: 1.0 m, y: 1.0 m)
  └─ Point(x: 0.0 m, y: 1.0 m)
  2 elements
  ├─ Triangle(1, 2, 3)
  └─ Triangle(1, 3, 4)

tmesh = refine(mesh, TriRefinement())

viz(tmesh, showsegments = true)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

qmesh = refine(mesh, QuadRefinement())

viz(qmesh, showsegments = true)

┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

We will have the chance to see more algorithms in action as we advance in the chapters of the book.

4.5 Domains

For most purposes, domains can be manipulated as if they were Julia vectors of geometries. As an example, we can compute the centroid of each geometry in the domain by broadcasting the function:

grid = CartesianGrid(2, 3)

centroid.(grid)

6-element Vector{Point{𝔼{2}, Cartesian2D{NoDatum, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}:
 Point(x: 0.5 m, y: 0.5 m)
 Point(x: 1.5 m, y: 0.5 m)
 Point(x: 0.5 m, y: 1.5 m)
 Point(x: 1.5 m, y: 1.5 m)
 Point(x: 0.5 m, y: 2.5 m)
 Point(x: 1.5 m, y: 2.5 m)

In some cases, however, it is important to know which geometries are adjacent to a given geometry; or which geometries make up the boundary of a given geometry. This is where the topology is useful:

topo = topology(grid)

2×3 GridTopology(aperiodic, aperiodic)

We can create an Adjacency topological relation to find which quadrangles are adjacent to quadrangle 1 in the domain:

A = Adjacency{2}(topo)

A(1)

(2, 3)

The number 2 that appears in the Adjacency relation is known as the parametric dimension. In this case, we are interested in the adjacency of 2D geometries. The quadrangle 1 is the first quadrangle in the bottom-left corner of the grid, and it is adjacent to quadrangles 2 and 3. As another example, the quadrangle 3 is adjacent to quadrangles 4, 1 and 5:

A(3)

(4, 1, 5)

In order to query the adjacency of vertices, we create a relation with parametric dimension 0:

A = Adjacency{0}(topo)

A(1)

(2, 4)

We can also query the vertices that are on the boundary of a given quadrangle with the Boundary topological relation. In this case, we specify the parametric dimension of the input and output geometries:

B = Boundary{2,0}(topo)

B(1)

(1, 2, 5, 4)

Topological relations are an advanced feature of the framework. They are mostly used by developers of geostatistical algorithms, or in geospatial queries that will be covered in Part III of the book.

4.6 Congratulations!

Congratulations on finishing Part I of the book. Let’s quickly review what we learned so far:

• The main concept in geospatial data science is the concept of geospatial data, represented in the GeoStats.jl framework as geotables over geospatial domains.
• The geotable representation generalizes traditional GIS representations (“raster” vs. “vector”), and enables an unified approach to visualization and manipulation of geospatial data.
• It is still possible to interface with existing GIS technology via input and output of files using the GeoIO.jl module. A typical workflow will load GIS data at the beginning of a script, and save the result of the analysis at the end of the script.
• Geometric processing doesn’t need to be complicated. It should be fun and read like math. If it feels “computer sciency”, that is a limitation of the software and programming language.

We are finally ready to learn the advanced features of the framework. Let’s get it started.

Cheng, Siu-Wing, Tamal K. Dey, and Jonathan Shewchuk. 2012. Delaunay Mesh Generation. Chapman; Hall/CRC. https://people.eecs.berkeley.edu/~jrs/meshbook.html.

Devadoss, Satyan L, and Joseph O’Rourke. 2011. Discrete and Computational Geometry. Princeton, NJ: Princeton University Press.

Held, M. 2001. “FIST: Fast Industrial-Strength Triangulation of Polygons.” Algorithmica 30 (4): 563–96. https://doi.org/10.1007/s00453-001-0028-4.