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{2,Float64}
A B geometry
Continuous Continuous Quadrangle
[NoUnits] [NoUnits]
0.150724 0.623096 Quadrangle((0.0, 0.0), ..., (0.0, 1.0))
0.587183 0.01164 Quadrangle((1.0, 0.0), ..., (1.0, 1.0))
0.741848 0.0420371 Quadrangle((2.0, 0.0), ..., (2.0, 1.0))
0.942143 0.751364 Quadrangle((3.0, 0.0), ..., (3.0, 1.0))
0.715103 0.839959 Quadrangle((4.0, 0.0), ..., (4.0, 1.0))
0.280057 0.64694 Quadrangle((5.0, 0.0), ..., (5.0, 1.0))
0.273594 0.842941 Quadrangle((6.0, 0.0), ..., (6.0, 1.0))
0.0753848 0.83328 Quadrangle((7.0, 0.0), ..., (7.0, 1.0))
0.69266 0.139107 Quadrangle((8.0, 0.0), ..., (8.0, 1.0))
0.971226 0.688456 Quadrangle((9.0, 0.0), ..., (9.0, 1.0))
⋮ ⋮ ⋮
Note

The georef function will create 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{2,Float64}
  minimum: Point(0.0, 0.0)
  maximum: Point(10.0, 10.0)
  spacing: (1.0, 1.0)
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.8548247150016492 0.7111251506993508 … 0.5792091709762329 0.3262983600857501; 0.3152671159374969 0.8055572673094858 … 0.24799866479785282 0.02056708190125156; … ; 0.8100722105251407 0.10974406814872684 … 0.638797355102476 0.9261368750069429; 0.41846936491885844 0.0017386844094416931 … 0.19820150227370403 0.050802445040011746],)

or a semicolon before the key/value pair:

(; A=rand(10, 10))
(A = [0.6814769881089838 0.9316879481007608 … 0.6490089889723978 0.9884589653709208; 0.7609969293071576 0.10042365589058455 … 0.947472865604312 0.4653743853377539; … ; 0.09412350312241224 0.6127788244058747 … 0.628253522737035 0.35528580174484126; 0.9896436172069042 0.17222464484106748 … 0.23695450349466274 0.11414404072072448],)

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/ND0gA/src/scenes.jl:220

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/ND0gA/src/scenes.jl:220

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/ND0gA/src/scenes.jl:220

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/ND0gA/src/scenes.jl:220

Other aesthetic options are available in the official documentation. To really see the benefits of the framework, let’s load data from a GeoJSON file and visualize it:

using GeoIO

gis = GeoIO.load("data/geotable.geojson")

viz(gis.geometry, color = 1:4)
┌ 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/ND0gA/src/scenes.jl:220

Note

The “data” folder is stored on GitHub. Check the Preface for download instructions.

As already mentioned, the viz! function can be used to add more geometries to an existing scene. 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, color = 1:100)
viz!(gis.geometry, color = 1:4)

# 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/ND0gA/src/scenes.jl:220

Let’s add an additional set of points:

pts = [Point(-20, -10), Point(-20, 0), Point(-40, 10)]

viz!(pts, color = 1:3)

Mke.current_figure()

And a set of line segments to conclude the example:

seg = [Segment((-40, -10), (0, 0)), Segment((-40, 0), (-20, 10))]

viz!(seg, color = 1:2)

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/ND0gA/src/scenes.jl:220

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/ND0gA/src/scenes.jl:220

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(3, 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/ND0gA/src/scenes.jl:220

Note

The georef function will create a PointSet whenever its second argument is a Matrix of coordinates. In this case, the points are represented in the columns of the matrix:

geotable.geometry
1000 PointSet{3,Float64}
├─ Point(0.6583055509202254, 0.30928339367274393, 0.8175641246561246)
├─ Point(0.5902222198825507, 0.19733738935214606, 0.9679625824672432)
├─ Point(0.3923231985421536, 0.6775797803845387, 0.8660743803177298)
├─ Point(0.23475535732993547, 0.35823462446665044, 0.4335690672464808)
├─ Point(0.7661694789683925, 0.8992031351415611, 0.718519943057125)
⋮
├─ Point(0.12892418108124237, 0.602383962696433, 0.11721040560214246)
├─ Point(0.9461735759394971, 0.6391558377979862, 0.9038309683490976)
├─ Point(0.24729587036805822, 0.8774046396803857, 0.36005124204643557)
├─ Point(0.5460654001269168, 0.023764089412898648, 0.06802761467435803)
└─ Point(0.926172000698097, 0.9556836125907876, 0.7546844928124977)

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/ND0gA/src/scenes.jl:220

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/ND0gA/src/scenes.jl:220

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/ND0gA/src/scenes.jl:220

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.