Histogram

You can use any Makie compatible recipe in GeoMakie, and this includes histograms! Until we get the datashader recipe working, you can also consider this a replacement for that if you use it with the correct nbins.

We use the excellent FHist.jl to create the histogram.

using CairoMakie, GeoMakie
using FHist
import CairoMakie: Point2d

First, we generate random points in a normal distribution:

random_data = randn(Point2d, 100_000)

100000-element Vector{GeometryBasics.Point{2, Float64}}:
 [0.9033349230673873, -2.015623399736624]
 [0.5632144931149654, 1.21468569982998]
 [1.849994266420347, -1.621050923384055]
 [1.7671260288143578, 2.267841126826689]
 [1.989652563127333, -0.2983794548561548]
 [-0.7985130788045014, -0.6728666698673912]
 [0.691269035414941, -0.23601057063101316]
 [-0.15467939135543077, 0.5394859642079955]
 [-0.4009825672734962, 0.9483359232722881]
 [-0.8482748548599394, 0.2717141445284994]
 ⋮
 [0.8496665628116443, -0.8252630346567763]
 [0.4717778686149525, -0.22394797320547455]
 [0.2860693976363315, -0.5643002535117538]
 [-0.4118121854594221, -0.7615201664191353]
 [-1.2559930882075876, -1.9539893877102088]
 [-1.9106415467783302, 1.2900799105818288]
 [0.6914153808695385, -0.6552127832976461]
 [-0.7659177989275551, 1.3765484054292598]
 [-0.15864359556300728, 1.07281634237227]

then, we rescale them to be within the lat/long bounds of the Earth:

xmin, xmax = extrema(first, random_data)
ymin, ymax = extrema(last, random_data)
latlong_data = random_data .* (Point2d(1/(xmax - xmin), 1/(ymax - ymin)) * Point2d(360, 180),)

100000-element Vector{GeometryBasics.Point{2, Float64}}:
 [38.63681875043913, -39.60549090589013]
 [24.089422131729147, 23.867664685981133]
 [79.12667974611496, -31.85243712315562]
 [75.58229660007325, 44.561380432611784]
 [85.10004816027993, -5.862932920581319]
 [-34.153451070908204, -13.221326353817739]
 [29.566482759712322, -4.637431035004813]
 [-6.61584032192033, 10.600495336622432]
 [-17.150550009983515, 18.634090966495673]
 [-36.28182746050393, 5.338979534337296]
 ⋮
 [36.34135263384826, -16.215800837764185]
 [20.178557846793915, -4.400410025672906]
 [12.235563116491754, -11.088077545424206]
 [-17.61374697525225, -14.963301194908718]
 [-53.72047073758577, -38.39442870889208]
 [-81.72064342344271, 25.34910448702695]
 [29.572742146634987, -12.874440698450266]
 [-32.75930822469892, 27.04814568032398]
 [-6.785394532154751, 21.080038015567407]

finally, we can create the histogram.

h = Hist2D((first.(latlong_data), last.(latlong_data)); nbins = (360, 180))

• edges: ([-167.0, -166.0, -165.0, -164.0, -163.0, -162.0, -161.0, -160.0, -159.0, -158.0 … 185.0, 186.0, 187.0, 188.0, 189.0, 190.0, 191.0, 192.0, 193.0, 194.0], [-98.0, -97.0, -96.0, -95.0, -94.0, -93.0, -92.0, -91.0, -90.0, -89.0 … 74.0, 75.0, 76.0, 77.0, 78.0, 79.0, 80.0, 81.0, 82.0, 83.0])
• bin counts: [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
• maximum count: 33.0
• total count: 100000.0

This is what the histogram looks like without any projection,

plot(h)

It's simple to plot to GeoAxis:

plot(h; axis = (; type = GeoAxis))

The projection can also be arbitrary!

plot(h; axis = (; type = GeoAxis, dest = "+proj=tissot +lat_1=60 +lat_2=65"))

---

This page was generated using Literate.jl.