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.

Answer 2

I am not sure if I understand, so please just tell me if I am mistaken. If I understood well, you can achieve what you want using numpy only (I am sure translating this to xarray will be easy):

# ---- snipped code already in the question -----
import numpy as np
import matplotlib.pyplot as plt

#  fake a choropleth in a dirty, fast way
height, width = 2880, 5120
choropleth = np.empty((height, width, 3,), dtype=np.uint8)
CHUNKS = 10
x_size = width // CHUNKS
for x_step, x in enumerate(range(0, width, width // CHUNKS)):
    y_size = height // CHUNKS
    for y_step, y in enumerate(range(0, height, height // CHUNKS)):
        choropleth[y: y+y_size, x: x+x_size] = (255-x_step*255//CHUNKS,
                                                0, y_step*255//CHUNKS)
plt.figure("Fake Choropleth")
plt.imshow(choropleth)

# Option 1: play with alpha only
outimage = np.empty((height, width, 4,), dtype=np.uint8)  # RGBA image
outimage[:, :, 3] = img  # Set alpha channel
outimage[:, :, :3] = choropleth  # Set color
plt.figure("Alpha filter only")
plt.imshow(outimage)

# Option 2: clear the empty points
outimage[img == 0, :3] = 0  # White. use 0 for black
plt.figure("Points erased")
plt.imshow(outimage[:,:,:3])  # change to 'outimage' to see the image with alpha

Results: Dummy choroplet

Alpha filtered figure

Black background, no alpha filter Note that the images might seem different because of matplotlib's antialiasing.