01 - NLCD + Census Tract intersection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This script expects two input data sources:\n",
    "\n",
    "- NLCD data rasters located at `data/raw/nlcd/`\n",
    "  - These can be downloaded at https://www.mrlc.gov/data\n",
    "  - In this script I load compressed `*.tif` files created from the raw NLCD `*.img` files\n",
    "- Shapefiles of all US census tracts located at `data/raw/shapes/`\n",
    "  - These can be downloaded at https://geodata.socialexplorer.com/ or https://www.census.gov/cgi-bin/geo/shapefiles/index.php\n"
   ]
  },
  {
   "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 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": "markdown",
   "metadata": {},
   "source": [
    "## Find the set of census tracts that remain the same over all of our years of data\n",
    "\n",
    "Here we assume that the shape of census tracts with the same GEOIDs don't change. I'm not sure if this is a safe assumption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape_fns = [\n",
    "    \"data/raw/shapes/TRACT_2011_US_SL140_Coast_Clipped.shp\",\n",
    "    \"data/raw/shapes/TRACT_2013_US_SL140_Coast_Clipped.shp\",\n",
    "    \"data/raw/shapes/TRACT_2016_US_SL140_Coast_Clipped.shp\"\n",
    "]\n",
    "nlcd_fns = [\n",
    "    \"data/raw/nlcd/nlcd_2011.tif\",\n",
    "    \"data/raw/nlcd/nlcd_2013.tif\",\n",
    "    \"data/raw/nlcd/nlcd_2016.tif\"\n",
    "]\n",
    "\n",
    "\n",
    "tracts_per_year = [\n",
    "    set()\n",
    "    for fn in shape_fns\n",
    "]\n",
    "\n",
    "for i, fn in enumerate(shape_fns):\n",
    "    with fiona.open(fn, \"r\") as f:\n",
    "        for line in f:\n",
    "            tracts_per_year[i].add(line[\"properties\"][\"GEOID\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "73837\n",
      "72790\n",
      "73799\n"
     ]
    }
   ],
   "source": [
    "for tracts in tracts_per_year:\n",
    "    print(len(tracts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "72749\n"
     ]
    }
   ],
   "source": [
    "common_tracts = tracts_per_year[0]\n",
    "for tracts in tracts_per_year[1:]:\n",
    "    common_tracts = common_tracts.intersection(tracts)\n",
    "print(len(common_tracts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load shapes for the set of census tracts we will use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sanity check, make sure all of our raster data uses the same coordinate system\n",
    "with rasterio.open(\"data/raw/nlcd/nlcd_2016.tif\", \"r\") as f:\n",
    "    target_crs = f.crs.to_string()\n",
    "    \n",
    "for fn in nlcd_fns:\n",
    "    with rasterio.open(fn, \"r\") as f:\n",
    "        assert target_crs == f.crs.to_string()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the census tract shapes from the latest set of definitions and convert them to the coordinate system used by NLCD \n",
    "census_tract_geoids = []\n",
    "census_tract_shapes = []\n",
    "with fiona.open(shape_fns[-1], \"r\") as f:\n",
    "    for line in f:\n",
    "        geoid = line[\"properties\"][\"GEOID\"]\n",
    "        if geoid in common_tracts:\n",
    "            geom = line[\"geometry\"]\n",
    "            geom = fiona.transform.transform_geom(f.crs, target_crs, geom)\n",
    "\n",
    "            census_tract_geoids.append(geoid)\n",
    "            census_tract_shapes.append(geom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write out a new shapefile with the set of shapes that we will be using\n",
    "os.makedirs(\"data/processed/\", exist_ok=True)\n",
    "with fiona.open(shape_fns[-1], 'r') as source:\n",
    "    with fiona.open(\n",
    "            \"data/processed/shapes/output_tracts.shp\", \"w\",\n",
    "            crs=source.crs,\n",
    "            driver=source.driver,\n",
    "            schema=source.schema,\n",
    "            ) as sink:\n",
    "\n",
    "        for f in source:\n",
    "            geoid = f[\"properties\"][\"GEOID\"]\n",
    "            if geoid in common_tracts:\n",
    "                sink.write(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mask the NLCD data with each shape we care about and count the values within the mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "NLCD_CLASS_VALS = [11, 12, 21, 22, 23, 24, 31, 41, 42, 43, 52, 71, 81, 82, 90, 95]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/raw/nlcd/nlcd_2011.tif\n",
      "5000 72749 886.6057615280151\n",
      "10000 72749 1387.7164342403412\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "bad_geoids = set() # keep track of all the census tracts that do not have any intersections\n",
    "\n",
    "for i, fn in enumerate(nlcd_fns):\n",
    "    print(fn)\n",
    "    \n",
    "    counts = np.zeros((len(census_tract_geoids), len(NLCD_CLASS_VALS)), dtype=int)\n",
    "    \n",
    "    tic = float(time.time())\n",
    "    with rasterio.io.MemoryFile(open(fn, \"rb\")) as memfile: # store the raster data in memory to avoid disk I/O\n",
    "\n",
    "        for j, (geoid, shape) in enumerate(zip(census_tract_geoids, census_tract_shapes)):\n",
    "            if j % 5000 == 0 and j != 0:\n",
    "                print(j, len(census_tract_geoids), time.time()-tic)\n",
    "                \n",
    "            with memfile.open() as f:\n",
    "                try:\n",
    "                    out_data, out_transform = rasterio.mask.mask(f, [shape], crop=True, all_touched=True)\n",
    "                    for k, val in enumerate(NLCD_CLASS_VALS):\n",
    "                        counts[j,k] = (out_data == val).sum()\n",
    "                    \n",
    "                    if counts[j].sum() == 0:\n",
    "                        bad_geoids.add(geoid)\n",
    "                except (rasterio.errors.WindowError, ValueError):\n",
    "                    bad_geoids.add(geoid)\n",
    "                    \n",
    "    np.save(\"data/processed/\" + os.path.basename(fn)[:-4] + \"_per_output_tracts.npy\", counts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"data/processed/output_tracts_geoids.txt\",\"w\")\n",
    "for geoid in census_tract_geoids:\n",
    "    f.write(\"%s\\n\" % (geoid))\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open(\"data/processed/output_tracts_bad_geoids.txt\",\"w\")\n",
    "for geoid in bad_geoids:\n",
    "    f.write(\"%s\\n\" % (geoid))\n",
    "f.close()"
   ]
  }
 ],
 "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
}