# Plot density map of points in GeoMakie
**URL:** https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107
**Category:** Geo
**Tags:** plotting
**Created:** [January 17, 2026, 1:58am UTC](https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107 "2026-01-17T01:58:52Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)
#### Post date: [January 17, 2026, 1:58am UTC](https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107/1 "2026-01-17T01:58:53Z")
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.)
```julia-auto
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
```
---
### Author:  [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [January 17, 2026, 12:08pm UTC](https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107/2 "2026-01-17T12:08:30Z")
Ah you’ve unfortunately encountered a recently introduced bug in Makie - `PlotSpec`s (which are what the KDE object gets converted to) are no longer transformed correctly.
I’d say the solution here is probably to manually deconstruct the KDE object into `density.x, density.y, density.density` and do a `GeoMakie.meshimage` or a `surface` plot.
I’m currently precompiling 😃 but will confirm that shortly.
---
### Author:  [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [January 17, 2026, 12:30pm UTC](https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107/3 "2026-01-17T12:30:40Z")
Here’s one way to do this, since your meshimage’s colormap is not transparent we have to put it beneath the polygons, whose colors also have to become transparent.
```julia
# 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 = :transparent, # this could be any vector of length `size(ita_df, 1)`
strokecolor = :blue, strokewidth = county_line_width,
)
poly!(
ga, state_borders_df.geom; alpha=0.5,
color = :transparent, # this could be any vector of length `size(ita_df, 1)`
strokecolor = :black, strokewidth = state_line_width,
)
mi_plt = meshimage!(ga, extrema(density.x), extrema(density.y), density.density; colormap=:Blues)
translate!(mi_plt, 0,0,-1) # get this below the outline polygons
fig
```
Otherwise you can also make a transparent colormap yourself, or find one, but they are a bit difficult to make correctly so I’d recommend just layering the plots in the right order.
---
### Author:  [@asinghvi17](https://discourse.julialang.org/u/asinghvi17)
#### Post date: [January 17, 2026, 12:36pm UTC](https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107/4 "2026-01-17T12:36:33Z")
Ah I just saw your comment about how you might want to color the polygons. In this case, here’s a way to do that, but you will have to play around with the colormap a bit.
CairoMakie is also not the best for projected data / heatmaps (although `heatmap` will probably work on a GeoAxis, but only in CairoMakie and not in the GL backends) so I’d suggest using e.g. GLMakie for this!
```julia
# First, design the colormap
blues_colormap = Makie.to_colormap(:Blues)
blues_with_alpha_colormap = map(enumerate(blues_colormap)) do (i, color)
RGBAf(color.r, color.g, color.b, (i) / length(blues_colormap))
end
# 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 = [RGBAf(rand(RGBf), .3) for _ in eachindex(state_counties_df.geom)], # this could be any vector of length `size(ita_df, 1)`
strokecolor = :blue, strokewidth = county_line_width,
)
poly!(
ga, state_borders_df.geom; alpha=0.5,
color = [RGBAf(rand(RGBf), .3) for _ in eachindex(state_borders_df.geom)], # this could be any vector of length `size(ita_df, 1)`
strokecolor = :black, strokewidth = state_line_width, shading = NoShading
)
mi_plt = meshimage!(ga, extrema(density.x), extrema(density.y), density.density; colormap=blues_with_alpha_colormap)
fig
```
---
### Author:  [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)
#### Post date: [January 18, 2026, 2:01am UTC](https://discourse.julialang.org/t/plot-density-map-of-points-in-geomakie/135107/5 "2026-01-18T02:01:09Z")
oooh yeah, that looks really nice! thanks! I think I can tweak it from here!