04.2 - Block Group Inundation Calculation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys, os, time, math, csv\n",
    "import itertools\n",
    "import collections\n",
    "import subprocess\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import fiona\n",
    "from fiona.transform import transform_geom\n",
    "\n",
    "from shapely.geometry import mapping, shape\n",
    "from shapely.ops import transform\n",
    "\n",
    "from rtree import index\n",
    "\n",
    "import rasterio\n",
    "import rasterio.mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge SLR Rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merged_slr_masks(slr_amount):\n",
    "\n",
    "    output_fn = \"data/processed/digital_coast/%dft.tif\" % (slr_amount)\n",
    "    input_fn_str = \"data/intermediate/digital_coast/slr_%dft/*.tif\" % (slr_amount)\n",
    "    \n",
    "    command = [\n",
    "        \"gdal_merge.py\",\n",
    "        \"-o\", output_fn,\n",
    "        \"-n\", \"-1\",\n",
    "        \"-a_nodata\", \"-1\",\n",
    "        \"-co\", \"COMPRESS=DEFLATE\",\n",
    "        \"-co\", \"PREDICTOR=1\",\n",
    "        \"-co\", \"TILED=NO\",\n",
    "        \"-co\", \"NUM_THREADS=ALL_CPUS\",\n",
    "        input_fn_str\n",
    "    ]\n",
    "    print(\"Running: \", \" \".join(command))\n",
    "    subprocess.call(\" \".join(command), shell=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running:  gdal_merge.py -o data/processed/digital_coast/0ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_0ft/*.tif\n",
      "Running:  gdal_merge.py -o data/processed/digital_coast/1ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_1ft/*.tif\n",
      "Running:  gdal_merge.py -o data/processed/digital_coast/2ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_2ft/*.tif\n",
      "Running:  gdal_merge.py -o data/processed/digital_coast/3ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_3ft/*.tif\n",
      "Running:  gdal_merge.py -o data/processed/digital_coast/4ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_4ft/*.tif\n",
      "Running:  gdal_merge.py -o data/processed/digital_coast/5ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_5ft/*.tif\n",
      "Running:  gdal_merge.py -o data/processed/digital_coast/6ft.tif -n -1 -a_nodata -1 -co COMPRESS=DEFLATE -co PREDICTOR=1 -co TILED=NO -co NUM_THREADS=ALL_CPUS data/intermediate/digital_coast/slr_6ft/*.tif\n"
     ]
    }
   ],
   "source": [
    "for slr_amount in range(7):\n",
    "    merged_slr_masks(slr_amount)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mask Merged SLR Rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_merged_slr_mask(slr_amount):\n",
    "    f = rasterio.open(\"data/processed/digital_coast/%dft.tif\" % (slr_amount), \"r\")\n",
    "    source_profile = f.profile\n",
    "    data = f.read()\n",
    "    f.close()\n",
    "    \n",
    "    dest_profile = source_profile.copy()\n",
    "    dest_profile[\"dtype\"] = rasterio.int16\n",
    "    #dest_profile[\"compress\"] = \"lzw\"\n",
    "    #dest_profile[\"nodata\"] = -1\n",
    "    #del dest_profile[\"transform\"]\n",
    "\n",
    "    dest_data = data.copy()\n",
    "    dest_data[dest_data>0] = 1\n",
    "    dest_data[dest_data<1] = 0\n",
    "    dest_data = dest_data.astype(np.int16)\n",
    "    \n",
    "    f = rasterio.open(\"data/processed/digital_coast/%dft_masked.tif\" % (slr_amount), \"w\", **dest_profile)\n",
    "    f.write(dest_data)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "for slr_amount in range(7):\n",
    "    process_merged_slr_mask(slr_amount)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Take Union of Masked Merged SLR Rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n"
     ]
    }
   ],
   "source": [
    "f = rasterio.open(\"data/processed/digital_coast/0ft_masked.tif\", \"r\")\n",
    "mask = f.read()\n",
    "f.close()\n",
    "\n",
    "for slr_amount in range(7):\n",
    "    print(slr_amount)\n",
    "    f = rasterio.open(\"data/processed/digital_coast/%dft_masked.tif\" % (slr_amount), \"r\")\n",
    "    profile = f.profile\n",
    "    mask = np.clip(mask + f.read(), None, 1)\n",
    "    f.close()\n",
    "\n",
    "    f = rasterio.open(\"data/processed/digital_coast/%dft_masked_union.tif\" % (slr_amount), \"w\", **profile)\n",
    "    f.write(mask)\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load block groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished loading 216330 block group geometries in 28.8010 seconds\n"
     ]
    }
   ],
   "source": [
    "load_time = float(time.time())\n",
    "\n",
    "block_group_geoids = []\n",
    "block_group_geoms = []\n",
    "f = fiona.open(\"data/processed/boundary_shapefiles/tl_2012_all_bg.shp\", \"r\")\n",
    "for s in f:\n",
    "    geoid = s[\"properties\"][\"GEOID\"]\n",
    "    #geom = shape(s['geometry'])\n",
    "    geom = s['geometry']\n",
    "    \n",
    "    block_group_geoids.append(geoid)\n",
    "    block_group_geoms.append(geom)\n",
    "f.close()\n",
    "\n",
    "num_block_groups = len(block_group_geoids)\n",
    "\n",
    "print(\"Finished loading %d block group geometries in %0.4f seconds\" % (num_block_groups, time.time() - load_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_block_group_intersections(slr_amount, block_group_geoms, block_group_geoids):\n",
    "\n",
    "    area = []\n",
    "    area_flooded = []\n",
    "\n",
    "    f = rasterio.open(\"data/processed/digital_coast/%dft_masked_union.tif\" % (slr_amount), \"r\")\n",
    "\n",
    "    for i in range(num_block_groups):\n",
    "        if i % 100000 == 0:\n",
    "            print(\"\\t%d/%d\" % (i+1, num_block_groups))\n",
    "        out_image, out_transform = rasterio.mask.mask(f, [block_group_geoms[i]], crop=True, all_touched=False)\n",
    "\n",
    "        num_outside_shape = np.sum(out_image==-1)\n",
    "        num_non_flooded = np.sum(out_image==0)\n",
    "        num_flooded = np.sum(out_image==1)\n",
    "\n",
    "        area.append(num_non_flooded + num_flooded)\n",
    "        area_flooded.append(num_flooded)\n",
    "\n",
    "        assert np.abs(out_image.size - (num_outside_shape+num_non_flooded+num_flooded)) < 1e-5\n",
    "    f.close()\n",
    "    \n",
    "    \n",
    "    f = open(\"data/processed/slr_%dft_bg_intersection.csv\" % (slr_amount),\"w\")\n",
    "    f.write(\"GEOID,Total Area,Area Flooded,Percent Flooded\\n\")\n",
    "    for i in range(num_block_groups):\n",
    "        flooded = 0.0\n",
    "        if area[i] == 0:\n",
    "            #print(\"Block group %s\" % (block_group_geoids[i]))\n",
    "            flooded = 0.0\n",
    "        else:\n",
    "            flooded = area_flooded[i] / float(area[i])\n",
    "        f.write(\"%s,%d,%d,%f\\n\" % (block_group_geoids[i], area[i], area_flooded[i], flooded))\n",
    "    f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n",
      "1\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n",
      "2\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n",
      "3\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n",
      "4\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n",
      "5\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n",
      "6\n",
      "\t1/216330\n",
      "\t100001/216330\n",
      "\t200001/216330\n"
     ]
    }
   ],
   "source": [
    "for i in range(0, 7):\n",
    "    print(i)\n",
    "    process_block_group_intersections(i, block_group_geoms, block_group_geoids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}