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}}:
 [1.4699307579613206, -0.23902817876398513]
 [1.4275868175059316, -0.6541203321569774]
 [-0.4702316350351161, 0.8570563747514128]
 [2.006033549598473, -1.7991574574458076]
 [-0.6686415015400734, -0.615251989664724]
 [-0.43977084682626594, 0.13900035170727623]
 [-0.38250708608576783, 0.772009564225611]
 [0.825990116412275, 0.29344128838703043]
 [-0.5783295053623603, -0.5591652494754065]
 [0.5919637209964014, 1.5858924333564006]
 ⋮
 [1.1746163184796394, 0.3519668910732268]
 [-1.2490675569346097, 0.18913953403043876]
 [1.4505621582021768, 1.5813755968475751]
 [2.4316787929824955, -1.1181539808595118]
 [0.49176651685012335, 1.752502405309367]
 [-0.6458670616099822, 0.10306251936003169]
 [-1.8441761220482291, -0.6132240396279558]
 [-1.27129105055938, -1.6377499167007448]
 [-2.6972055520897036, -2.044073412220044]

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}}:
 [60.816902406012595, -4.8911705651305954]
 [59.064964581585386, -13.38509179647282]
 [-19.45535958157215, 17.537718500465207]
 [82.99761464840894, -36.81568442431712]
 [-27.66436767414211, -12.58973915773126]
 [-18.195075194927178, 2.8443275279469598]
 [-15.825844855678934, 15.797428052480933]
 [34.174507898485764, 6.004611672876524]
 [-23.92779992910776, -11.442050989220432]
 [24.491901848267567, 32.4516984968339]
 ⋮
 [48.59856535321632, 7.202205641275965]
 [-51.67891024606185, 3.8703123888403543]
 [60.0155461278933, 32.359271662919525]
 [100.60825724926997, -22.880489935304986]
 [20.34634359464549, 35.86099440030237]
 [-26.7220981943872, 2.1089411452184077]
 [-76.3009268475061, -12.548241751112979]
 [-52.5983848781682, -33.5128444982267]
 [-111.59415907311141, -41.82733499671856]

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" >-181.0</text> <text x="237.5" y="185.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >180.0</text> <text x="18.75" y="170.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >-96.0</text> <text x="18.75" y="10.0" fill="black" dominant-baseline="middle" text-anchor="middle" font-size="85%" font-family="sans-serif" >85.0</text> 

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