Example California water system boundary overlaps.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "681ac8f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Pip install the census_area package from here: https://github.com/datamade/census_area\n",
    "\n",
    "#Example script: https://github.com/datamade/chicago-community-area-ward-demographics/blob/master/scripts/calculate_demographics.py\n",
    "#Documentation: https://census-area.readthedocs.io/en/latest/\n",
    "\n",
    "#Description of the summary file (sf1) variables: \n",
    "#main page: https://www.census.gov/data/datasets/2010/dec/summary-file-1.html\n",
    "#List of tables: https://www2.census.gov/programs-surveys/decennial/2010/technical-documentation/complete-tech-docs/summary-file/sf1.pdf\n",
    "\n",
    "#Description of the ACS 5-year (acs5) variables:\n",
    "#main page: https://www.census.gov/data/developers/data-sets/acs-5year.html\n",
    "#variables: https://api.census.gov/data/2020/acs/acs5/variables.html\n",
    "\n",
    "from census_area import Census\n",
    "\n",
    "import pandas as pd\n",
    "import json\n",
    "\n",
    "#Load up the census data\n",
    "#API KEY REDACTED\n",
    "api_key='######################'\n",
    "c = Census(api_key, year=2020)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b85db228",
   "metadata": {},
   "outputs": [],
   "source": [
    "#BACKGROUND: We are attempting to obtain demographic information for California water systems. \n",
    "#These system boundaries are available in the SABL dataset: \n",
    "#https://gispublic.waterboards.ca.gov/portal/apps/webappviewer/index.html?id=272351aa7db14435989647a86e6d3ad8\n",
    "#Note: The shapefile is only available to download by request.\n",
    "\n",
    "#NECESSARY CONVERSIONS\n",
    "#You only need to do this once to produce the geojson file in the proper format. Comment out afterwards\n",
    "import geopandas as gpd\n",
    "# SABL to geojson\n",
    "sfile = gpd.read_file(\"California_Drinking_Water_System_Area_Boundaries.shp\")\n",
    "sfile.to_file(\"SABL.json\", driver=\"GeoJSON\")\n",
    "\n",
    "#NOTE: I had to manually add the ogr2ogr script to my Library folder in my environment. \n",
    "#Get it here: http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr2ogr.py\n",
    "\n",
    "#To get help on usage: ogr2ogr.Usage()\n",
    "\n",
    "#This is the step to adjust the geojson to the correct projection\n",
    "import os\n",
    "\n",
    "command=\"ogr2ogr SABL_layer.geojson -t_srs \\\"EPSG:4326\\\" SABL.json\"\n",
    "os.system(command)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb195ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Open the geojson\n",
    "with open('SABL_layer.geojson') as infile:\n",
    "    my_shape_geojson = json.load(infile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40f61d9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's say we want only a subset of water systems. Make a list.\n",
    "pwsid_list=['CA0010005','CA4010011','CA3010037','CA1910065','CA2410005']\n",
    "\n",
    "#We're going to create a smaller geojson of just these shapes\n",
    "subset_json={'type':my_shape_geojson['type'],'crs':my_shape_geojson['crs']}\n",
    "subset_list=[]\n",
    "for f in my_shape_geojson['features']:\n",
    "    pwsid=f['properties']['SABL_PWSID']\n",
    "    if pwsid in pwsid_list:\n",
    "        subset_list+=[f]\n",
    "        \n",
    "subset_json['features']=subset_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aaa8f535",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Here we make a list of the ACS fields that we will want to retrieve from the census\n",
    "acs_keys={'B19001_001E':'n_hh_income','B19001_002E':'n_hh_income_lt_10k','B19001_003E':'n_hh_income_10k_15k',\n",
    "         'B19001_004E':'n_hh_income_15k_20k','B19001_005E':'n_hh_income_20k_25k','B19001_006E':'n_hh_income_25k_30k',\n",
    "         'B19001_007E':'n_hh_income_30k_35k','B19001_008E':'n_hh_income_35k_40k','B19001_009E':'n_hh_income_40k_45k',\n",
    "         'B19001_010E':'n_hh_income_45k_50k','B19001_011E':'n_hh_income_50k_60k','B19001_012E':'n_hh_income_60k_75k',\n",
    "         'B19001_013E':'n_hh_income_75k_100k','B19001_014E':'n_hh_income_100k_125k','B19001_015E':'n_hh_income_125k_150k',\n",
    "         'B19001_016E':'n_hh_income_150k_200k','B19001_017E':'n_hh_income_gt_200k','B07412_002E':'n_100pct_pov_lvl',\n",
    "         'B07412_003E':'n_101_149pct_pov_lvl','B07412_004E':'n_150pct_pov_lvl','B25032_002E':'n_hh_own',\n",
    "         'B25032_013E':'n_hh_rent','B08202_001E':'hh_size','B08202_003E':'hh_1worker','B08202_004E':'hh_2worker',\n",
    "         'B08202_005E':'hh_3+worker','B08202_013E':'n_hh_3ppl','B08202_018E':'n_hh_4+ppl','B11011_001E':'hh_type'}\n",
    "\n",
    "vlist=tuple(['NAME']+list(acs_keys.keys()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81ea87ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "#OPTION 1: Store as a Pandas dataframe to export to excel\n",
    "\n",
    "df_store=pd.DataFrame()\n",
    "\n",
    "#Iterate over each water system boundary in the subset geojson\n",
    "for n in range(0,len(subset_list)):\n",
    "    system_area=subset_json['features'][n]\n",
    "    #Generate the overlap\n",
    "    overlap_features=c.acs5.geo_tract(vlist, system_area['geometry'],2020)\n",
    "    #Iterate over each overlapped census tract\n",
    "    for tract_geojson, tract_data, tract_proportion in overlap_features:\n",
    "        #Include the tract properties obtained from the overlap\n",
    "        tract_geojson['properties'].update(tract_data)\n",
    "        #add the proportion of overlap between the census tract and the water system\n",
    "        tract_geojson['properties'].update({'proportion':tract_proportion})\n",
    "        #Add information about the water system ID, name, regulating agency, classification, and total population\n",
    "        for k in ['SABL_PWSID','WATER_SY_1','REGULATING','STATE_CLAS','POPULATION']:\n",
    "            tract_geojson['properties'].update({k:system_area['properties'][k]})\n",
    "        #Append these properties as a new row to the dataframe\n",
    "        df_store=df_store.append(tract_geojson['properties'],ignore_index=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0b78fb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Rename the columns using the acs_keys dictonary\n",
    "df_final=df_store.rename(columns=acs_keys)\n",
    "#Output to excel\n",
    "df_final.to_excel('acs_SABL_overlaps.xlsx',index=False)\n",
    "#Or, output to csv\n",
    "df_final.to_csv('acs_SABL_overlaps.csv',index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fb0b19e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#OPTION 2 (not fully tested): Create a list of geojsons (possibly able to re-export to shapefile?)\n",
    "\n",
    "#List to store the resulting overlap jsons\n",
    "store_jsons=[]\n",
    "for n in range(0,len(subset_list)):\n",
    "    system_area=subset_json['features'][n]\n",
    "    #Generate the overlap\n",
    "    overlap_features=c.acs5.geo_tract(vlist, system_area['geometry'],2020)\n",
    "    #Create the features list for the overlaps\n",
    "    features=[]\n",
    "    for tract_geojson, tract_data, tract_proportion in overlap_features:\n",
    "        tract_geojson['properties'].update(tract_data)\n",
    "        tract_geojson['properties'].update({'proportion':tract_proportion})\n",
    "        for k in ['SABL_PWSID','WATER_SY_1','REGULATING','STATE_CLAS','POPULATION']:\n",
    "            tract_geojson['properties'].update({k:system_area['properties'][k]})\n",
    "        features.append(tract_geojson)\n",
    "    #Create a new geojson and add the features to it\n",
    "    overlap_geojson = {\n",
    "        'type': 'FeatureCollection',\n",
    "        'crs':my_shape_geojson['crs'],\n",
    "        'features': features\n",
    "    }\n",
    "    #Store the new geojson in the list\n",
    "    store_jsons+=[overlap_geojson]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:root] *",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}