2  Scientific visualization

The visualization ecosystem in Julia is evolving very quickly. Among the various visualization projects, Makie.jl by Danisch and Krumbiegel (2021) is the most advanced for scientific visualization.

Makie.jl is currently organized in backend modules:

In this book, we use CairoMakie.jl:

import CairoMakie as Mke
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 by Breloff (2023), 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:

using GeoStats

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

img = georef((A=rand(10, 10), B=rand(10, 10)))
100×3 GeoTable over 10×10 CartesianGrid
A B geometry
Continuous Continuous Quadrangle
[NoUnits] [NoUnits] 🖈 Cartesian{NoDatum}
0.92962 0.414363 Quadrangle((x: 0.0 m, y: 0.0 m), ..., (x: 0.0 m, y: 1.0 m))
0.400857 0.849582 Quadrangle((x: 1.0 m, y: 0.0 m), ..., (x: 1.0 m, y: 1.0 m))
0.620492 0.5517 Quadrangle((x: 2.0 m, y: 0.0 m), ..., (x: 2.0 m, y: 1.0 m))
0.766519 0.332543 Quadrangle((x: 3.0 m, y: 0.0 m), ..., (x: 3.0 m, y: 1.0 m))
0.153079 0.079491 Quadrangle((x: 4.0 m, y: 0.0 m), ..., (x: 4.0 m, y: 1.0 m))
0.982223 0.704892 Quadrangle((x: 5.0 m, y: 0.0 m), ..., (x: 5.0 m, y: 1.0 m))
0.18714 0.656235 Quadrangle((x: 6.0 m, y: 0.0 m), ..., (x: 6.0 m, y: 1.0 m))
0.0378243 0.289888 Quadrangle((x: 7.0 m, y: 0.0 m), ..., (x: 7.0 m, y: 1.0 m))
0.478947 0.0562625 Quadrangle((x: 8.0 m, y: 0.0 m), ..., (x: 8.0 m, y: 1.0 m))
0.448199 0.416113 Quadrangle((x: 9.0 m, y: 0.0 m), ..., (x: 9.0 m, y: 1.0 m))
⋮ ⋮ ⋮
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:

img.geometry
10×10 CartesianGrid
├─ minimum: Point(x: 0.0 m, y: 0.0 m)
├─ maximum: Point(x: 10.0 m, y: 10.0 m)
└─ spacing: (1.0 m, 1.0 m)
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:

(A=rand(10, 10),)
(A = [0.013176935096539055 0.26631715348159757 … 0.4413747884197794 0.509945239342051; 0.6873096741986273 0.9610479656463154 … 0.69849228766682 0.42180237888233696; … ; 0.41672288524546897 0.8368862780682905 … 0.6733883201159057 0.28817716996154963; 0.5187068483466101 0.629162672630517 … 0.6339632114085998 0.1043714045788382],)

or a semicolon before the key/value pair:

(; A=rand(10, 10))
(A = [0.9487255303727864 0.8191666884620601 … 0.5800690579356103 0.7696318845777067; 0.41841519504425084 0.6271038781835595 … 0.1482048313592934 0.036254284646667356; … ; 0.1824675418387599 0.3663345261170172 … 0.5964591494532647 0.8433675224981445; 0.8270225884627757 0.03396326268812433 … 0.5686773773922941 0.0892944970273628],)

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

viz(img.geometry)
┌ 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 pass any vector of 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:

viz(img.geometry, color = img.A)
┌ 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

but any vector with the same length can be passed:

viz(img.geometry, color = 1:length(img.A))
┌ 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 alpha option can be used to control the transparency of each geometry in the domain:

viz(img.geometry, color = img.B, alpha = rand(length(img.B)))
┌ 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 aesthetic options are available in the official documentation.

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

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)
┌ 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 create a scene with the geometries from the first geotable (“raster data”), and then add the geometries from the second geotable (“vector data”):

viz(img.geometry)
viz!(gis.geometry)

# display current figure
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

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

viz!([Point(-20, -10), Point(-20, 0), Point(-40, 10)])
viz!([Segment((-40, -10), (0, 0)), Segment((-40, 0), (-20, 10))])

Mke.current_figure()

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:

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:

viz(CartesianGrid(20, 10), color = 1:200, axis = (; aspect = Mke.DataAspect()))
┌ 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

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.

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:

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

Below is a preview of this feature in action:

using GeoIO

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

viz(world.geometry, color = 1:nrow(world))
┌ 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

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

geotable = georef((A=rand(1000), B=rand(1000)), rand(Point, 1000))

viewer(geotable)
┌ 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

Note

The rand(Point, 1000) function call creates 1000 random Points 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:

geotable = georef((A=rand(1000), B=rand(1000)), CartesianGrid(10, 10, 10))

geotable |> viewer
┌ 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

Note

The viz/viz! and the viewer automatically select color schemes for variables based on data science traits from the DataScienceTraits.jl module. Additionally, the viewer recognizes units from the Unitful.jl module:

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

geotable |> viewer
┌ 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

2.3 cbar

Any vector of Julia objects implementing the getcolors function from the 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.

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