d4c-datapkg-statistical/experiments/duckdb_census_of_population_lonboard.ipynb

In [1]:

import duckdb
import geopandas as gpd
import jenkspy
from lonboard import BitmapTileLayer, Map, PolygonLayer
from lonboard.colormap import apply_categorical_cmap
import numpy as np
import pyarrow as pa

con = duckdb.connect()
con.install_extension("spatial")
con.load_extension("spatial")

1. Total private dwellings and private dwellings per square kilometer at Dissemination Area geographic level

These values are from the 2021 Census of Population

In [2]:

"""
Vancouver CMA is geo.cma_dguid = '2021S0503933'
"""

con.execute("""
DROP TABLE IF EXISTS geo_data;
CREATE TABLE geo_data AS
SELECT
    geo.da_dguid,
    cop.count_total_1,
    cop.count_total_4,
    cop.count_total_6,
    cop.count_total_7,
    CASE
        WHEN cop.count_total_7 = 0.0 THEN 0
        WHEN cop.count_total_4 = 0 THEN 0
        WHEN cop.count_total_4 IS NULL THEN 0
        ELSE 
            ROUND(
                (cop.count_total_4 / cop.count_total_7), 2
            ) 
    END AS count_total_4_per_square_km,
    geo.geom
FROM
    'https://data.dataforcanada.org/processed/statistics_canada/census_of_population/2021/tabular/da_2021.parquet' AS cop,
    'https://data.dataforcanada.org/processed/statistics_canada/boundaries/2021/digital_boundary_files/da_2021.parquet' AS geo
WHERE geo.csd_name = 'Vancouver' AND cop.da_dguid = geo.da_dguid;

""")

con.execute("""
COPY geo_data TO './da_2021_characteristic.parquet' (FORMAT PARQUET);
""")

Out [2]:

<duckdb.duckdb.DuckDBPyConnection at 0x7f05770247b0>

In [3]:

characteristic_values = con.execute("SELECT DISTINCT count_total_4_per_square_km FROM geo_data").fetchall()

values = np.array([v[0] for v in characteristic_values])

# Compute Jenks breaks
num_classes = 5
breaks = jenkspy.jenks_breaks(values, n_classes=num_classes)

# Create a bin range mapping: (lower, upper) for each bin
bin_ranges = [(breaks[i], breaks[i+1]) for i in range(len(breaks)-1)]

# Create a function to get the range string for a value
def jenks_range(value) -> str:
    for i, (low, high) in enumerate(bin_ranges):
        if low <= value <= high:
            return f"{int(low)}-{int(high)}"
    return "unknown"


characteristic_df = gpd.read_parquet('./da_2021_characteristic.parquet')
characteristic_df['category'] = characteristic_df["count_total_4_per_square_km"].apply(lambda v: jenks_range(v))
characteristic_df['category'] = characteristic_df['category'].astype('category')

# Categories to colors
cmap = {}
colors = [
    [255, 255, 255],
    [255, 191.25, 191.25],
    [255, 127.50, 127.50],
    [255, 63.75, 63.75],
    [255, 0, 0]
]
for index, value in enumerate(sorted(characteristic_df['category'].unique(), key=lambda x: int(x.split('-')[0]))):
    cmap[value] = colors[index]

In [17]:

# OpenStreetMap

# We set `max_requests < 0` because `tile.openstreetmap.org` supports HTTP/2.
basemap = BitmapTileLayer(
    data="https://tile.openstreetmap.org/{z}/{x}/{y}.png",
    tile_size=256,
    max_requests=-1,
    min_zoom=0,
    max_zoom=19,
)

In [22]:

# Google Satellite
basemap = BitmapTileLayer(
    data="http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z}",
    tile_size=256,
    max_requests=-1,
    min_zoom=0,
    max_zoom=19,
)

In [6]:

get_color = apply_categorical_cmap(pa.array(characteristic_df['category']), cmap)

