plotly-bivariate-choropleth.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8d9747ca-6fa7-460c-81e4-b87e060aa3ef",
   "metadata": {},
   "source": [
    "![Header image: Example of a bivariate map of EU countries](https://raw.githubusercontent.com/yotkadata/plotly-bivariate-choropleth/main/img/header.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2c7cd65-fc1a-400b-87a6-53c923f0ceef",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Bivariate choropleth map using Plotly Choropleth Mapbox\n",
    "\n",
    "Based on the example by [_empet_](https://chart-studio.plotly.com/~empet/15191/texas-bivariate-choropleth-assoc/#/)\n",
    "\n",
    "**Choropleth** maps are a great way to visualize values in maps. But sometimes we not only want to show one variable, but the relationship of two. That's where **bivariate maps** come in. They show values of two variables using colors or symbols for each one of them. In the case of this example, we are using two colors that are blended together. That way, we can see areas were both variables are high or low and also those where one value is high and the other one is low.\n",
    "\n",
    "When searching for ways to create bivariate maps using Plotly, I was surprised to find little information on how to do it. What helped me a lot, was the example that [_empet_ published on plotly.com](https://chart-studio.plotly.com/~empet/15191/texas-bivariate-choropleth-assoc/#/). I built on that and tried to make it more generic and reuseable. I hope it worked. This notebook walks you through the whole process."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f5752c2-7ef6-4a31-bae6-c110592f64cb",
   "metadata": {},
   "source": [
    "## Part I: Generic code\n",
    "\n",
    "First, we declare functions and other things that are independet from the data in our examples. We start importing the necessary libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03a9658d-9b36-4ae6-9ae1-504eb3158207",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json, requests\n",
    "import plotly.graph_objs as go\n",
    "import os\n",
    "import math\n",
    "from urllib.parse import urlparse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84f62df-3eea-4842-80c4-affeec3d64e7",
   "metadata": {},
   "source": [
    "### Define functions\n",
    "\n",
    "In the next step, we define all the functions we will use later on. The 'core' functions for bivariate maps are `set_interval_value()` and `prepare_df()`. Both of them took the functionality from empet's example. But the latter works a little different and, in my opinion, is simpler than the original."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34cad7dd-ffb7-4c21-b4be-6cee32d21feb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function to set default variables\n",
    "\"\"\"\n",
    "\n",
    "def conf_defaults():\n",
    "    # Define some variables for later use\n",
    "    conf = {\n",
    "        'plot_title': 'Bivariate choropleth map using Ploty',  # Title text\n",
    "        'plot_title_size': 20,  # Font size of the title\n",
    "        'width': 1000,  # Width of the final map container\n",
    "        'ratio': 0.8,  # Ratio of height to width\n",
    "        'center_lat': 0,  # Latitude of the center of the map\n",
    "        'center_lon': 0,  # Longitude of the center of the map\n",
    "        'map_zoom': 3,  # Zoom factor of the map\n",
    "        'hover_x_label': 'Label x variable',  # Label to appear on hover\n",
    "        'hover_y_label': 'Label y variable',  # Label to appear on hover\n",
    "        'borders_width': 0.5,  # Width of the geographic entity borders\n",
    "        'borders_color': '#f8f8f8',  # Color of the geographic entity borders\n",
    "\n",
    "        # Define settings for the legend\n",
    "        'top': 1,  # Vertical position of the top right corner (0: bottom, 1: top)\n",
    "        'right': 1,  # Horizontal position of the top right corner (0: left, 1: right)\n",
    "        'box_w': 0.04,  # Width of each rectangle\n",
    "        'box_h': 0.04,  # Height of each rectangle\n",
    "        'line_color': '#f8f8f8',  # Color of the rectagles' borders\n",
    "        'line_width': 0,  # Width of the rectagles' borders\n",
    "        'legend_x_label': 'Higher x value',  # x variable label for the legend \n",
    "        'legend_y_label': 'Higher y value',  # y variable label for the legend\n",
    "        'legend_font_size': 9,  # Legend font size\n",
    "        'legend_font_color': '#333',  # Legend font color\n",
    "    }\n",
    "\n",
    "    # Calculate height\n",
    "    conf['height'] = conf['width'] * conf['ratio']\n",
    "    \n",
    "    return conf\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Function to recalculate values in case width is changed\n",
    "\"\"\"\n",
    "def recalc_vars(new_width, variables, conf=conf_defaults()):\n",
    "    \n",
    "    # Calculate the factor of the changed width\n",
    "    factor = new_width / 1000\n",
    "    \n",
    "    # Apply factor to all variables that have been passed to th function\n",
    "    for var in variables:\n",
    "        if var == 'map_zoom':\n",
    "            # Calculate the zoom factor\n",
    "            # Mapbox zoom is based on a log scale. map_zoom needs to be set to value ideal for our map at 1000px.\n",
    "            # So factor = 2 ^ (zoom - map_zoom) and zoom = log(factor) / log(2) + map_zoom\n",
    "            conf[var] = math.log(factor) / math.log(2) + conf[var]\n",
    "        else:\n",
    "            conf[var] = conf[var] * factor\n",
    "\n",
    "    return conf\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Function to load GeoJSON file with geographical data of the entities\n",
    "\"\"\"\n",
    "\n",
    "def load_geojson(geojson_url, data_dir='data', local_file=False):\n",
    "    \n",
    "    # Make sure data_dir is a string\n",
    "    data_dir = str(data_dir)\n",
    "    \n",
    "    # Set name for the file to be saved\n",
    "    if not local_file:\n",
    "        # Use original file name if none is specified\n",
    "        url_parsed = urlparse(geojson_url)\n",
    "        local_file = os.path.basename(url_parsed.path)\n",
    "        \n",
    "    geojson_file = data_dir + '/' + str(local_file)\n",
    "\n",
    "    # Create folder for data if it does not exist\n",
    "    if not os.path.exists(data_dir):\n",
    "        os.makedirs(data_dir)\n",
    "\n",
    "    # Download GeoJSON in case it doesn't exist\n",
    "    if not os.path.exists(geojson_file):\n",
    "\n",
    "        # Make http request for remote file data\n",
    "        geojson_request = requests.get(geojson_url)\n",
    "\n",
    "        # Save file to local copy\n",
    "        with open(geojson_file, 'wb') as file:\n",
    "            file.write(geojson_request.content)\n",
    "\n",
    "    # Load GeoJSON file\n",
    "    geojson = json.load(open(geojson_file, 'r'))\n",
    "    \n",
    "    # Return GeoJSON object\n",
    "    return geojson\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Function that assigns a value (x) to one of three bins (0, 1, 2).\n",
    "The break points for the bins can be defined by break_a and break_b.\n",
    "\"\"\"\n",
    "\n",
    "def set_interval_value(x, break_1, break_2):\n",
    "    if x <= break_1: \n",
    "        return 0\n",
    "    elif break_1 < x <= break_2: \n",
    "        return 1\n",
    "    else: \n",
    "        return 2\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Function that adds a column 'biv_bins' to the dataframe containing the \n",
    "position in the 9-color matrix for the bivariate colors\n",
    "    \n",
    "Arguments:\n",
    "    df: Dataframe\n",
    "    x: Name of the column containing values of the first variable\n",
    "    y: Name of the column containing values of the second variable\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "def prepare_df(df, x='x', y='y'):\n",
    "    \n",
    "    # Check if arguments match all requirements\n",
    "    if df[x].shape[0] != df[y].shape[0]:\n",
    "        raise ValueError('ERROR: The list of x and y coordinates must have the same length.')\n",
    "    \n",
    "    # Calculate break points at percentiles 33 and 66\n",
    "    x_breaks = np.percentile(df[x], [33, 66])\n",
    "    y_breaks = np.percentile(df[y], [33, 66])\n",
    "    \n",
    "    # Assign values of both variables to one of three bins (0, 1, 2)\n",
    "    x_bins = [set_interval_value(value_x, x_breaks[0], x_breaks[1]) for value_x in df[x]]\n",
    "    y_bins = [set_interval_value(value_y, y_breaks[0], y_breaks[1]) for value_y in df[y]]\n",
    "    \n",
    "    # Calculate the position of each x/y value pair in the 9-color matrix of bivariate colors\n",
    "    df['biv_bins'] = [int(value_x + 3 * value_y) for value_x, value_y in zip(x_bins, y_bins)]\n",
    "    \n",
    "    return df\n",
    "   \n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Function to create a color square containig the 9 colors to be used as a legend\n",
    "\"\"\"\n",
    "\n",
    "def create_legend(fig, colors, conf=conf_defaults()):\n",
    "    \n",
    "    # Reverse the order of colors\n",
    "    legend_colors = colors[:]\n",
    "    legend_colors.reverse()\n",
    "\n",
    "    # Calculate coordinates for all nine rectangles\n",
    "    coord = []\n",
    "\n",
    "    # Adapt height to ratio to get squares\n",
    "    width = conf['box_w']\n",
    "    height = conf['box_h']/conf['ratio']\n",
    "    \n",
    "    # Start looping through rows and columns to calculate corners the squares\n",
    "    for row in range(1, 4):\n",
    "        for col in range(1, 4):\n",
    "            coord.append({\n",
    "                'x0': round(conf['right']-(col-1)*width, 4),\n",
    "                'y0': round(conf['top']-(row-1)*height, 4),\n",
    "                'x1': round(conf['right']-col*width, 4),\n",
    "                'y1': round(conf['top']-row*height, 4)\n",
    "            })\n",
    "\n",
    "    # Create shapes (rectangles)\n",
    "    for i, value in enumerate(coord):\n",
    "        # Add rectangle\n",
    "        fig.add_shape(go.layout.Shape(\n",
    "            type='rect',\n",
    "            fillcolor=legend_colors[i],\n",
    "            line=dict(\n",
    "                color=conf['line_color'],\n",
    "                width=conf['line_width'],\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper',\n",
    "            xanchor='right',\n",
    "            yanchor='top',\n",
    "            x0=coord[i]['x0'],\n",
    "            y0=coord[i]['y0'],\n",
    "            x1=coord[i]['x1'],\n",
    "            y1=coord[i]['y1'],\n",
    "        ))\n",
    "    \n",
    "        # Add text for first variable\n",
    "        fig.add_annotation(\n",
    "            xref='paper',\n",
    "            yref='paper',\n",
    "            xanchor='left',\n",
    "            yanchor='top',\n",
    "            x=coord[8]['x1'],\n",
    "            y=coord[8]['y1'],\n",
    "            showarrow=False,\n",
    "            text=conf['legend_x_label'] + ' 🠒',\n",
    "            font=dict(\n",
    "                color=conf['legend_font_color'],\n",
    "                size=conf['legend_font_size'],\n",
    "            ),\n",
    "            borderpad=0,\n",
    "        )\n",
    "        \n",
    "        # Add text for second variable\n",
    "        fig.add_annotation(\n",
    "            xref='paper',\n",
    "            yref='paper',\n",
    "            xanchor='right',\n",
    "            yanchor='bottom',\n",
    "            x=coord[8]['x1'],\n",
    "            y=coord[8]['y1'],\n",
    "            showarrow=False,\n",
    "            text=conf['legend_y_label'] + ' 🠒',\n",
    "            font=dict(\n",
    "                color=conf['legend_font_color'],\n",
    "                size=conf['legend_font_size'],\n",
    "            ),\n",
    "            textangle=270,\n",
    "            borderpad=0,\n",
    "        )\n",
    "    \n",
    "    return fig\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Function to create the map\n",
    "\n",
    "Arguments:\n",
    "    df: The dataframe that contains all the necessary columns\n",
    "    colors: List of 9 blended colors\n",
    "    x: Name of the column that contains values of first variable (defaults to 'x')\n",
    "    y: Name of the column that contains values of second variable (defaults to 'y')\n",
    "    ids: Name of the column that contains ids that connect the data to the GeoJSON (defaults to 'id')\n",
    "    name: Name of the column conatining the geographic entity to be displayed as a description (defaults to 'name')\n",
    "\"\"\"\n",
    "\n",
    "def create_bivariate_map(df, colors, geojson, x='x', y='y', ids='id', name='name', conf=conf_defaults()):\n",
    "    \n",
    "    if len(colors) != 9:\n",
    "        raise ValueError('ERROR: The list of bivariate colors must have a length eaqual to 9.')\n",
    "    \n",
    "    # Recalculate values if width differs from default\n",
    "    if not conf['width'] == 1000:             \n",
    "        conf = recalc_vars(conf['width'], ['height', 'plot_title_size', 'legend_font_size', 'map_zoom'], conf)\n",
    "        \n",
    "    # Prepare the dataframe with the necessary information for our bivariate map\n",
    "    df_plot = prepare_df(df, x, y)\n",
    "    \n",
    "    # Create the figure\n",
    "    fig = go.Figure(go.Choroplethmapbox(\n",
    "        geojson=geojson,\n",
    "        locations=df_plot[ids],\n",
    "        z=df_plot['biv_bins'],\n",
    "        marker_line_width=.5,\n",
    "        colorscale=[\n",
    "            [0/8, colors[0]],\n",
    "            [1/8, colors[1]],\n",
    "            [2/8, colors[2]],\n",
    "            [3/8, colors[3]],\n",
    "            [4/8, colors[4]],\n",
    "            [5/8, colors[5]],\n",
    "            [6/8, colors[6]],\n",
    "            [7/8, colors[7]],\n",
    "            [8/8, colors[8]],\n",
    "        ],\n",
    "        customdata=df_plot[[name, ids, x, y]],  # Add data to be used in hovertemplate\n",
    "        hovertemplate='<br>'.join([  # Data to be displayed on hover\n",
    "            '<b>%{customdata[0]}</b> (ID: %{customdata[1]})',\n",
    "            conf['hover_x_label'] + ': %{customdata[2]}',\n",
    "            conf['hover_y_label'] + ': %{customdata[3]}',\n",
    "            '<extra></extra>',  # Remove secondary information\n",
    "        ])\n",
    "    ))\n",
    "\n",
    "    # Add some more details\n",
    "    fig.update_layout(\n",
    "        title=dict(\n",
    "            text=conf['plot_title'],\n",
    "            font=dict(\n",
    "                size=conf['plot_title_size'],\n",
    "            ),\n",
    "        ),\n",
    "        mapbox_style='white-bg',\n",
    "        width=conf['width'],\n",
    "        height=conf['height'],\n",
    "        autosize=True,\n",
    "        mapbox=dict(\n",
    "            center=dict(lat=conf['center_lat'], lon=conf['center_lon']),  # Set map center\n",
    "            zoom=conf['map_zoom']  # Set zoom\n",
    "        ),\n",
    "    )\n",
    "\n",
    "    fig.update_traces(\n",
    "        marker_line_width=conf['borders_width'],  # Width of the geographic entity borders\n",
    "        marker_line_color=conf['borders_color'],  # Color of the geographic entity borders\n",
    "        showscale=False,  # Hide the colorscale\n",
    "    )\n",
    "\n",
    "    # Add the legend\n",
    "    fig = create_legend(fig, colors, conf)\n",
    "    \n",
    "    return fig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8a8f2c2-962c-4f2c-9f63-d576571597d1",
   "metadata": {},
   "source": [
    "### Define lists of blended colors \n",
    "\n",
    "To get bivariate maps working, we need a set of nine colors that are the result of blending two main colors. Here are some examples, but others can easily be added:\n",
    "\n",
    "1) \"pink-blue\" by [Joshua Stevens](http://www.joshuastevens.net/cartography/make-a-bivariate-choropleth-map/)\n",
    "2) \"teal-red\" by [Joshua Stevens](http://www.joshuastevens.net/cartography/make-a-bivariate-choropleth-map/)\n",
    "3) \"blue-orange\" from the [ArcGIS Website](https://pro.arcgis.com/de/pro-app/latest/help/mapping/layer-properties/bivariate-colors.htm)\n",
    "\n",
    "![Three examples of color sets](https://raw.githubusercontent.com/yotkadata/plotly-bivariate-choropleth/main/img/colors.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae18c48-a070-405f-b198-3c223458f58f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define sets of 9 colors to be used\n",
    "# Order: bottom-left, bottom-center, bottom-right, center-left, center-center, center-right, top-left, top-center, top-right\n",
    "color_sets = {\n",
    "    'pink-blue':   ['#e8e8e8', '#ace4e4', '#5ac8c8', '#dfb0d6', '#a5add3', '#5698b9', '#be64ac', '#8c62aa', '#3b4994'],\n",
    "    'teal-red':    ['#e8e8e8', '#e4acac', '#c85a5a', '#b0d5df', '#ad9ea5', '#985356', '#64acbe', '#627f8c', '#574249'],\n",
    "    'blue-organe': ['#fef1e4', '#fab186', '#f3742d',  '#97d0e7', '#b0988c', '#ab5f37', '#18aee5', '#407b8f', '#5c473d']\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e902def9-3887-4e3f-9b99-1d6a4fe7f62e",
   "metadata": {},
   "source": [
    "## Part II: Prepare our data\n",
    "\n",
    "The next step is where we load our data and \"connect\" it to the map to actually create the map. Basically, we build a dataframe that contains at least one column with the **id** that will map the data to the geospatial data, two columns with the values of our two variables (**x** and **y**) and ideally a **name** or description of the geographic entities.\n",
    "\n",
    "We then load the geospatial data and pass both dataframe and the geospatial data to the function `create_bivariate_map()`. If the columns of the dataframe have exactly the names mentioned above (id, x, y, name), the function will work. If not, we have to pass the column names as keyword arguments to the function (see example 2 below).\n",
    "\n",
    "### First example: Housing and population ratios in Texas, USA (from empet's example)\n",
    "\n",
    "Because I wrote this notebook based on empet's example, we reproduce it. I am not exactly sure, though, what the data actually stand for. One variable is about housing, the other one about population. In both cases, a ratio of the values from two years (2010 and 2016) is calculated. The results are saved as the values of our x and y columns.\n",
    "\n",
    "#### Load data\n",
    "\n",
    "We load the data from CSV files and do some minor transformations. In this example we make sure that the data is saved in columns with the default names (id, name, x, y)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6caa25af-8809-4e95-872b-7914de42fbcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get data and write it to a dataframe containing\n",
    "    id: Id of the geographic entity (needs to be the same as references in the geospatial data)\n",
    "    name: Name of the geographic entity to be displayed as a description\n",
    "    x: Values of the first variable\n",
    "    y: Values of the second variable\n",
    "\"\"\"\n",
    "\n",
    "# Get data for first variable (x)\n",
    "url_x = 'https://raw.githubusercontent.com/empet/Datasets/master/Texas-housing.csv'\n",
    "df_x = pd.read_csv(url_x, usecols=['GEO.id2', 'GEO.display-label', 'hucen42010', 'huest72016'])\n",
    "\n",
    "# Rename columns\n",
    "df_x = df_x.rename(columns={'GEO.id2': 'id', 'GEO.display-label': 'name', 'hucen42010': 'x_2010', 'huest72016': 'x_2016'})\n",
    "\n",
    "# Calculate the ratio between two years\n",
    "df_x['x'] = df_x['x_2016'] / df_x['x_2010']\n",
    "\n",
    "\n",
    "# Get data for second variable (y)\n",
    "url_y = 'https://raw.githubusercontent.com/empet/Datasets/master/Texas-population.csv'\n",
    "df_y = pd.read_csv(url_y, usecols=['GEO.id2', 'rescen42010', 'respop72016'], skiprows=[1])\n",
    "\n",
    "# Rename columns\n",
    "df_y = df_y.rename(columns={'GEO.id2': 'id', 'rescen42010': 'y_2010', 'respop72016': 'y_2016'})\n",
    "\n",
    "# Calculate the ratio between two years\n",
    "df_y['y'] = df_y['y_2016'] / df_y['y_2010']\n",
    "\n",
    "\n",
    "# Merge both dataframes\n",
    "df = pd.merge(df_x, df_y[['id', 'y_2010', 'y_2016', 'y']], how=\"left\", left_on=['id'], right_on=['id']);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc0699c3-3dbc-4eaf-a69e-0e5e92ef8cee",
   "metadata": {},
   "source": [
    "#### Load geospatial data (GeoJSON)\n",
    "\n",
    "Next we load the GeoJSON file for the geospatial representation. In this case, we have to transform it a little, because it does not contain the necessary 'id' property."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e5e03f2-af48-4ef5-9846-a2e5bff35007",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define URL of the GeoJSON file\n",
    "geojson_url = 'https://data.texas.gov/api/geospatial/48ag-x9aa?method=export&format=GeoJSON'\n",
    "\n",
    "# Load GeoJSON file\n",
    "geojson = load_geojson(geojson_url, local_file='texas-counties.geojson')\n",
    "\n",
    "# Iterate over JSON to add necessary 'id' property (Alternatively, use featureidkey when building the fig)\n",
    "for i in range(len(geojson['features'])):\n",
    "    # Assign the county id to a new 'id' property for later linking to dataframe\n",
    "    # In our dataframe all ids start with the state code 48, so we'll add it\n",
    "    geojson['features'][i]['id'] = '48' + geojson['features'][i]['properties']['fips_code']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "913aea92-5849-4601-bb34-4fb920232b3c",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Create the bivariate map\n",
    "\n",
    "In the last step we create the map. To do so, we first override some variables defined in `conf_defaults()`. Then we pass the dataframe, the colorset, the geojson, and the configuration to `create_bivariate_map()`. The result looks pretty much the same as the one from the example:\n",
    "\n",
    "![Example of a bivariate map of Texas counties](https://raw.githubusercontent.com/yotkadata/plotly-bivariate-choropleth/main/img/example-map-texas-300.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2fdb3d1-1c6f-49f6-b44f-49d343e2cda7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load conf defaults\n",
    "conf = conf_defaults()\n",
    "\n",
    "# Override some variables\n",
    "conf['plot_title'] = 'Texas housing and population'\n",
    "conf['width'] = 1000  # Width of the final map container\n",
    "conf['ratio'] = 0.8  # Ratio of height to width\n",
    "conf['height'] = conf['width'] * conf['ratio']  # Width of the final map container\n",
    "conf['center_lat'] = 31.968599  # Latitude of the center of the map\n",
    "conf['center_lon'] = -99.901810  # Longitude of the center of the map\n",
    "conf['map_zoom'] = 5  # Zoom factor of the map\n",
    "conf['hover_x_label'] = 'Housing ratio'  # Label to appear on hover\n",
    "conf['hover_y_label'] = 'Population ratio'  # Label to appear on hover\n",
    "\n",
    "# Define settings for the legend\n",
    "conf['line_width'] = 0.5  # Width of the rectagles' borders\n",
    "conf['legend_x_label'] = 'Higher housing ratio'  # x variable label for the legend \n",
    "conf['legend_y_label'] = 'Higher pop ratio'  # y variable label for the legend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e73fda2f-cedc-4a40-b7c1-1efefeca6859",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create our bivariate map\n",
    "fig = create_bivariate_map(df, color_sets['pink-blue'], geojson, conf=conf)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da76e016-73f0-4327-a86a-22e7616b0bf8",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Second example: European data on population density and COVID-19 cases\n",
    "\n",
    "In a second example we use [data from Eurostat](https://ec.europa.eu/eurostat/databrowser/view/demo_r_d3dens/default/table?lang=en) on population density in the so-called NUTS regions and combine them with data about cumulated detected COVID-19 cases in the same regions. The latter comes from the [COVID19-European-Regional-Tracker by Asjad Naqvi](https://github.com/asjadnaqvi/COVID19-European-Regional-Tracker). I used it to visualize the spread of Sars-CoV-2 through Europe in my [Covid Waves](https://github.com/yotkadata/covid-waves) project. \n",
    "\n",
    "In this example we'll just use the cumulated detected cases in each NUTS region between 15th April 2020 and 15th April 2022 as our second variable. It has been discussed often, that the magnitude of the spread could be related to population density. It seems totally logical because if more people live in the same geographical space, more contacts between them are likely to happen which could lead to faster spread. Yet the data suggest that this assumption is not necessarily correct in every setting.\n",
    "\n",
    "Data on detected COVID-19 cases has many [caveats](https://github.com/yotkadata/covid-waves#caveats) and this example is not meant to draw conclusions but rather to show how to build a bivariate map. You'll also notice the missing data for many countries.\n",
    "\n",
    "#### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2a8f62e-27e2-487a-abba-227ccca5d969",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data for population density by NUTS region\n",
    "url_pop='https://ec.europa.eu/eurostat/api/dissemination/sdmx/2.1/data/DEMO_R_D3DENS$DEFAULTVIEW/?format=SDMX-CSV&startPeriod=2019&endPeriod=2019&i'\n",
    "df_pop = pd.read_csv(url_pop, usecols=['geo', 'OBS_VALUE'])\n",
    "\n",
    "# Rename column\n",
    "df_pop = df_pop.rename(columns={'OBS_VALUE': 'pop_2019'})\n",
    "\n",
    "\n",
    "# Load data for cumulated COVID-19 cases by NUTS region\n",
    "file_covid = 'data/covid-waves-data-cumulated.csv'\n",
    "\n",
    "# Import CSV\n",
    "df_covid = pd.read_csv(file_covid, usecols=['nuts_id', 'nuts_name', 'cumulated_pop'], header=0)\n",
    "\n",
    "# Merge dataframes\n",
    "df = pd.merge(df_pop, df_covid, how='left', left_on=['geo'], right_on=['nuts_id']);\n",
    "\n",
    "# Filter to use just NUTS level 3 (at least 3 characters)\n",
    "df = df.loc[df['geo'].str.len() > 3]\n",
    "\n",
    "# Remove rows with missing values\n",
    "df = df.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd359c01-20dd-4c75-89be-bfa53877aa06",
   "metadata": {},
   "source": [
    "#### Load geospatial data (GeoJSON)\n",
    "\n",
    "We'll use GeoJSON files provided by Eurostat. I included them in the Covid Waves project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d08f1a83-ed77-419b-9206-f03023b48e00",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define URL of the GeoJSON file\n",
    "geojson_url = 'https://github.com/yotkadata/covid-waves/raw/main/data/NUTS_RG_10M_2016_4326.geojson'\n",
    "\n",
    "# Load GeoJSON file\n",
    "geojson = load_geojson(geojson_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0638513a-b49a-4864-9f6c-5fefacbb9c4f",
   "metadata": {},
   "source": [
    "#### Create the bivariate map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7621b04d-21e5-4359-972e-3dd3cff8cef0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load conf defaults\n",
    "conf = conf_defaults()\n",
    "\n",
    "# Override some variables\n",
    "conf['plot_title'] = 'Population density and cumulated COVID-19 cases (per population) in EU-Countries'\n",
    "conf['hover_x_label'] = 'Population density'  # Label to appear on hover\n",
    "conf['hover_y_label'] = 'Cumulated cases per million'  # Label to appear on hover\n",
    "conf['width'] = 1000\n",
    "conf['center_lat'] = 48  # Latitude of the center of the map\n",
    "conf['center_lon'] = 9  # Longitude of the center of the map\n",
    "conf['map_zoom'] = 3.7  # Zoom factor of the map\n",
    "conf['borders_width'] = 0  # Width of the geographic entity borders\n",
    "\n",
    "# Define settings for the legend\n",
    "conf['top'] = 0.5  # Vertical position of the top right corner (0: bottom, 1: top)\n",
    "conf['right'] = 0.1  # Horizontal position of the top right corner (0: left, 1: right)\n",
    "conf['line_width'] = 0  # Width of the rectagles' borders\n",
    "conf['legend_x_label'] = 'Higher pop density'  # x variable label for the legend \n",
    "conf['legend_y_label'] = 'More COVID-19 cases'  # y variable label for the legend"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eb16c47-44b6-4489-a3c9-f06e8288428c",
   "metadata": {},
   "source": [
    "Since the column names of the dataframe differ from the defaults, this time we have to declare them when we call `create_bivariate_map()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ed50ac3-43b6-4ac9-a6a7-5b45a9ebf274",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create our bivariate map\n",
    "fig = create_bivariate_map(df, color_sets['teal-red'], geojson, x='pop_2019', y='cumulated_pop', ids='geo', name='nuts_name', conf=conf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3a61b75-d7fb-449d-8b60-90fd90337b8f",
   "metadata": {},
   "source": [
    "After creating the figure, we can still modify it. In this case we add another layer to draw the countries' borders ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bb40a2e-965c-417a-a851-2c65e992a1b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get geo data for country borders\n",
    "geojson_url = 'https://github.com/yotkadata/covid-waves/raw/main/data/CNTR_RG_10M_2016_4326.geojson'\n",
    "\n",
    "# Load GeoJSON file\n",
    "geo_countries = load_geojson(geojson_url)\n",
    "\n",
    "# Add country borders as thin lines\n",
    "fig.update_layout(\n",
    "    mapbox_layers = [{\n",
    "        'source': geo_countries,\n",
    "        'type': 'line',\n",
    "        'line': {'width': 0.5}\n",
    "    },],\n",
    ")\n",
    "\n",
    "# Add attribution\n",
    "fig.add_annotation(\n",
    "    xref='paper',\n",
    "    yref='paper',\n",
    "    x=0.01,\n",
    "    y=0,\n",
    "    showarrow=False,\n",
    "    text='<b>Data:</b> COVID19-European-Regional-Tracker/Eurostat, '\n",
    "         '<b>Graph:</b> Jan Kühn (https://yotka.org)'\n",
    ")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a7d3c6d-dfcb-489d-9b0a-09bb987adda3",
   "metadata": {},
   "source": [
    "#### Export map as image file\n",
    "\n",
    "We can also use Plotly to export the map as a PNG file or in other formats. In the code above, the function `recalc_vars()` makes sure that the resizing works fine, changing the font size, the height, and the zoom factor based on the width we set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5d6ac79-3aba-4acd-b444-f55130358239",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export  as png file\n",
    "\n",
    "# Set the width\n",
    "conf['width'] = 2000\n",
    "\n",
    "# Create figure\n",
    "fig = create_bivariate_map(df, color_sets['teal-red'], geojson, x='pop_2019', y='cumulated_pop', ids='geo', name='nuts_name', conf=conf)\n",
    "\n",
    "# Get geo data for country borders\n",
    "geojson_url = 'https://github.com/yotkadata/covid-waves/raw/main/data/CNTR_RG_10M_2016_4326.geojson'\n",
    "\n",
    "# Load GeoJSON file into\n",
    "geo_countries = load_geojson(geojson_url)\n",
    "\n",
    "# Add country borders as thin lines\n",
    "fig.update_layout(\n",
    "    mapbox_layers = [{\n",
    "        'source': geo_countries,\n",
    "        'type': 'line',\n",
    "        'line': {'width': 1}\n",
    "    },],\n",
    ")\n",
    "\n",
    "# Add attribution\n",
    "fig.add_annotation(\n",
    "    xref='paper',\n",
    "    yref='paper',\n",
    "    x=0.01,\n",
    "    y=0,\n",
    "    showarrow=False,\n",
    "    text='<b>Data:</b> COVID19-European-Regional-Tracker/Eurostat, '\n",
    "         '<b>Graph:</b> Jan Kühn (https://yotka.org)'\n",
    ")\n",
    "\n",
    "# Export as image\n",
    "fig.write_image('export.png', width=conf['width'], height=conf['height'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5aa020e-a9e0-4266-92df-eef0cc523cba",
   "metadata": {},
   "source": [
    "On the final map, we can see that there are some regions like Belgium, where high population density correlates with higher cumulated cases (dark red). In other regions, for example in France, there are high numbers of cumulated cases in lesser populated regions, too (teal color). Then again, there are also densely populated regions like in Spain, where density doesn't seem to correlate with cumulated cases. Of course, these findings are influenced by different national testing policies etc. and should be taken with a grain of salt.\n",
    "\n",
    "![Example of a bivariate map of EU countries](https://raw.githubusercontent.com/yotkadata/plotly-bivariate-choropleth/main/img/example-map-covid.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bivariate",
   "language": "python",
   "name": "bivariate"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}