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 GLMakie, GeoMakie
using FHist
import GLMakie: Point2d

First, we generate random points in a normal distribution:

random_data = randn(Point2d, 100_000)

100000-element Vector{Point{2, Float64}}:
 [0.2672400289850586, 0.7870431069957928]
 [0.5970819465023851, 0.843497004454412]
 [0.4417078402074984, -0.4389923845071178]
 [-0.4704379944924307, 0.1296014104240895]
 [1.6187291450272185, -0.4365283282484533]
 [-1.5041127750917016, 1.40663179904623]
 [0.14271103342928237, -0.02174060983283337]
 [-0.2656804527851519, -0.9681278243033627]
 [-0.8014762232000493, -2.35499925029525]
 [-0.09046608017143691, 0.7323485310989386]
 ⋮
 [-1.6717451897535722, -1.3152141061125093]
 [-0.8636939233213432, 0.15008076292151834]
 [-0.4923871555428924, -1.7596593819476012]
 [-1.3882977523879811, 0.3016021506206063]
 [1.4000935713305003, 1.7598553876154077]
 [-1.671761318508211, 1.0286917014055315]
 [-2.2704307492371134, -0.4784415207457223]
 [-1.6759540241266424, 1.647646285374616]
 [0.2963377722781525, -0.48968514683206626]

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{Point{2, Float64}}:
 [11.317793305019201, 16.099363966260434]
 [25.28681830464233, 17.25415693048509]
 [18.706621368292172, -8.979810780090368]
 [-19.923362546824055, 2.65106225874013]
 [68.55425794484105, -8.929407265735222]
 [-63.70017830254943, 28.77336336227136]
 [6.043907362353294, -0.44471514639555293]
 [-11.251744213717245, -19.803543250404182]
 [-33.943052122527455, -48.17269820898941]
 [-3.831298778046697, 14.9805588124934]
 ⋮
 [-70.79952276666187, -26.90336831573108]
 [-36.57801318191966, 3.069977749796217]
 [-20.852924143307764, -35.99472073995402]
 [-58.795334916792655, 6.169424206486313]
 [59.29489570924537, 35.998730134811076]
 [-70.80020583016834, 21.042407922502623]
 [-96.15425514963631, -9.786762771429368]
 [-70.97776970695138, 33.70343631768158]
 [12.550102122992566, -10.016756817569409]

finally, we can create the histogram.

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

• edges: ([-191.0, -190.0, -189.0, -188.0, -187.0, -186.0, -185.0, -184.0, -183.0, -182.0 … 161.0, 162.0, 163.0, 164.0, 165.0, 166.0, 167.0, 168.0, 169.0, 170.0], [-93.0, -92.0, -91.0, -90.0, -89.0, -88.0, -87.0, -86.0, -85.0, -84.0 … 79.0, 80.0, 81.0, 82.0, 83.0, 84.0, 85.0, 86.0, 87.0, 88.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: 34.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.