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

Question

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")

Answer 1

Datashader will let you combine data of many types into a common raster shape where you can do whatever making or filtering you like using xarray operations based on NumPy. E.g. you can render the choropleth as polygons, then mask out uninhabited regions. How to normalize by area is up to you, and could get very complex, but should be doable once you define precisely what you are intending to do. See the transform code at https://examples.pyviz.org/nyc_taxi/nyc_taxi.html for examples of how to do this, as in:

def transform(overlay):
    picks = overlay.get(0).redim(pickup_x='x', pickup_y='y')
    drops = overlay.get(1).redim(dropoff_x='x', dropoff_y='y')
    pick_agg = picks.data.Count.data
    drop_agg = drops.data.Count.data
    more_picks = picks.clone(picks.data.where(pick_agg>drop_agg))
    more_drops = drops.clone(drops.data.where(drop_agg>pick_agg))
    return (hd.shade(more_drops, cmap=['lightcyan', "blue"]) *
            hd.shade(more_picks, cmap=['mistyrose', "red"]))

picks = hv.Points(df, ['pickup_x',  'pickup_y'])
drops = hv.Points(df, ['dropoff_x', 'dropoff_y'])
((hd.rasterize(picks) * hd.rasterize(drops))).apply(transform).opts(
    bgcolor='white', xaxis=None, yaxis=None, width=900, height=500)

Here it's not really masking anything, but hopefully you can see how masking would work; just get some rasterized object then do a mathematical operation using some other rasterized object. Here the steps are all done in a function using HoloViews objects so that you can have a live interactive plot, but you would probably want to work out the approach using the more basic code at datashader.org where you only have to deal with xarray objects and not a HoloViews pipeline; you can then translate what you did for a single xarray into the HoloViews pipeline that would then allow full interactive usage with pan, zoom, axes, etc.