How to use a continuous color scale in a thematic map using Julia and Plots.jl

I would like to plot a thematic map (choropleth) using Plots.jl. The DataFrame df I use has 3 columns: geometry, value1 and value2. The column geometry is of element type Polygon (data was read using the Shapefile package). The columns value1 and value2 are of element type Int64 with values >= 0.

plot(df.geometry, color = :lightblue) creates a map with lightblue polygons without issues (and without colorbar to the side).

However, I would like the value2 column to color the polygons (shapes) based on a colormap/colorscheme like viridis with a corresponding colorbar to the side as legend.

Any suggestions how to do this?

plot(shp.geometry, palette = ColorSchemes.viridis) does plot a coloured map, however it is not clear on which variable this colours are based. I would like to specify that manually (value2 column) and have the colorbar added.

Here's the full code for reproducibility and possible soltuions.

# Load modules
using CSV
using DataFrames
using Dates
using Shapefile
using ZipFile
using Plots
using ColorSchemes

# Create directorties
for dir in ["data", "downloads", "shapefiles"]
    path = joinpath(pwd(), dir)
    if !ispath(path)
        mkpath(path)
    end
end

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


# Extract shape files
r = ZipFile.Reader(zip_loc)
for f in r.files
    file_name = split(f.name, "/")[end]
    if startswith(file_name, "gemeente")
        println("Extracting: $(file_name)")
        open(joinpath(pwd(), "shapefiles", file_name), "w") do io
            write(io, read(f))
        end
    end
end

# Read shapefiles
name_shapefile = "gemeente_2021_v1.shp"
path_shapefile = joinpath(pwd(), "shapefiles", name_shapefile)
table = Shapefile.Table(path_shapefile)
df = table |> DataFrame
row_filter = df.H2O .== "NEE"
df = df[row_filter, [:geometry]]

# Add columns value1 and value2
df.value1 = rand(1:100, nrow(df))
df.value2 = rand(1:5000, nrow(df))

# Plot map
plot(df.geometry, palette = ColorSchemes.viridis, axis=false, ticks=false, size=(500, 600))

Try

plot(df.geometry, fill = palette(:viridis), 
                  fill_z = reshape(df.value2, 1, nrow(df)), 
                  axis=false,
                  ticks=false,
                  size=(500, 600))