02 - Join NLCD counts per tract with ACS data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os, time\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import rasterio\n",
    "import rasterio.mask\n",
    "import rasterio.io\n",
    "import rasterio.errors\n",
    "import fiona\n",
    "import fiona.transform\n",
    "import shapely\n",
    "import shapely.geometry\n",
    "\n",
    "import collections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "NLCD_CLASS_VALS = [11, 12, 21, 22, 23, 24, 31, 41, 42, 43, 52, 71, 81, 82, 90, 95]\n",
    "NLCD_CLASS_NAMES = [\n",
    "    \"Open Water\",\n",
    "    \"Perennial Ice/Snow\",\n",
    "    \"Developed, Open Space\",\n",
    "    \"Developed, Low Intensity\",\n",
    "    \"Developed, Medium Intensity\",\n",
    "    \"Developed High Intensity\",\n",
    "    \"Barren Land (Rock/Sand/Clay)\",\n",
    "    \"Deciduous Forest\",\n",
    "    \"Evergreen Forest\",\n",
    "    \"Mixed Forest\",\n",
    "    \"Shrub/Scrub\",\n",
    "    \"Grassland/Herbaceous\",\n",
    "    \"Pasture/Hay\",\n",
    "    \"Cultivated Crops\",\n",
    "    \"Woody Wetlands\",\n",
    "    \"Emergent Herbaceous Wetlands\"\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "acs_fns = [\n",
    "    \"data/raw/acs5yr/acs_5yr_2011_data.csv\",\n",
    "    \"data/raw/acs5yr/acs_5yr_2013_data.csv\",\n",
    "    \"data/raw/acs5yr/acs_5yr_2016_data.csv\"\n",
    "]\n",
    "\n",
    "# These are created by `01 - NLCD + Census Tract intersection.ipynb`\n",
    "processed_nlcd_fns = [\n",
    "    \"data/processed/nlcd_2011_counts_per_output_tracts.npy\",\n",
    "    \"data/processed/nlcd_2013_counts_per_output_tracts.npy\",\n",
    "    \"data/processed/nlcd_2016_counts_per_output_tracts.npy\"\n",
    "]\n",
    "\n",
    "output_fns = [\n",
    "    \"data/processed/acs5yr/acs_5yr_2011_nlcd_joined_data.csv\",\n",
    "    \"data/processed/acs5yr/acs_5yr_2013_nlcd_joined_data.csv\",\n",
    "    \"data/processed/acs5yr/acs_5yr_2016_nlcd_joined_data.csv\"\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "72749\n",
      "72257\n",
      "data/raw/acs5yr/acs_5yr_2011_data.csv\n",
      "data/raw/acs5yr/acs_5yr_2013_data.csv\n",
      "data/raw/acs5yr/acs_5yr_2016_data.csv\n"
     ]
    }
   ],
   "source": [
    "with open(\"data/processed/output_tracts_geoids.txt\", \"r\") as f:\n",
    "    output_tracts_geoids = f.read().strip().split(\"\\n\")\n",
    "\n",
    "with open(\"data/processed/output_tracts_bad_geoids.txt\", \"r\") as f:\n",
    "    output_tracts_bad_geoids = set(f.read().strip().split(\"\\n\"))\n",
    "\n",
    "print(len(output_tracts_geoids))\n",
    "output_tracts_geoids_mask = [\n",
    "    geoid not in output_tracts_bad_geoids\n",
    "    for geoid in output_tracts_geoids\n",
    "]\n",
    "print(sum(output_tracts_geoids_mask))\n",
    "\n",
    "for i in range(len(acs_fns)):\n",
    "    print(acs_fns[i])\n",
    "    df = pd.read_csv(acs_fns[i], sep=\"\\t\", encoding=\"ISO-8859-1\", dtype={\n",
    "        \"Geo_FIPS\":str,\n",
    "        \"Geo_STATE\":str,\n",
    "        \"Geo_COUNTY\":str\n",
    "    })\n",
    "    df[\"Geo_StateCounty\"] = df[\"Geo_STATE\"]+df[\"Geo_COUNTY\"]\n",
    "    df.set_index(\"Geo_FIPS\", inplace=True)\n",
    "    \n",
    "    for geoid in output_tracts_geoids: # sanity checks\n",
    "        assert geoid in df.index\n",
    "        \n",
    "    df = df.loc[output_tracts_geoids] # only keep the geoids that are in our intersected census tract set\n",
    "        \n",
    "    counts = np.load(processed_nlcd_fns[i])\n",
    "    assert counts.shape[1] == len(NLCD_CLASS_NAMES) # sanity check\n",
    "    \n",
    "    #df = df.iloc[output_tracts_geoids_mask] # only keep the rows that have an intersection with NLCD\n",
    "    #counts = counts[output_tracts_geoids_mask] # same as previous line\n",
    "    \n",
    "    df = df.dropna(axis=1, how='all') # get rid of columns that are all n/a\n",
    "    \n",
    "    for j, class_name in enumerate(NLCD_CLASS_NAMES):\n",
    "        class_name = class_name.replace(\",\", \"\").replace(\" \", \"_\")\n",
    "        class_name = class_name.lower()\n",
    "        df[\"nlcd-\"+class_name] = counts[:,j]\n",
    "        \n",
    "    os.makedirs(os.path.dirname(output_fns[i]), exist_ok=True)\n",
    "    df.to_csv(output_fns[i], sep=\"\\t\", index_label=\"Geo_FIPS\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.5",
   "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}