pythonrasteriodatashadergeoviews

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

  • 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 Choroplet

    Alpha filtered figure Only alpha filter

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