cop_layer = PolygonLayer.from_geopandas(gdf=characteristic_df,
                                        stroked=True,
                                        get_fill_color=get_color,
                                        get_line_color=[255, 255, 255],
                                        get_line_width=5,
                                        line_width_min_pixels=0.2,
                                        line_width_units="meters",
                                        opacity=0.4,
                                        auto_highlight = True
                                       )

m = Map([basemap, cop_layer])

m

Out [6]:

Map(custom_attribution='', layers=(BitmapTileLayer(data='http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={…

2. Percentage of people with income $100,000 and over

These values are from the 2021 Census of Population

In [24]:

con.execute("""
DROP TABLE IF EXISTS geo_data_2;
CREATE TABLE geo_data_2 AS
SELECT
    geo.da_dguid,
    cop.count_total_1,
    cop.count_total_155,
    cop.count_total_168,
    CASE
        WHEN cop.count_total_168 = 0.0 THEN 0
        WHEN cop.count_total_155 = 0 THEN 0
        WHEN cop.count_total_168 IS NULL THEN 0
        WHEN cop.count_total_155 IS NULL THEN 0
        ELSE 
            ((cop.count_total_168/cop.count_total_155) * 100) 
    END AS percentage_over_100k,
    geo.geom
FROM
    'https://data.dataforcanada.org/processed/statistics_canada/census_of_population/2021/tabular/da_2021.parquet' AS cop,
    'https://data.dataforcanada.org/processed/statistics_canada/boundaries/2021/digital_boundary_files/da_2021.parquet' AS geo
WHERE geo.cma_dguid = '2021S0503933' AND cop.da_dguid = geo.da_dguid;
""")

con.execute("""
COPY geo_data_2 TO './da_2021_characteristic_2.parquet' (FORMAT PARQUET);
""")

Out [24]:

<duckdb.duckdb.DuckDBPyConnection at 0x7f05770247b0>

In [25]:

characteristic_values = con.execute("SELECT DISTINCT percentage_over_100k FROM geo_data_2").fetchall()

values = np.array([v[0] for v in characteristic_values])

# Compute Jenks breaks
num_classes = 5
breaks = jenkspy.jenks_breaks(values, n_classes=num_classes)

# Create a bin range mapping: (lower, upper) for each bin
bin_ranges = [(breaks[i], breaks[i+1]) for i in range(len(breaks)-1)]

# Create a function to get the range string for a value
def jenks_range(value) -> str:
    for i, (low, high) in enumerate(bin_ranges):
        if low <= value <= high:
            return f"{int(low)}-{int(high)}"
    return "unknown"


characteristic_df = gpd.read_parquet('./da_2021_characteristic_2.parquet')
characteristic_df['category'] = characteristic_df["percentage_over_100k"].apply(lambda v: jenks_range(v))
characteristic_df['category'] = characteristic_df['category'].astype('category')

# Categories to colors
cmap = {}
colors = [
    [255, 255, 255],
    [255, 191.25, 191.25],
    [255, 127.50, 127.50],
    [255, 63.75, 63.75],
    [255, 0, 0]
]
for index, value in enumerate(sorted(characteristic_df['category'].unique(), key=lambda x: int(x.split('-')[0]))):
    cmap[value] = colors[index]

In [26]:

get_color = apply_categorical_cmap(pa.array(characteristic_df['category']), cmap)

cop_layer = PolygonLayer.from_geopandas(gdf=characteristic_df,
                                        stroked=True,
                                        get_fill_color=get_color,
                                        get_line_color=[255, 255, 255],
                                        get_line_width=5,
                                        line_width_min_pixels=0.2,
                                        line_width_units="meters",
                                        opacity=0.4,
                                        auto_highlight = True
                                       )

In [27]:

m = Map([basemap, cop_layer])

m

Out [27]:

Map(custom_attribution='', layers=(BitmapTileLayer(data='http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={…

In [ ]: