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.7974131098347723, 0.04415788239162358]
 [-0.5373964339507625, -0.4368780442419598]
 [0.13239892518290133, 0.6252416101947327]
 [0.8774756670893418, 0.07415792913979354]
 [0.5876868570485041, -2.4336325204281217]
 [0.3194848378970301, 1.1015228866700126]
 [-0.6240003283077714, 0.42885107785375604]
 [0.9422681354586031, 1.4230456261330702]
 [-0.8072876542810622, -1.8888307708861458]
 [1.022837037009336, 0.663614026494272]
 ⋮
 [0.21923802665202527, 0.21599881297505547]
 [0.11139936208054696, -0.6146587979487436]
 [0.660210935280807, -0.8837012910095953]
 [1.1746796289002195, -0.6602393666273839]
 [1.3562055110080429, 0.5694305084645944]
 [0.03449911567128472, 0.5650581221239018]
 [0.015211789037207582, 0.5671750941752667]
 [-0.00012756059025262277, -0.8362436518196229]
 [0.46811762903300136, 1.546252664664989]

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}}:
 [29.4226920926041, 0.8895786601396426]
 [-19.828680533071484, -8.801087465981198]
 [4.885212897810779, 12.595748793562851]
 [32.376814543306956, 1.4939419118390962]
 [21.68428036678657, -49.02652571957066]
 [11.788248647739934, 22.190630541262983]
 [-23.02416313331951, 8.639380934375442]
 [34.76750617257246, 28.667838058583378]
 [-29.787039853149597, -38.05127092584854]
 [37.740311552029716, 13.368776865319752]
 ⋮
 [8.08937409432516, 4.3513847184518095]
 [4.110377782085371, -12.38255369840851]
 [24.360250446550026, -17.80252511766981]
 [43.34295060153583, -13.30079295757279]
 [50.040834132950955, 11.471411248168428]
 [1.2729372584226355, 11.383327731208293]
 [0.5612796924195894, 11.42597500184289]
 [-0.004706689573904085, -16.846471502834]
 [17.272453502817072, 31.149894405510658]

finally, we can create the histogram.

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

• edges: ([-173.0, -172.0, -171.0, -170.0, -169.0, -168.0, -167.0, -166.0, -165.0, -164.0 … 179.0, 180.0, 181.0, 182.0, 183.0, 184.0, 185.0, 186.0, 187.0, 188.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: 35.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.