How to add a colorbar to a thematic map plot in Julia?

I would like to add a gradient (vertical) colorbar to my thematic map plot in Julia with values from the variable values. The map below is a population density choropleth.

Any suggestions how to add a gradient legend bar?

My code (including downloading en unzipping open data shapefiles):

using DataFrames
using Plots
using Shapefile
using ZipFile

# make directory downloads
downloads = joinpath(pwd(), "downloads")
if ~ispath(downloads)
    mkpath(downloads)
end

# make directory shapefiles
shapefiles = joinpath(pwd(), "shapefiles")
if ~ispath(shapefiles)
    mkpath(shapefiles)
end

# download shapefiles
url = "https://www.cbs.nl/-/media/cbs/dossiers/nederland-regionaal/wijk-en-buurtstatistieken/wijkbuurtkaart_2020_v1.zip"
name_zipfile = split(url, "/")[end]
path_zipfile = joinpath(pwd(), "downloads", name_zipfile)
if ~isfile(path_zipfile)
    download(url, path_zipfile)
end

# extract shapefiles
path_shapefile = joinpath(pwd(), "shapefiles", "gemeente_2020_v1.shp")
if ~isfile(path_shapefile)
    r = ZipFile.Reader(path_zipfile)
    for file in r.files
        open(joinpath(pwd(), "shapefiles", file.name), "w") do io
            write(io, read(file))
        end
    end
end

# read shapefile
table = Shapefile.Table(path_shapefile)
df = table |> DataFrame

# filter for land (i.e. not water)
row_filter = df.H2O .== "NEE"

# filter data and shapes
municipality_data = df[row_filter, :]
municipality_shape = Shapefile.shapes(table)[row_filter]

function normalize(array)
    """
    Normalize array to values between 0 and 1
    """
    return [(x - minimum(array))/(maximum(array) - minimum(array)) for x in array]
end

# select variable to plot
var = "BEV_DICHTH" # population density

# values to plot
values = municipality_data[:, var]
normalized_values = normalize(values)

# colors
colormap = :heat
colors = Array([cgrad(colormap)[value] for value in normalized_values])

# plot thematic map
p = plot(size=(500, 600), axis=false, ticks=false)
for i = 1:nrow(municipality_data)
    plot!(municipality_shape[i], color=colors[i])
end
p

Colorbar example:

x = y = 0:10
plot(x, y, (x,y)->x*y, st=:contourf, fill=(true, cgrad(:heat)))