data_prep_radon_uranium_geodata.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "be869160",
   "metadata": {},
   "source": [
    "## Align, combine, subset EPA radon, uranium data and US Census GIS data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06819c51-5a2c-4c93-a38f-65090a52766a",
   "metadata": {},
   "source": [
    "### Notebook setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2fe0ea6e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/mitzi/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/libpysal/cg/alpha_shapes.py:39: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  def nb_dist(x, y):\n",
      "/Users/mitzi/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/libpysal/cg/alpha_shapes.py:165: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  def get_faces(triangle):\n",
      "/Users/mitzi/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/libpysal/cg/alpha_shapes.py:199: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  def build_faces(faces, triangles_is, num_triangles, num_faces_single):\n",
      "/Users/mitzi/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/libpysal/cg/alpha_shapes.py:261: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  def nb_mask_faces(mask, faces):\n"
     ]
    }
   ],
   "source": [
    "# import all libraries used in this notebook\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "# GIS libs\n",
    "import geopandas as gpd\n",
    "import libpysal as sa\n",
    "# plotting libs\n",
    "import matplotlib.pyplot as plt\n",
    "import splot as splt\n",
    "import plotnine as p9\n",
    "get_ipython().run_line_magic('matplotlib', 'inline')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68172a6e",
   "metadata": {},
   "source": [
    "### State/EPA Residiental Survey (SRRS) datasets\n",
    "\n",
    "The rawest form of the radon data was collected and archived by Phil Price and is available here:\n",
    "http://www.stat.columbia.edu/~gelman/arm/examples/radon_complete\n",
    "\n",
    "\n",
    "* The documentation is in file http://www.stat.columbia.edu/~gelman/arm/examples/radon_complete/SRRSdoc.pdf\n",
    "\n",
    "* There are 5 files, srrs1.dat through srrs5.dat  - but data is duplicated between them.\n",
    "\n",
    "* This directory also contains data from both national survey - NRRS - and state surveys - cf. https://link.springer.com/article/10.1007/BF02034901.   This is in a different format and is\n",
    "not used in Gelman and Hill analysis.\n",
    "\n",
    "* README notes that files are old backups, things may be missing.\n",
    "\n",
    "The combined de-duplicated SRRS dataset is in file  [srrs_all.csv](data/srrs_all.csv)\n",
    "\n",
    "*State counties and tribal lands*\n",
    "\n",
    "The SRRS dataset contains observations taken from Indian lands.\n",
    "The county-level information for these entries doesn't line up with US FIPS data -\n",
    "the names and county codes don't align.\n",
    "Indian lands have column 'STATE' code R5, R6, R7, RB, RC, RN.\n",
    "The regions cross state boundaries - for example,\n",
    "EPA region 5 covers Indian lands in MN, WI, and MI:\n",
    "https://www.epa.gov/sites/default/files/2015-08/documents/r5-tribal-land-map.pdf.\n",
    "\n",
    "Data in state counties in file [radon_all_states.csv](data/radon_all_states.csv).\n",
    "\n",
    "Data from Indian lands is in file [radon_indg_lands.csv](data/radon_indg.csv).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3798cf-5e30-4531-afad-d348fc5b1449",
   "metadata": {},
   "source": [
    "### US census county boundaries GIS files\n",
    "\n",
    "The US Census provides shapefiles for the US, including Alaska, Hawaii, and territories.  We can use these to visualize radon and uranium levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4f07daae-74a1-431a-8255-04bd1f29670d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3220\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>STATEFP</th>\n",
       "      <th>COUNTYFP</th>\n",
       "      <th>COUNTYNS</th>\n",
       "      <th>AFFGEOID</th>\n",
       "      <th>GEOID</th>\n",
       "      <th>NAME</th>\n",
       "      <th>LSAD</th>\n",
       "      <th>ALAND</th>\n",
       "      <th>AWATER</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>37</td>\n",
       "      <td>017</td>\n",
       "      <td>01026336</td>\n",
       "      <td>0500000US37017</td>\n",
       "      <td>37017</td>\n",
       "      <td>Bladen</td>\n",
       "      <td>06</td>\n",
       "      <td>2265887723</td>\n",
       "      <td>33010866</td>\n",
       "      <td>POLYGON ((-78.90200 34.83527, -78.79960 34.850...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>37</td>\n",
       "      <td>167</td>\n",
       "      <td>01025844</td>\n",
       "      <td>0500000US37167</td>\n",
       "      <td>37167</td>\n",
       "      <td>Stanly</td>\n",
       "      <td>06</td>\n",
       "      <td>1023370459</td>\n",
       "      <td>25242751</td>\n",
       "      <td>POLYGON ((-80.49737 35.20210, -80.29542 35.502...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>39</td>\n",
       "      <td>153</td>\n",
       "      <td>01074088</td>\n",
       "      <td>0500000US39153</td>\n",
       "      <td>39153</td>\n",
       "      <td>Summit</td>\n",
       "      <td>06</td>\n",
       "      <td>1069181981</td>\n",
       "      <td>18958267</td>\n",
       "      <td>POLYGON ((-81.68699 41.13596, -81.68495 41.277...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  STATEFP COUNTYFP  COUNTYNS        AFFGEOID  GEOID    NAME LSAD       ALAND  \\\n",
       "0      37      017  01026336  0500000US37017  37017  Bladen   06  2265887723   \n",
       "1      37      167  01025844  0500000US37167  37167  Stanly   06  1023370459   \n",
       "2      39      153  01074088  0500000US39153  39153  Summit   06  1069181981   \n",
       "\n",
       "     AWATER                                           geometry  \n",
       "0  33010866  POLYGON ((-78.90200 34.83527, -78.79960 34.850...  \n",
       "1  25242751  POLYGON ((-80.49737 35.20210, -80.29542 35.502...  \n",
       "2  18958267  POLYGON ((-81.68699 41.13596, -81.68495 41.277...  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shpfile = os.path.join('geo_data','cb_2018_us_county_20m', 'cb_2018_us_county_20m.shp')\n",
    "us_geodata = gpd.read_file(shpfile)\n",
    "# GEOID should be numeric\n",
    "us_geodata = us_geodata.astype({'GEOID': 'int32'}, copy=False)\n",
    "print(us_geodata.shape[0])\n",
    "us_geodata.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6548b6aa",
   "metadata": {},
   "source": [
    "### EPA/State Residential Radon Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2c67d0f1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(59395, 5)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>state</th>\n",
       "      <th>stfips</th>\n",
       "      <th>floor</th>\n",
       "      <th>activity</th>\n",
       "      <th>cntyfips</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AK</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>0.9</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AK</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1.1</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AK</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  state  stfips  floor  activity  cntyfips\n",
       "0    AK       2      0       0.9        20\n",
       "1    AK       2      0       1.1        20\n",
       "2    AK       2      0       1.0        20"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "us_radon = pd.read_csv(os.path.join('data','radon_all_counties.csv'),\n",
    "                     usecols=['state', 'stfips', 'floor', 'activity', 'cntyfips'],\n",
    "                     skipinitialspace=True,    # CSV file has spaces after delimiter, ignore them\n",
    "    ).convert_dtypes()\n",
    "print(us_radon.shape)\n",
    "us_radon.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c89b204",
   "metadata": {},
   "source": [
    "**datacleanup**\n",
    "\n",
    "Colorado and CT have data with cntyfips codes '0' and '999'.   Dropping for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d62fa643",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(57792, 5)\n"
     ]
    }
   ],
   "source": [
    "us_radon.drop(us_radon[us_radon.cntyfips==0].index, inplace=True)\n",
    "us_radon.drop(us_radon[us_radon.cntyfips==999].index, inplace=True)\n",
    "print(us_radon.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e78557c8-db6a-4bbd-a6a4-071e9d013ec8",
   "metadata": {},
   "source": [
    "### US county soil uranium levels\n",
    "\n",
    "Also distributed from Gelman website."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c7cdc9bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3111\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>stfips</th>\n",
       "      <th>ctfips</th>\n",
       "      <th>st</th>\n",
       "      <th>Uppm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>AL</td>\n",
       "      <td>1.78331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>AL</td>\n",
       "      <td>1.38323</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>AL</td>\n",
       "      <td>2.10105</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   stfips  ctfips  st     Uppm\n",
       "0       1       1  AL  1.78331\n",
       "1       1       3  AL  1.38323\n",
       "2       1       5  AL  2.10105"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "us_uranium = pd.read_csv(os.path.join('data','raw_uranium.csv'),\n",
    "                        usecols=['st', 'stfips', 'ctfips', 'Uppm'],\n",
    "                        skipinitialspace=True,\n",
    "                        ).drop_duplicates().convert_dtypes()\n",
    "print(us_uranium.shape[0])\n",
    "us_uranium.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93402d8a-4b4d-4e92-bd69-6fa008045cc3",
   "metadata": {},
   "source": [
    "### Join and merge tables using US FIPS codes\n",
    "\n",
    "To join or merge tables, we need to create a common key in both, then\n",
    "use the [DataFrame.merge](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.merge.html) method.\n",
    "\n",
    "\n",
    "We have three datasets:  SRRS survey data, soil uranium measurements, and geodata.\n",
    "All files use different capitalization and punctuation for county names.\n",
    "Therefore we rely on \n",
    "[FIPS code](https://transition.fcc.gov/oet/info/maps/census/fips/fips.txt),\n",
    "which uniquely identify geographic areas. \n",
    "The US census datasets have \"GEOID\" code, the first 2 digits of which are the state FIPS code, the last 3 are the county-level FIPS code.\n",
    "The other datasets have separate columns for stats and county codes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "41668fc1-a0c9-41a0-9ae7-9b44cb4127dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create merge column\n",
    "us_uranium['FIPS'] = us_uranium.stfips*1000 + us_uranium.ctfips\n",
    "us_radon['FIPS'] = us_radon.stfips*1000 + us_radon.cntyfips"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06c27212",
   "metadata": {},
   "source": [
    "### County level information:   uranium, number of homes in radon survey, census county name\n",
    "\n",
    "We create a new table which contains county-level information from across the three datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "29923eaf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3105\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>state</th>\n",
       "      <th>uranium</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>county</th>\n",
       "      <th>count</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AL</td>\n",
       "      <td>1.78331</td>\n",
       "      <td>1001</td>\n",
       "      <td>Autauga</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AL</td>\n",
       "      <td>1.38323</td>\n",
       "      <td>1003</td>\n",
       "      <td>Baldwin</td>\n",
       "      <td>31.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AL</td>\n",
       "      <td>2.10105</td>\n",
       "      <td>1005</td>\n",
       "      <td>Barbour</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  state  uranium  FIPS   county  count\n",
       "0    AL  1.78331  1001  Autauga    9.0\n",
       "1    AL  1.38323  1003  Baldwin   31.0\n",
       "2    AL  2.10105  1005  Barbour    9.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "us_counties = us_uranium.merge(us_geodata[['GEOID', 'NAME']],\n",
    "                               how='inner', left_on='FIPS', right_on='GEOID')\n",
    "\n",
    "homes = us_radon.value_counts(subset=['FIPS'], sort=False).to_frame().reset_index()\n",
    "homes.rename(columns={0:'homes'}, inplace=True)\n",
    "\n",
    "us_counties = us_counties.merge(homes, how='left', on='FIPS')\n",
    "us_counties.fillna(0, inplace=True)\n",
    "\n",
    "us_counties.drop(columns=['stfips', 'ctfips', 'GEOID'], inplace=True)\n",
    "us_counties.rename(columns={'st': 'state', 'NAME':'county', 'Uppm':'uranium'}, inplace=True)\n",
    "\n",
    "print(us_counties.shape[0])\n",
    "us_counties.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b5a51fd7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(87, 5)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "us_counties[us_counties.state=='MN'].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3a1505f",
   "metadata": {},
   "source": [
    "#### Put data on log scale\n",
    "\n",
    "Following Gelman and Hill chapter 4, section 4, we work with data on the log scale,\n",
    "for two reasons\n",
    "\n",
    "+ the outcome variable log_radon is always positive.\n",
    "+ it provides modeling flexibility.\n",
    "\n",
    "We know from geology that both radon measurements and soil uranium levels are always greater than zero,\n",
    "however a few radon measurements in the EPA dataset are 0.\n",
    "In order to be able to work with these measurements on the log scale, we replace 0 with 0.1,\n",
    "which corresponds to a low radon level (following Gelman and Hill)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4137f5ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>state</th>\n",
       "      <th>floor</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>radon</th>\n",
       "      <th>log_radon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AK</td>\n",
       "      <td>0</td>\n",
       "      <td>2020</td>\n",
       "      <td>0.9</td>\n",
       "      <td>-0.105361</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AK</td>\n",
       "      <td>0</td>\n",
       "      <td>2020</td>\n",
       "      <td>1.1</td>\n",
       "      <td>0.095310</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AK</td>\n",
       "      <td>0</td>\n",
       "      <td>2020</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  state  floor  FIPS  radon  log_radon\n",
       "0    AK      0  2020    0.9  -0.105361\n",
       "1    AK      0  2020    1.1   0.095310\n",
       "2    AK      0  2020    1.0   0.000000"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "us_radon['radon'] = us_radon.activity.apply(lambda x: x if x > 0.1 else 0.1)\n",
    "us_radon['log_radon'] = np.log(us_radon['radon'])\n",
    "us_radon.drop(columns=['activity', 'stfips', 'cntyfips'], inplace=True)\n",
    "us_radon.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "20f5742d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>state</th>\n",
       "      <th>uranium</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>county</th>\n",
       "      <th>count</th>\n",
       "      <th>log_uranium</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AL</td>\n",
       "      <td>1.78331</td>\n",
       "      <td>1001</td>\n",
       "      <td>Autauga</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.578471</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AL</td>\n",
       "      <td>1.38323</td>\n",
       "      <td>1003</td>\n",
       "      <td>Baldwin</td>\n",
       "      <td>31.0</td>\n",
       "      <td>0.324421</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AL</td>\n",
       "      <td>2.10105</td>\n",
       "      <td>1005</td>\n",
       "      <td>Barbour</td>\n",
       "      <td>9.0</td>\n",
       "      <td>0.742437</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  state  uranium  FIPS   county  count  log_uranium\n",
       "0    AL  1.78331  1001  Autauga    9.0     0.578471\n",
       "1    AL  1.38323  1003  Baldwin   31.0     0.324421\n",
       "2    AL  2.10105  1005  Barbour    9.0     0.742437"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "us_counties.uranium.fillna(0.1, inplace=True)\n",
    "us_counties['u'] = us_counties.uranium.apply(lambda x: x if x > 0.1 else 0.1)\n",
    "us_counties['log_uranium'] = np.log(us_counties['u'])\n",
    "us_counties.drop(columns=['u'], inplace=True)\n",
    "us_counties.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b137c20",
   "metadata": {},
   "source": [
    "### Restrict dataset to Minnesota\n",
    "\n",
    "In order to work with just the data from Minnesota, we use a \n",
    "use a conditional expression to [filter specific rows of a dataframe](https://pandas.pydata.org/docs/getting_started/intro_tutorials/03_subset_data.html#how-do-i-filter-specific-rows-from-a-dataframe), combined with operation [reset_index(drop=True)](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.reset_index.html?highlight=reset_index#pandas.DataFrame.reset_index) so that the rows are indexed starting from 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e5df7c39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>floor</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>radon</th>\n",
       "      <th>log_radon</th>\n",
       "      <th>county</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>27001</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>Aitkin</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>27001</td>\n",
       "      <td>2.2</td>\n",
       "      <td>0.788457</td>\n",
       "      <td>Aitkin</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>27001</td>\n",
       "      <td>2.9</td>\n",
       "      <td>1.064711</td>\n",
       "      <td>Aitkin</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   floor   FIPS  radon  log_radon  county\n",
       "0      0  27001    1.0   0.000000  Aitkin\n",
       "1      0  27001    2.2   0.788457  Aitkin\n",
       "2      0  27001    2.9   1.064711  Aitkin"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mn_radon = us_radon[us_radon['state']=='MN'].reset_index(drop=True)\n",
    "mn_radon.drop(columns=['state'], inplace=True)\n",
    "mn_radon = mn_radon.merge(us_counties[['FIPS', 'county']], on='FIPS')\n",
    "mn_radon = mn_radon.sort_values(by='county', axis=0).reset_index(drop=True)\n",
    "mn_radon.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "8643e3d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>uranium</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>county</th>\n",
       "      <th>count</th>\n",
       "      <th>log_uranium</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.502054</td>\n",
       "      <td>27001</td>\n",
       "      <td>Aitkin</td>\n",
       "      <td>4.0</td>\n",
       "      <td>-0.689048</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.428565</td>\n",
       "      <td>27003</td>\n",
       "      <td>Anoka</td>\n",
       "      <td>52.0</td>\n",
       "      <td>-0.847313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.892741</td>\n",
       "      <td>27005</td>\n",
       "      <td>Becker</td>\n",
       "      <td>3.0</td>\n",
       "      <td>-0.113459</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    uranium   FIPS  county  count  log_uranium\n",
       "0  0.502054  27001  Aitkin    4.0    -0.689048\n",
       "1  0.428565  27003   Anoka   52.0    -0.847313\n",
       "2  0.892741  27005  Becker    3.0    -0.113459"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mn_counties = us_counties[us_counties['state']=='MN'].reset_index(drop=True)\n",
    "mn_counties.drop(columns=['state'], inplace=True)\n",
    "mn_counties.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c5f573-e008-4aa4-bd15-74603df7edec",
   "metadata": {},
   "source": [
    "#### Unique county ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "20e68449",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>uranium</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>county</th>\n",
       "      <th>count</th>\n",
       "      <th>log_uranium</th>\n",
       "      <th>county_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.502054</td>\n",
       "      <td>27001</td>\n",
       "      <td>Aitkin</td>\n",
       "      <td>4.0</td>\n",
       "      <td>-0.689048</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.428565</td>\n",
       "      <td>27003</td>\n",
       "      <td>Anoka</td>\n",
       "      <td>52.0</td>\n",
       "      <td>-0.847313</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.892741</td>\n",
       "      <td>27005</td>\n",
       "      <td>Becker</td>\n",
       "      <td>3.0</td>\n",
       "      <td>-0.113459</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    uranium   FIPS  county  count  log_uranium  county_id\n",
       "0  0.502054  27001  Aitkin    4.0    -0.689048          1\n",
       "1  0.428565  27003   Anoka   52.0    -0.847313          2\n",
       "2  0.892741  27005  Becker    3.0    -0.113459          3"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# super clunky to index offset from 1\n",
    "mn_counties.reset_index(inplace=True)\n",
    "mn_counties['county_id'] = mn_counties.index + 1\n",
    "mn_counties.drop(columns=['index'], inplace=True)\n",
    "mn_counties.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d91947f1",
   "metadata": {},
   "source": [
    "Add county ids to radon data as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9d68e3a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>floor</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>radon</th>\n",
       "      <th>log_radon</th>\n",
       "      <th>county</th>\n",
       "      <th>county_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>27001</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>Aitkin</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>27001</td>\n",
       "      <td>2.2</td>\n",
       "      <td>0.788457</td>\n",
       "      <td>Aitkin</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>27001</td>\n",
       "      <td>2.9</td>\n",
       "      <td>1.064711</td>\n",
       "      <td>Aitkin</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   floor   FIPS  radon  log_radon  county  county_id\n",
       "0      0  27001    1.0   0.000000  Aitkin          1\n",
       "1      0  27001    2.2   0.788457  Aitkin          1\n",
       "2      0  27001    2.9   1.064711  Aitkin          1"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mn_radon = mn_radon.merge(mn_counties[['FIPS', 'county_id']], on='FIPS')\n",
    "mn_radon.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394c9602",
   "metadata": {},
   "source": [
    "**Save as CSV files**\n",
    "\n",
    "These files are already part of this notebook, therefore calls to the  [pandas.to_csv](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.to_csv.html?highlight=to_csv#pandas.DataFrame.to_csv) method have been commented out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b04aea68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment as needed\n",
    "# mn_radon.to_csv(os.path.join('data', 'mn_radon.csv'), index=False)\n",
    "# mn_counties.to_csv(os.path.join('data', 'mn_counties.csv'), index=False)\n",
    "\n",
    "# us_radon.to_csv(os.path.join('data', 'us_radon.csv'), index=False)\n",
    "# us_counties.to_csv(os.path.join('data', 'us_counties.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62e3a42d-cdb0-44e3-adf1-b033bd39ebd6",
   "metadata": {},
   "source": [
    "### Add GeoSpatial Information for US Counties\n",
    "\n",
    "If we want to build a spatial model which allows for local pooling of information between nearby counties, we need to create a neighbor graph over all the counties in Minnesota.\n",
    "\n",
    "**GIS Data**\n",
    "\n",
    "Geographic information systems (GIS) data is any item which has a geographic location, either a single point or a set of bounding polygons.  In order to manage, analyze, and visualize GIS data, we use specialized packages which can do the geographic math.  In this notebook we use the following packages:\n",
    "\n",
    "- GeoPandas - manages a set of GIS records in tabular format\n",
    "- libpysal - spatial analysis package which can analyze distance between locations\n",
    "\n",
    "Cartographic data (maps) are encoded as a set of records, one per map region.  The [shapefile format](https://en.wikipedia.org/wiki/Shapefile) is an open specification used to insure interoperatility among GIS software packages.  When items in a dataset contain location labels, it is necessary to obtain a set of shapefiles for the corresponding map.\n",
    "\n",
    "The shapefiles for US counties are available from the [US Census Bureau](https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html).\n",
    "For this analysis, we are using shapefiles where the boundary information is specified with the lowest possible resolution; this greatly speeds up analysis and plotting.\n",
    "These can be downloaded via URL: https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_20m.zip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "375d07b9-d9d9-4d3d-a62c-209dffe009df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'geopandas.geodataframe.GeoDataFrame'>\n",
      "RangeIndex: 3220 entries, 0 to 3219\n",
      "Data columns (total 3 columns):\n",
      " #   Column    Non-Null Count  Dtype   \n",
      "---  ------    --------------  -----   \n",
      " 0   STATEFP   3220 non-null   string  \n",
      " 1   GEOID     3220 non-null   int32   \n",
      " 2   geometry  3220 non-null   geometry\n",
      "dtypes: geometry(1), int32(1), string(1)\n",
      "memory usage: 63.0 KB\n"
     ]
    }
   ],
   "source": [
    "shpfile = os.path.join('geo_data','cb_2018_us_county_20m', 'cb_2018_us_county_20m.shp')\n",
    "us_geodata = gpd.read_file(shpfile).convert_dtypes()\n",
    "us_geodata = us_geodata.astype({'GEOID': 'int32'}, copy=False)\n",
    "us_geodata.drop(columns=['COUNTYFP', 'COUNTYNS', 'AFFGEOID', 'NAME', 'LSAD', 'ALAND', 'AWATER'], inplace=True)\n",
    "us_geodata.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "925d0dfe-4c00-43cb-8ca4-527dd63c185e",
   "metadata": {},
   "source": [
    "For this analysis, we restrict our analysis to the counties in Minnesota."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8ce6122a-7306-4596-bd88-35dbfdcf3462",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>geometry</th>\n",
       "      <th>FIPS</th>\n",
       "      <th>county</th>\n",
       "      <th>county_id</th>\n",
       "      <th>log_uranium</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>POLYGON ((-93.81041 46.25095, -93.77790 46.589...</td>\n",
       "      <td>27001</td>\n",
       "      <td>Aitkin</td>\n",
       "      <td>1</td>\n",
       "      <td>-0.689048</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>POLYGON ((-93.51007 45.41480, -93.01956 45.411...</td>\n",
       "      <td>27003</td>\n",
       "      <td>Anoka</td>\n",
       "      <td>2</td>\n",
       "      <td>-0.847313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>POLYGON ((-96.19467 47.15115, -96.06707 47.151...</td>\n",
       "      <td>27005</td>\n",
       "      <td>Becker</td>\n",
       "      <td>3</td>\n",
       "      <td>-0.113459</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            geometry   FIPS  county  \\\n",
       "0  POLYGON ((-93.81041 46.25095, -93.77790 46.589...  27001  Aitkin   \n",
       "1  POLYGON ((-93.51007 45.41480, -93.01956 45.411...  27003   Anoka   \n",
       "2  POLYGON ((-96.19467 47.15115, -96.06707 47.151...  27005  Becker   \n",
       "\n",
       "   county_id  log_uranium  \n",
       "0          1    -0.689048  \n",
       "1          2    -0.847313  \n",
       "2          3    -0.113459  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get MN subset, add county-level information\n",
    "mn_geodata = us_geodata[us_geodata['STATEFP']=='27'].copy()\n",
    "mn_geodata = mn_geodata.merge(mn_counties[['FIPS', 'county', 'county_id', 'log_uranium']],\n",
    "                               how='inner', left_on='GEOID', right_on='FIPS')\n",
    "mn_geodata.drop(columns=['STATEFP', 'GEOID'], inplace=True)\n",
    "mn_geodata = mn_geodata.sort_values(by='county_id', axis=0).reset_index(drop=True)\n",
    "mn_geodata.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6e95681-8f4c-4937-81b8-c47ac166d20a",
   "metadata": {},
   "source": [
    "To check our work, we plot the counties in Minnesota, colored by log uranium level and labeled by county name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "98a35a71-0780-4052-a82a-05061eb7f1fa",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/db/4jnggnf549s42z50bd61jskm0000gq/T/ipykernel_33423/2611589911.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
      "\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "'Figure' object has no attribute 'set_layout_engine'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "File \u001b[0;32m~/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/IPython/core/formatters.py:708\u001b[0m, in \u001b[0;36mPlainTextFormatter.__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    701\u001b[0m stream \u001b[38;5;241m=\u001b[39m StringIO()\n\u001b[1;32m    702\u001b[0m printer \u001b[38;5;241m=\u001b[39m pretty\u001b[38;5;241m.\u001b[39mRepresentationPrinter(stream, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mverbose,\n\u001b[1;32m    703\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmax_width, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnewline,\n\u001b[1;32m    704\u001b[0m     max_seq_length\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mmax_seq_length,\n\u001b[1;32m    705\u001b[0m     singleton_pprinters\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msingleton_printers,\n\u001b[1;32m    706\u001b[0m     type_pprinters\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mtype_printers,\n\u001b[1;32m    707\u001b[0m     deferred_pprinters\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mdeferred_printers)\n\u001b[0;32m--> 708\u001b[0m \u001b[43mprinter\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpretty\u001b[49m\u001b[43m(\u001b[49m\u001b[43mobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    709\u001b[0m printer\u001b[38;5;241m.\u001b[39mflush()\n\u001b[1;32m    710\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m stream\u001b[38;5;241m.\u001b[39mgetvalue()\n",
      "File \u001b[0;32m~/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/IPython/lib/pretty.py:410\u001b[0m, in \u001b[0;36mRepresentationPrinter.pretty\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    407\u001b[0m                         \u001b[38;5;28;01mreturn\u001b[39;00m meth(obj, \u001b[38;5;28mself\u001b[39m, cycle)\n\u001b[1;32m    408\u001b[0m                 \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mcls\u001b[39m \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mobject\u001b[39m \\\n\u001b[1;32m    409\u001b[0m                         \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mcallable\u001b[39m(\u001b[38;5;28mcls\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__dict__\u001b[39m\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m__repr__\u001b[39m\u001b[38;5;124m'\u001b[39m)):\n\u001b[0;32m--> 410\u001b[0m                     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_repr_pprint\u001b[49m\u001b[43m(\u001b[49m\u001b[43mobj\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcycle\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    412\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m _default_pprint(obj, \u001b[38;5;28mself\u001b[39m, cycle)\n\u001b[1;32m    413\u001b[0m \u001b[38;5;28;01mfinally\u001b[39;00m:\n",
      "File \u001b[0;32m~/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/IPython/lib/pretty.py:778\u001b[0m, in \u001b[0;36m_repr_pprint\u001b[0;34m(obj, p, cycle)\u001b[0m\n\u001b[1;32m    776\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"A pprint that just redirects to the normal repr function.\"\"\"\u001b[39;00m\n\u001b[1;32m    777\u001b[0m \u001b[38;5;66;03m# Find newlines and replace them with p.break_()\u001b[39;00m\n\u001b[0;32m--> 778\u001b[0m output \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mrepr\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    779\u001b[0m lines \u001b[38;5;241m=\u001b[39m output\u001b[38;5;241m.\u001b[39msplitlines()\n\u001b[1;32m    780\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m p\u001b[38;5;241m.\u001b[39mgroup():\n",
      "File \u001b[0;32m~/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/plotnine/ggplot.py:114\u001b[0m, in \u001b[0;36mggplot.__repr__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    110\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__repr__\u001b[39m(\u001b[38;5;28mself\u001b[39m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m \u001b[38;5;28mstr\u001b[39m:\n\u001b[1;32m    111\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    112\u001b[0m \u001b[38;5;124;03m    Print/show the plot\u001b[39;00m\n\u001b[1;32m    113\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 114\u001b[0m     figure \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdraw\u001b[49m\u001b[43m(\u001b[49m\u001b[43mshow\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m    116\u001b[0m     dpi \u001b[38;5;241m=\u001b[39m figure\u001b[38;5;241m.\u001b[39mget_dpi()\n\u001b[1;32m    117\u001b[0m     W \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mint\u001b[39m(figure\u001b[38;5;241m.\u001b[39mget_figwidth() \u001b[38;5;241m*\u001b[39m dpi)\n",
      "File \u001b[0;32m~/miniconda3/envs/jupyterlab/lib/python3.9/site-packages/plotnine/ggplot.py:241\u001b[0m, in \u001b[0;36mggplot.draw\u001b[0;34m(self, show)\u001b[0m\n\u001b[1;32m    239\u001b[0m     \u001b[38;5;66;03m# Artist object theming\u001b[39;00m\n\u001b[1;32m    240\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mtheme\u001b[38;5;241m.\u001b[39mapply()\n\u001b[0;32m--> 241\u001b[0m     \u001b[43mfigure\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mset_layout_engine\u001b[49m(PlotnineLayoutEngine(\u001b[38;5;28mself\u001b[39m))\n\u001b[1;32m    243\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfigure\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'Figure' object has no attribute 'set_layout_engine'"
     ]
    }
   ],
   "source": [
    "county_pts = mn_geodata.geometry.centroid.to_crs(mn_geodata.crs)\n",
    "(p9.ggplot()\n",
    " + p9.geom_map(data=mn_geodata, mapping=p9.aes(fill='log_uranium'), alpha=0.7)\n",
    " + p9.geom_text(mapping=p9.aes(x=county_pts.x, y=county_pts.y, label=list(mn_geodata.county)),\n",
    "                size=6, ha='center')\n",
    " + p9.theme(figure_size=(10, 8))\n",
    " + p9.scale_fill_cmap(cmap_name='viridis')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8b3fd47-c5a0-4c36-ae28-c27c278804f6",
   "metadata": {},
   "source": [
    "We need to find the set of pairs of adjacent counties; to do this we use `libpysal` which parses the bounding polygons to find all counties which have a common border of non-zero length, (\"Rook\" metric).  To check our work, we plot the resulting neighbor graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b68cc050-7e65-4f70-9782-0b0e52cc8cd7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# neighbors are counties which have a common line boundary\n",
    "from splot.libpysal import plot_spatial_weights\n",
    "from libpysal.weights.contiguity import Rook\n",
    "mn_nbs = Rook(mn_geodata['geometry'])\n",
    "plot_spatial_weights(mn_nbs, mn_geodata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29378142-32c5-43a9-95c3-853cb3bc6b45",
   "metadata": {},
   "source": [
    "As a sanity check, we get the edges list labeled by county name check that the upper rightmost county \"Lake\" is adjacent to \"Cook\" county."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39dcaeea-2209-4013-b2ab-dde42b5f2fc0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mn_nbs_names = Rook(mn_geodata['geometry'], ids=mn_geodata['county'].tolist())\n",
    "mn_nbs_edges_names =  mn_nbs_names.to_adjlist(remove_symmetric=True).reset_index(drop=True)\n",
    "mn_nbs_edges_names[mn_nbs_edges_names['focal']=='Cook']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efe44d25-434a-43aa-a47a-4387bde6a37c",
   "metadata": {},
   "source": [
    "We extract neighbor relationships as an edgelist and output this as a json file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62c01420-75aa-45ee-b741-41e6c0757582",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mn_nbs_edges =  mn_nbs.to_adjlist(remove_symmetric=True).reset_index(drop=True)\n",
    "mn_nbs_edges.head(3)\n",
    "\n",
    "node1 = (mn_nbs_edges['focal'] + 1).tolist()\n",
    "node2 = (mn_nbs_edges['neighbor'] + 1).tolist()\n",
    "\n",
    "mn_nbs_dict = { 'node1' : node1, 'node2' : node2, 'J_edges' : len(node1) }\n",
    "\n",
    "from cmdstanpy import write_stan_json\n",
    "write_stan_json(\"mn_nbs.json\", mn_nbs_dict)"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.16"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}