I am trying to figure out how to plot a heat map / density map of geolocations (lon, lat) and overlay it on a geographical map. I think I am nearly there using GeoMakie.jl, but I just can’t get the final step. I am able to produce this very nice map with points on it.
But now I would like to replace the individual points with some sort of density map. I found that KernelDensity.jl produces pretty much what I need, but it doesn’t plot correctly. Right now I have this
where I’ve fed the density estimator a list of (lon, lat) pairs. But the output lacks the information needed to plot it (datum, projection, etc.).
Is there a way I can attach the geographic information to the output of KernelDensity.kde such that the plots line up? Or, is there a better way to achieve what I want altogether?
Here is the code I used to make the above plots. (Other tips for making the code better are also appreciated.)
using CairoMakie
using DataFrames
using GADM
using GeoMakie
using GeoInterface
using GeometryOps
using KernelDensity
using Random
# Get USA state borders
usa = GADM.get("USA"; depth=1) |> DataFrame
# filter down to region of interest
state_list = ["Maine", "Massachusetts", "New Hampshire", "Vermont"]
state_borders_df = filter(row -> row["NAME_1"] in state_list, usa)
"""
Take a list of US states, download the county-level data, and combine into a single DataFrame.
Also, output the calculated centroid for each state
"""
function combine_states_df(state_list::Vector{String})
state_counties = DataFrame()
state_centroids = zeros(length(state_list), 2)
for (i, state) in enumerate(state_list)
state_df = GADM.get("USA", state; depth=1) |> DataFrame
append!(state_counties, state_df)
state_centroids[i, :] .= GeometryOps.centroid(state_df.geom)
end
return state_counties, state_centroids
end
state_counties_df, state_centroids = combine_states_df(state_list)
#= Create fake data
The following creates a random lon,lat pair for each state
centroid and only keeps the point if it is within the borders
of the states of interest.
=#
nrand = 1000
δ = 0.3
slons = zeros(nrand)
slats = zeros(nrand)
counter = 0
i = 1
while slons[end] == 0
keep = false
j = rand(1:size(state_centroids)[1])
lon = state_centroids[j, 1] .+ randn()*δ
lat = state_centroids[j, 2] .+ randn()*δ
if within((lon, lat), state_borders_df.geom[1]) keep = true end
if within((lon, lat), state_borders_df.geom[2]) keep = true end
if within((lon, lat), state_borders_df.geom[3]) keep = true end
if within((lon, lat), state_borders_df.geom[4]) keep = true end
if keep
slons[i] = lon
slats[i] = lat
i += 1
end
counter += 1
if counter % 100 == 0
println(counter)
end
# This just helps prevent infinite loops in case bad parameters are chosen.
if 2000 < counter
break
end
end
state_color = [:white]
county_line_width = 0.5
state_line_width = 1.0
# Create a nice figure with the points overlaid on the states
# the KDE output does not work on this plot
fig = Figure()
ga = GeoAxis(fig[1, 1]; xgridvisible=false, ygridvisible=false, xticklabelsvisible=false, yticklabelsvisible=false, dest="+proj=ortho +datum=WGS84 +lon_0=$(state_centroids[1, 1]) +lat_0=$(state_centroids[1, 2])")
poly!(
ga, state_counties_df.geom;
color = repeat(state_color, size(state_counties_df, 1)), # this could be any vector of length `size(ita_df, 1)`
strokecolor = :blue, strokewidth = county_line_width, shading = NoShading
)
poly!(
ga, state_borders_df.geom; alpha=0.5,
color = repeat(state_color, size(state_borders_df, 1)), # this could be any vector of length `size(ita_df, 1)`
strokecolor = :black, strokewidth = state_line_width, shading = NoShading
)
scatter!(ga, slons, slats; alpha=0.4)
plot!(ga, density; colormap=:Blues)
fig
# Calculate the density of points using KernelDensity.jl
density = kde((slons, slats))
# The KDE does work in GeoMakie, but the bounds are weird.
# How to put this plot into above plot??
fig = Figure()
ax = GeoAxis(fig[1, 1])
plot!(ax, density; colormap=:Blues)
fig