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.7621609626765424, -0.1100425267249384]
 [0.6517294699899076, -1.1348767261038886]
 [2.330245431799896, -0.4283206400685533]
 [-0.7636666036298797, -0.49826806867872087]
 [0.36679844077485896, 0.9441768673231247]
 [-0.394878775997052, 1.3047224835057118]
 [-0.3977664937142015, -0.1560771790385185]
 [2.0949624367746065, -0.8982550219324746]
 [1.220310697169507, 0.8862556326643849]
 [0.5600980708032204, 0.33747692066811086]
 ⋮
 [-0.9811503983942305, -1.2590667725199889]
 [-0.3423133344670288, 0.11414792382705621]
 [0.5619032730364749, 0.03346563402933181]
 [-0.14009347227895683, -1.252781305384135]
 [0.025092297758777492, -0.07056021773120452]
 [0.5116711320673775, -0.2784132161335003]
 [-0.22269960649321632, 1.183986386088327]
 [-1.1435041755535809, 2.6722949484037946]
 [0.45699463416890584, -0.5772601214661285]

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}}:
 [-32.727026780936185, -2.2324837961991597]
 [27.98512238593496, -23.023770692247357]
 [100.06023450066759, -8.68953955338427]
 [-32.79167867760287, -10.108595491170252]
 [15.750245659772762, 19.15495417797659]
 [-16.956009176716876, 26.469510376147333]
 [-17.080007150495337, -3.166410146423896]
 [89.95723361716601, -18.223316394136468]
 [52.39987722157648, 17.979879216687277]
 [24.050489936870456, 6.846550869065103]
 ⋮
 [-42.13038575422496, -25.543271784458582]
 [-14.698860494287048, 2.3157718919946184]
 [24.128004930061994, 0.6789328446334463]
 [-6.015583378879438, -25.415755596399663]
 [1.0774578349734343, -1.4314878750007811]
 [21.97104766871627, -5.648298091331216]
 [-9.562672903388465, 24.020081149804877]
 [-49.101821806815295, 54.214087485364225]
 [19.623250682726457, -11.711143916082339]

finally, we can create the histogram.

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

<rect x="38" y="10" width="200" height="160.0" fill="none" stroke="#000" stroke-width="1" /> <text x="37.5" y="185.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >-166.0</text> <text x="237.5" y="185.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >195.0</text> <text x="18.75" y="170.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >-90.0</text> <text x="18.75" y="10.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >91.0</text> 

• edges: ([-166.0, -165.0, -164.0, -163.0, -162.0, -161.0, -160.0, -159.0, -158.0, -157.0 … 186.0, 187.0, 188.0, 189.0, 190.0, 191.0, 192.0, 193.0, 194.0, 195.0], [-90.0, -89.0, -88.0, -87.0, -86.0, -85.0, -84.0, -83.0, -82.0, -81.0 … 82.0, 83.0, 84.0, 85.0, 86.0, 87.0, 88.0, 89.0, 90.0, 91.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: 31.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.