02-geoviz.qmd

# Scientific visualization

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

The visualization ecosystem in Julia is evolving very quickly. Among the various
visualization projects, [Makie.jl](https://github.com/MakieOrg/Makie.jl) by
@Danisch2021 is the most advanced for
[scientific visualization](https://en.wikipedia.org/wiki/Scientific_visualization).

Makie.jl is currently organized in backend modules:

- **GLMakie.jl** is the preferred backend for interactive *high-performance* visualization.
- **WGLMakie.jl** is the preferred backend for interactive visualization on the *web browser*.
- **CairoMakie.jl** is the preferred backend for *publication-quality* static visualization.

In this book, we use **CairoMakie.jl**:

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

::: {.callout-note}

We import the backend as `Mke` to avoid polluting the Julia
session with names from the visualization stack.

:::

Makie.jl provides a plot recipe system developed after
[Plots.jl](https://github.com/JuliaPlots/Plots.jl) by @Breloff2023,
which enables automatic visualization of custom Julia types.
The GeoStats.jl framework is integrated with this system,
and provides powerful visualization functions for geospatial
data.

Julia will automatically trigger the compilation of these
visualization functions whenever GeoStats.jl and Makie.jl
are loaded in the same session:

```{julia}
using GeoStats
```

## viz/viz!

The main visualization function that the framework provides is the
`viz`/`viz!` function. The `viz` function *creates a scene* and displays
geometries within a geospatial domain. On the other hand, the `viz!`
function adds more geometries to an *existing scene*.

Let's create a small geotable over a Cartesian grid for illustration purposes:

```{julia}
img = georef((A=rand(10, 10), B=rand(10, 10)))
```

::: {.callout-note}

The `georef` function creates a `CartesianGrid` starting at the origin
whenever the domain is omitted. The size of the grid is taken as the
size of the first array in the named tuple:

```{julia}
img.geometry
```

:::

::: {.callout-note}

## Tip for all users

To create a named tuple with a single key in Julia, we need an
extra comma after the key/value pair:

```{julia}
(A=rand(10, 10),)
```

or a semicolon before the key/value pair:

```{julia}
(; A=rand(10, 10))
```

:::

By default, all geometries are displayed with a single color:

```{julia}
viz(img.geometry)
```

We can pass any vector of [Colors.jl](https://github.com/JuliaGraphics/Colors.jl)
or numbers (automatically converted to colors) to the function via the `color`
option. It is common to pass colors from another column of the geotable:

```{julia}
viz(img.geometry, color = img.A)
```

but any vector with the same length can be passed:

```{julia}
viz(img.geometry, color = 1:length(img.A))
```

The `alpha` option can be used to control the transparency of each geometry
in the domain:

```{julia}
viz(img.geometry, color = img.B, alpha = rand(length(img.B)))
```

Other aesthetic options are available in the official documentation.

As another example, consider the visualization of data over a `GeometrySet`:

```{julia}
gset = GeometrySet([
  Triangle((12, 12), (15, 15), (12, 15)),
  Quadrangle((5, 12), (10, 15), (10, 18), (5, 15))
])

gis = georef((A=[0.1, 0.2], B=[0.3, 0.4]), gset)

viz(gis.geometry, color = gis.A)
```

We can create a scene with the geometries from the first geotable ("raster data"),
and then add the geometries from the second geotable ("vector data"):

```{julia}
viz(img.geometry)
viz!(gis.geometry)

# display current figure
Mke.current_figure()
```

Let's add an additional set of points and line segments to conclude the example:

```{julia}
viz!([Point(-20, -10), Point(-20, 0), Point(-40, 10)])
viz!([Segment((-40, -10), (0, 0)), Segment((-40, 0), (-20, 10))])

Mke.current_figure()
```

::: {.callout-note}

## Tip for all users

Makie.jl can set the aspect ratio of the axis after the visualization is created.
The following code can be used to adjust the aspect ratio for the data in the scene:

```{julia}
ax = Mke.current_axis()
ax.aspect = Mke.DataAspect()
Mke.current_figure()
```

Alternatively, it is possible to set the aspect ratio in the `viz/viz!` call directly:

```{julia}
viz(CartesianGrid(20, 10), color = 1:200, axis = (; aspect = Mke.DataAspect()))
```

:::

::: {.callout-note}

## Tip for advanced users

Makie.jl dispatches the `viz` and `viz!` functions whenever it encounters
a geospatial domain, a vector of geometries or a single geometry from Meshes.jl.
This means that you can replace `viz` with `Mke.plot` and `viz!` with `Mke.plot!`
in scripts and the result will be the same.

:::

::: {.callout-note}

## Tip for advanced users

In the case of `Mesh` domains, it is also possible to specify a color for each vertex
of the mesh. In this case, the `viz` and `viz!` functions fill in the domain with using
the interpolation routine from the graphics library:

```{julia}
grid = CartesianGrid(10, 10)

fig = Mke.Figure()
viz(fig[1,1], grid, color = 1:nelements(grid))
viz(fig[1,2], grid, color = 1:nvertices(grid))
fig
```

:::

The `viz` and `viz!` functions are aware of coordinate reference systems, which is a quite
unique feature of the framework, explored in the chapter [Map projections](06-projections.qmd).

Below is a preview of this feature in action:

```{julia}
using GeoIO

world = GeoIO.load("data/countries.geojson")

viz(world.geometry, color = 1:nrow(world))
```

## viewer

As geospatial data scientists we are often interested in quick inspection
of intermediate results from multivariate geostatistical analysis. Visualizing
all the variables manually with `viz`/`viz!` can be very time consuming.
To address this issue, the framework provides a basic `viewer` that displays
all variables stored in a geotable:

```{julia}
geotable = georef((A=rand(1000), B=rand(1000)), rand(Point, 1000))

viewer(geotable)
```

::: {.callout-note}

The `rand(Point, 1000)` function call creates `1000` random `Point`s in 3D space.

:::

It adds interactive elements to the scene, including a menu to select the
variable used as color, and a color bar that automatically updates upon
menu selection. The `viewer` will be particularly useful when we start to
work with geospatial transforms in **Part II** of the book. The pipe
operator (`|>`) in Julia will be preferred for reasons that will become
clear later:

```{julia}
geotable = georef((A=rand(1000), B=rand(1000)), CartesianGrid(10, 10, 10))

geotable |> viewer
```

::: {.callout-note}

The `viz/viz!` and the `viewer` automatically select color schemes
for variables based on data science traits from the
[DataScienceTraits.jl](https://github.com/JuliaML/DataScienceTraits.jl)
module. Additionally, the `viewer` recognizes units from the
[Unitful.jl](https://github.com/PainterQubits/Unitful.jl)
module:

```{julia}
geotable = georef((; A=[1,2,3,4]u"m"), CartesianGrid(2, 2))

geotable |> viewer
```

:::

## cbar

Any vector of Julia objects implementing the `getcolors`
function from the [Colorfy.jl](https://github.com/JuliaGraphics/Colorfy.jl)
module can be passed to the `color` option of `viz`/`viz!`.
The final vector of colors will be a function of the `colormap`
and `colorrange` options.

To reproduce this behavior in the colorbar, the framework provides
the `cbar` function. It is similar to Makie's colorbar, but addresses
various practical challenges with `missing` values, units, etc.

```{julia}
grid = CartesianGrid(2, 2)
vals = [1, missing, 3, 4]
cmap = "cividis"

fig = Mke.Figure()
viz(fig[1,1], grid, color = vals, colormap = cmap)
cbar(fig[1,2], vals, colormap = cmap)
fig
```

We are now equipped with a set of visualization functions that can
really improve the speed at which we explore and analyze geospatial data.
These functions provide a consistent set of aesthetic options that we
will cover in more detail with future examples.

Before we start learning the advanced features of the framework, we would
like to say a few words about integration with existing GIS technology.