CensusAgg.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Script for aggregating quadtree data to census tract level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import packages\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import convert_quadtree_to_latlon as con\n",
    "from shapely.geometry import Polygon\n",
    "from pygeotile.tile import Tile\n",
    "from shapely.geometry import Point\n",
    "import os\n",
    "import os.path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/projects/mpi/shared/Data'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#set working directory\n",
    "os.chdir('/projects/mpi/shared/Data')\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#grab file strings\n",
    "chicago = 'chicago_all.csv'\n",
    "nyc = 'nyc_all.csv'\n",
    "la = 'la_all.csv'\n",
    "sf = 'sf_all.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_list = ['qtid_origin','qtid_destination','OD_counts','num_dates','avg_travel_distance (mile)',\n",
    "            'avg_travel_time (second)','time_period','dominant_mode_1','dominant_mode_2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qt2geom(qt):\n",
    "    bbox = Tile.from_quad_tree(qt).bounds\n",
    "    # min_lng, min_lat, max_lng, max_lat\n",
    "    poly = Polygon.from_bounds(bbox[0][1], bbox[0][0], bbox[1][1], bbox[1][0]) \n",
    "    return poly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qtid_to_geo(filepath,col_list):\n",
    "    '''Function for assigning lat-lon geometry to qtid files'''\n",
    "    #read data into a pandas df\n",
    "    df = pd.read_csv(filepath,header = None,names =col_list,\n",
    "                       dtype={'qtid_origin': object, 'qtid_destination': object})\n",
    "    \n",
    "    #extract qtid for origin and destination and apply polygon mapping function\n",
    "    qcode_list = set(list(df['qtid_origin'])+list(df['qtid_destination']))\n",
    "    qt_poly = list(map(qt2geom,qcode_list))\n",
    "    \n",
    "    #set as dataframes for merging\n",
    "    qcode_df = pd.DataFrame(list(qcode_list),columns=['qtid'])\n",
    "    qt_poly_df = pd.DataFrame(qt_poly,columns=['geometry'])\n",
    "    \n",
    "    #merge polygon and qtid dfs together and merge back with OD data to grab the origin polygon\n",
    "    merged = pd.merge(left=qcode_df,right=qt_poly_df,left_index=True,right_index=True)\n",
    "    \n",
    "    origin = pd.merge(left=df,right=merged,left_on='qtid_origin',right_on='qtid')\n",
    "    origin.drop(['qtid'],axis=1,inplace=True) #clean up\n",
    "    origin[\"trip_id\"] = origin.index + 1\n",
    "    \n",
    "    dest = pd.merge(left=df,right=merged,left_on='qtid_destination',right_on='qtid')\n",
    "    dest.drop(['qtid'],axis=1,inplace=True)\n",
    "    dest[\"trip_id\"] = dest.index + 1\n",
    "    \n",
    "    return origin, dest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generate dataframes with polygons for each city\n",
    "origin_chic, dest_chic = qtid_to_geo(chicago,col_list)\n",
    "origin_nyc, dest_nyc = qtid_to_geo(nyc,col_list)\n",
    "origin_la, dest_la = qtid_to_geo(la,col_list)\n",
    "origin_sf, dest_sf = qtid_to_geo(sf,col_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "#list of census tract files. source: https://www.census.gov/geo/maps-data/data/cbf/cbf_tracts.html\n",
    "\n",
    "census_tract_zips = ['cb_2017_06_tract_500k.zip','cb_2017_09_tract_500k.zip','cb_2017_17_tract_500k.zip',\n",
    "                    'cb_2017_34_tract_500k.zip','cb_2017_36_tract_500k.zip','cb_2017_18_tract_500k.zip',\n",
    "                    'cb_2017_55_tract_500k.zip','cb_2017_26_tract_500k.zip','cb_2017_10_tract_500k.zip',\n",
    "                    'cb_2017_25_tract_500k.zip','cb_2017_42_tract_500k.zip','cb_2017_44_tract_500k.zip']\n",
    "\n",
    "#move and unzip the files\n",
    "for file in census_tract_zips:\n",
    "    if os.path.isfile(os.getcwd()+'/shapefiles/'+file):\n",
    "        os.system('unzip '+os.getcwd()+'/shapefiles/'+file)\n",
    "    else:\n",
    "        os.system('mv '+file+' '+os.getcwd()+'/shapefiles')\n",
    "        os.system('unzip '+os.getcwd()+'/shapefiles/'+file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#list of census tracts to loop over and import and concatenate\n",
    "folder = os.getcwd()+'/shapefiles/'\n",
    "\n",
    "census_tracts = ['cb_2017_06_tract_500k.shp','cb_2017_09_tract_500k.shp','cb_2017_17_tract_500k.shp',\n",
    "                    'cb_2017_34_tract_500k.shp','cb_2017_36_tract_500k.shp','cb_2017_18_tract_500k.shp',\n",
    "                    'cb_2017_55_tract_500k.shp','cb_2017_26_tract_500k.shp','cb_2017_10_tract_500k.shp',\n",
    "                    'cb_2017_25_tract_500k.shp','cb_2017_42_tract_500k.shp','cb_2017_44_tract_500k.shp']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#function for reading in the shapefiles and concatenating them\n",
    "def geo_agg(tract_file_list,path):\n",
    "    gdf_list = []\n",
    "    for i in tract_file_list:\n",
    "        ct = gpd.read_file(path+i)\n",
    "        gdf_list.append(ct)\n",
    "    \n",
    "    all_tracts = pd.concat(gdf_list)\n",
    "    \n",
    "    return all_tracts\n",
    "\n",
    "all_tracts = geo_agg(census_tracts,folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['06', '09', '17', '34', '36', '18', '55', '26', '10', '25', '42',\n",
       "       '44'], dtype=object)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_tracts['STATEFP'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "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>STATEFP</th>\n",
       "      <th>COUNTYFP</th>\n",
       "      <th>TRACTCE</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>06</td>\n",
       "      <td>001</td>\n",
       "      <td>400600</td>\n",
       "      <td>1400000US06001400600</td>\n",
       "      <td>06001400600</td>\n",
       "      <td>4006</td>\n",
       "      <td>CT</td>\n",
       "      <td>297856</td>\n",
       "      <td>0</td>\n",
       "      <td>POLYGON ((-122.26807 37.844136, -122.26514 37....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>06</td>\n",
       "      <td>001</td>\n",
       "      <td>400900</td>\n",
       "      <td>1400000US06001400900</td>\n",
       "      <td>06001400900</td>\n",
       "      <td>4009</td>\n",
       "      <td>CT</td>\n",
       "      <td>420877</td>\n",
       "      <td>0</td>\n",
       "      <td>POLYGON ((-122.285576 37.839778, -122.283186 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>06</td>\n",
       "      <td>001</td>\n",
       "      <td>401400</td>\n",
       "      <td>1400000US06001401400</td>\n",
       "      <td>06001401400</td>\n",
       "      <td>4014</td>\n",
       "      <td>CT</td>\n",
       "      <td>758204</td>\n",
       "      <td>0</td>\n",
       "      <td>POLYGON ((-122.278611 37.826878, -122.268563 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>06</td>\n",
       "      <td>001</td>\n",
       "      <td>403000</td>\n",
       "      <td>1400000US06001403000</td>\n",
       "      <td>06001403000</td>\n",
       "      <td>4030</td>\n",
       "      <td>CT</td>\n",
       "      <td>352394</td>\n",
       "      <td>0</td>\n",
       "      <td>POLYGON ((-122.274757 37.79883299999999, -122....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>06</td>\n",
       "      <td>001</td>\n",
       "      <td>405902</td>\n",
       "      <td>1400000US06001405902</td>\n",
       "      <td>06001405902</td>\n",
       "      <td>4059.02</td>\n",
       "      <td>CT</td>\n",
       "      <td>487280</td>\n",
       "      <td>0</td>\n",
       "      <td>POLYGON ((-122.247175 37.789913, -122.243512 3...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  STATEFP COUNTYFP TRACTCE              AFFGEOID        GEOID     NAME LSAD  \\\n",
       "0      06      001  400600  1400000US06001400600  06001400600     4006   CT   \n",
       "1      06      001  400900  1400000US06001400900  06001400900     4009   CT   \n",
       "2      06      001  401400  1400000US06001401400  06001401400     4014   CT   \n",
       "3      06      001  403000  1400000US06001403000  06001403000     4030   CT   \n",
       "4      06      001  405902  1400000US06001405902  06001405902  4059.02   CT   \n",
       "\n",
       "    ALAND  AWATER                                           geometry  \n",
       "0  297856       0  POLYGON ((-122.26807 37.844136, -122.26514 37....  \n",
       "1  420877       0  POLYGON ((-122.285576 37.839778, -122.283186 3...  \n",
       "2  758204       0  POLYGON ((-122.278611 37.826878, -122.268563 3...  \n",
       "3  352394       0  POLYGON ((-122.274757 37.79883299999999, -122....  \n",
       "4  487280       0  POLYGON ((-122.247175 37.789913, -122.243512 3...  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_tracts.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29732, 10)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_tracts.shape #check total number of census tracts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save aggregated CT file\n",
    "all_tracts.to_file('all_cts.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qtid_ct_join(df1,df2,shp):\n",
    "    '''function to geojoin cts to qtid\n",
    "    input: df= dataframe of OD pair data; shp = tract shapefile\n",
    "    output: qtid level data frame with census tracts tagged'''\n",
    "    #turn dataframes into geodata frames and assign geometry for origin and destination\n",
    "    origin_df = gpd.GeoDataFrame(df1, crs={'init': 'epsg:4269'}, geometry=df1['geometry'])\n",
    "    dest_df = gpd.GeoDataFrame(df2, crs={'init': 'epsg:4269'}, geometry=df2['geometry'])\n",
    "    \n",
    "    #join each origin and destination geo_df to the census tract file\n",
    "    origin_geo = gpd.sjoin(origin_df,shp, how='left', op='intersects')\n",
    "    dest_geo = gpd.sjoin(dest_df,shp, how='left', op='intersects')\n",
    "    \n",
    "    #clean up geodataframes\n",
    "    origin_geo.drop(['index_right','COUNTYFP','AFFGEOID','TRACTCE', 'NAME', 'LSAD', 'ALAND', 'AWATER'],axis=1,inplace=True)\n",
    "    origin_geo.rename(columns={'STATEFP':'origin_st','GEOID':'origin_ct'},inplace=True)\n",
    "    \n",
    "    dest_geo.drop(['index_right','COUNTYFP','AFFGEOID','TRACTCE', 'NAME', 'LSAD', 'ALAND', 'AWATER'],axis=1,inplace=True)\n",
    "    dest_geo.rename(columns={'STATEFP':'dest_st','GEOID':'dest_ct'},inplace=True)\n",
    "    \n",
    "    df_od = pd.merge(left=origin_geo,right=dest_geo,how='left',left_on='trip_id',right_on='trip_id') #merge dataframes\n",
    "    \n",
    "    df_od = df_od[~df_od.index.duplicated(keep='last')] #drop the duplicate indices\n",
    "    \n",
    "    df_od = df_od[['qtid_origin_x','qtid_destination_x','OD_counts_x',\n",
    "                'num_dates_x','time_period_x','origin_st', 'origin_ct','dest_st','dest_ct']] #just keep the columns needed\n",
    "    \n",
    "    df_od.rename({'qtid_origin_x':'qtid_origin','qtid_destination_x':'qtid_dest','OD_counts_x':'od_counts',\n",
    "                'num_dates_x':'num_dates','time_period_x':'time_period'},axis=1,inplace=True)\n",
    "    \n",
    "    return df_od"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert data and output to csv in the cwd at the qtid level\n",
    "new_york = qtid_ct_join(origin_nyc,dest_nyc,all_tracts)\n",
    "new_york.to_csv('ny_qtid_ct.csv')\n",
    "\n",
    "chicago = qtid_ct_join(origin_chic,dest_chic,all_tracts)\n",
    "chicago.to_csv('chic_qtid_ct.csv')\n",
    "\n",
    "la = qtid_ct_join(origin_la,dest_la,all_tracts)\n",
    "la.to_csv('la_qtid_ct.csv')\n",
    "\n",
    "sf = qtid_ct_join(origin_sf,dest_sf,all_tracts)\n",
    "sf.to_csv('sf_qtid_ct.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#aggregate to the CT level\n",
    "chic = chicago.groupby(['origin_ct','dest_ct','time_period','origin_st','dest_st'])[['od_counts','num_dates']].sum()\n",
    "ny = new_york.groupby(['origin_ct','dest_ct','time_period','origin_st','dest_st'])[['od_counts','num_dates']].sum()\n",
    "la = la.groupby(['origin_ct','dest_ct','time_period','origin_st','dest_st'])[['od_counts','num_dates']].sum()\n",
    "sf = sf.groupby(['origin_ct','dest_ct','time_period','origin_st','dest_st'])[['od_counts','num_dates']].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "#output files to cwd\n",
    "chic.to_csv(\"chic_ct.csv\")\n",
    "ny.to_csv(\"ny_ct.csv\")\n",
    "la.to_csv(\"la_ct.csv\")\n",
    "sf.to_csv(\"sf_ct.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "deeplearn",
   "language": "python",
   "name": "deeplearn"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}