04-geoproc.qmd

# Geometric processing

```{julia}
#| echo: false
#| output: false
import Pkg
Pkg.activate(".")
```

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:

```{julia}
using GeoStats
```

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

```{julia}
import CairoMakie as Mke
```

## 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:

```{julia}
box = Box((0, 0, 0), (1, 1, 1))
ball = Ball((0, 0, 2), 0.5)

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

Other examples include the `Cylinder`:

```{julia}
cyl = Cylinder(1.0)

viz(cyl)
```

And the `Torus`:

```{julia}
torus = Torus((1, 1, 0), (-1, -1, 0), (1, -1, 0), 0.5)

viz(torus)
```

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

```{julia}
subtypes(Primitive)
```

::: {.callout-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:

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

```{julia}
subtypes(Polyhedron)
```

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

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

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:

```{julia}
subtypes(Polygon)
```

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:

```{julia}
t = Triangle((0, 0), (1, 0), (1, 1))
q = Quadrangle((1, 1), (2, 1), (2, 2), (1, 2))

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

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:

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

viz(poly)
```

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.

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

```{julia}
subtypes(Chain)
```

These are 1-dimensional polytopes connecting `Point`s in sequence. We've
seen the `Ring`s in the `PolyArea` and `Ngon` geometries:

```{julia}
r = rings(poly)

viz(r)
```

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

```{julia}
viz(open.(r))
```

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

```{julia}
collect(segments(first(r)))
```

The `Segment` geometry is a `Chain` with just 2 vertices:

```{julia}
viz(Segment((0, 0), (1, 1)))
```

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):

```{julia}
Multi([Point(1, 2), Point(2, 3)])
```

```{julia}
Multi(r)
```

```{julia}
Multi([t, q])
```

## 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:

```{julia}
p = Point(0.0, 0.0)
b = Ball((0.5, 0.5), 1.0)

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

```{julia}
p ∈ b
```

```{julia}
b1 = Box((0, 0), (1, 1))
b2 = Box((0.5, 0.5), (2, 2))

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

```{julia}
b1 ⊆ b2
```

::: {.callout-note}

## 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:

```julia
?∈
```

:::

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

```{julia}
q = Quadrangle((0, 0), (1, 0), (1, 1), (0.6, 0.4))

viz(q)
```

```{julia}
issimple(q)
```

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

```{julia}
isconvex(q)
```

A very useful predicate is `intersects` (with a "s" at the end):

```{julia}
outer = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
hole1 = [(0.2, 0.2), (0.4, 0.2), (0.4, 0.4), (0.2, 0.4)]
hole2 = [(0.6, 0.2), (0.8, 0.2), (0.8, 0.4), (0.6, 0.4)]
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"])
```

```{julia}
intersects(poly, ball1)
```

```{julia}
intersects(poly, ball2)
```

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.

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

```{julia}
s1 = Segment((0.0, 0.0), (1.0, 0.0))
s2 = Segment((0.5, 0.0), (2.0, 0.0))

s1 ∩ s2
```

::: {.callout-note}

## Tip for advanced users

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

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

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:

```{julia}
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()
```

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

```{julia}
boundary(poly)
```

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

```{julia}
area(poly)
```

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

```{julia}
perimeter(poly)
```

All `Polytope` geometries have `vertices`:

```{julia}
vertices(poly)
```

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

## 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:

```{julia}
points = rand(Point{2}, 100)
```

```{julia}
viz(boundingbox(points))
viz!(points, color = "black")
Mke.current_figure()
```

```{julia}
viz(convexhull(points))
viz!(points, color = "black")
Mke.current_figure()
```

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:

```{julia}
subtypes(Mesh)
```

```{julia}
subtypes(Grid)
```

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):

```{julia}
ball = Ball((0, 0), 1)

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

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

```{julia}
viz(mesh, showsegments = true)
```

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):

```{julia}
tmesh = simplexify(ball)
```

```{julia}
viz(tmesh, showsegments = true)
```

Various triangulation algorithms are provided such as `Dehn1899` [@Devadoss2011],
and `FIST` [@Held2001]. After the mesh is created, it is possible to `refine` the
elements:

```{julia}
mesh = discretize(Box((0, 0), (1, 1)))
```

```{julia}
tmesh = refine(mesh, TriRefinement())

viz(tmesh, showsegments = true)
```

```{julia}
qmesh = refine(mesh, QuadRefinement())

viz(qmesh, showsegments = true)
```

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

## 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:

```{julia}
grid = CartesianGrid(2, 3)

centroid.(grid)
```

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:

```{julia}
topo = topology(grid)
```

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

```{julia}
A = Adjacency{2}(topo)

A(1)
```

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

```{julia}
A(3)
```

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

```{julia}
A = Adjacency{0}(topo)

A(1)
```

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:

```{julia}
B = Boundary{2,0}(topo)

B(1)
```

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.

## 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.