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.5055878152143495, -0.12613077058505626]
 [-0.025407429232147692, -0.3614521207856617]
 [1.107353031566484, 0.04793668156066028]
 [-1.523598035088602, -0.5083141012227388]
 [2.1654303069208627, -0.43133392500949086]
 [0.1070595058515665, 0.5059165965378665]
 [1.8710228802291373, -0.3103938662239685]
 [-0.5543673708047137, 1.7331001489367812]
 [-0.37453843239558743, 0.42280507770518355]
 [0.40757603540555, 0.1465845500474848]
 ⋮
 [1.1173636347441886, -0.3309948300573684]
 [0.5423392916778367, -1.2002535994742785]
 [-0.0282792028050363, 2.4531441362641675]
 [-0.8856044638475629, 0.23763765353799826]
 [1.6727075498272304, -0.15933275517703674]
 [-0.14421594791943276, 0.544154383138599]
 [0.7092768695868212, -0.20332489638844364]
 [-0.21040008769395016, -1.755431252450707]
 [-0.9198766966935891, 0.496272213017]

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}}:
 [-22.107168180174835, -2.588011607808536]
 [-1.1109569775190218, -7.416447865348529]
 [48.419758085522616, 0.9835878092503201]
 [-66.62035157316362, -10.42983237377503]
 [94.68490049992891, -8.850316220128095]
 [4.681249092491616, 10.380639223479951]
 [81.81173722441146, -6.368810125032045]
 [-24.240087144476178, 35.56057956467048]
 [-16.37694553892266, 8.675317243095055]
 [17.821537010537625, 3.0076920587150373]
 ⋮
 [48.85747846045034, -6.791510575410483]
 [23.714151273119022, -24.627378659029894]
 [-1.2365272136694208, 50.33478722781466]
 [-38.72365241850467, 4.875963279666095]
 [73.14026566206147, -3.2692658420271457]
 [-6.305939579589276, 11.165220457072493]
 [31.01360944708693, -4.171917681696175]
 [-9.199885724701849, -36.01877985030372]
 [-40.22222891231514, 10.182751139657231]

finally, we can create the histogram.

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

• edges: ([-184.0, -183.0, -182.0, -181.0, -180.0, -179.0, -178.0, -177.0, -176.0, -175.0 … 168.0, 169.0, 170.0, 171.0, 172.0, 173.0, 174.0, 175.0, 176.0, 177.0], [-84.0, -83.0, -82.0, -81.0, -80.0, -79.0, -78.0, -77.0, -76.0, -75.0 … 88.0, 89.0, 90.0, 91.0, 92.0, 93.0, 94.0, 95.0, 96.0, 97.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.