How can a choropleth map be combined with a shaded raster in Python?

I want to plot characteristics of areas on a map, but with very uneven population density, the larger tiles misleadingly attract attention. Think of averages (of test scores, say) by ZIP codes.

High-resolution maps are available to separate inhabited locales and even density within them. The Python code below does produce a raster colored according to the average such density for every pixel.

However, what I would really need is coloring from a choropleth map of the same area (ZIP codes of Hungary in this case) but the coloring affecting only points that would show up on the raster anyway. The raster could only determine the gamma of the pixel (or maybe height in some 3D analog). What is a good way to go about this?

A rasterio.mask.mask somehow?

(By the way, an overlay with the ZIP code boundaries would also be nice, but I have a better understanding of how that could work with GeoViews.)

import rasterio
import os
import datashader as ds
from datashader import transfer_functions as tf
import xarray as xr
from matplotlib.cm import viridis

# download a GeoTIFF from this location: https://data.humdata.org/dataset/hungary-high-resolution-population-density-maps-demographic-estimates
data_path = '~/Downloads/'
file_name = 'HUN_youth_15_24.tif'  # young people
file_path = os.path.join(data_path, file_name)
src = rasterio.open(file_path)
da = xr.open_rasterio(file_path)
cvs = ds.Canvas(plot_width=5120, plot_height=2880)
img = tf.shade(cvs.raster(da,layer=1), cmap=viridis)
ds.utils.export_image(img, "map", export_path=data_path, fmt=".png")

from How can a choropleth map be combined with a shaded raster in Python?