Wanying Li - Coding Sample.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f85fb1fc",
   "metadata": {},
   "source": [
    "# README\n",
    "\n",
    "This Python script assists in preprocessing data for the regression analysis to evaluate the efficacy of the NYC bike lane program in enhancing cyclist safety.\n",
    "\n",
    "Data transformation involves the following cleaning, wrangling, and harmonization steps:\n",
    "\n",
    "    1. Clean the Motor Vehicle Collisions Crashes data (i.e., crash data) that will be used in the geospatial analysis (geospatial data wrangling and integration) in ArcGIS.\n",
    "    2. Clean the data generated from ArcGIS after the geospatial analysis (i.e., merged ArcGIS data).\n",
    "    3. Combine multiple years of census data with the data from part 2.\n",
    "    4. Manipulate the dataset into the desired format using 'group by' operations.\n",
    "    5. Calculate the ridability score (requiring the merging and processing of another dataset).\n",
    "\n",
    "It takes in 4 sets of data:\n",
    "\n",
    "    1. Motor Vehicle Collisions Crashes data downloaded from the NYC OpenData Portal as a CSV file (used in Step 1).\n",
    "    2. 4 American Community Survey (ACS) 5-Year Estimates as CSV files: 2014-2018, 2015-2019, 2016-2020, 2016-2021 (used in Step 3).\n",
    "    3. Data outputted from ArcGIS after merging the crash data and the bike lane location as an Excel file (used in Step 2).\n",
    "    4. Data outputted from ArcGIS after the geospatial analysis for the calculation of rideability score as an Excel file (used in Step 5).\n",
    "\n",
    "It outputs 2 datasets:\n",
    "\n",
    "    1. Step 1 outputs a dataset that is used in geospatial analysis before proceeding to Step 2.\n",
    "    2. The final step outputs the clean dataset used for the regression analysis in STATA.\n",
    "\n",
    "Last updated on Oct 12, 2023 by Wanying Li"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dea777a2-3d18-4792-aaee-5fbe5d488c81",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from functools import reduce\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d05a4b4d",
   "metadata": {},
   "source": [
    "# Clean Crash Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c23d6e-fdc2-4ae1-9d19-3394f4d7967f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reading in the raw data\n",
    "# The Motor Vehicle Collisions Crashes data is downloaded from NYC OpenData Portal\n",
    "# https://data.cityofnewyork.us/Public-Safety/Motor-Vehicle-Collisions-Crashes/h9gi-nx95\n",
    "df_crash = pd.read_csv(\"Motor_Vehicle_Collisions_Crashes_data.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d88c4eb-08cb-496e-b4aa-357d490b430b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are the first five lines:\n",
    "df_crash.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4df2f9b5-064c-4f99-b0a7-74be514f6bbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are all the column headers:\n",
    "df_crash.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f44ae80",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are the date types for all columns:\n",
    "df_crash.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1c06ab3-c8c9-45bb-b5c8-3c9171580058",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Are there any duplicate IDs? No. That's good. \n",
    "(df_crash['COLLISION_ID'].value_counts() > 1).any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afe89874-059d-4fb4-9372-c84f9382a6e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What is the shape of the data? # rows and # cols\n",
    "df_crash.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91dfe353-62da-491d-8800-c7ad7b13522f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What was the number of cells not filled in for each question?\n",
    "df_crash.isna().sum()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3cd813ea",
   "metadata": {},
   "source": [
    "## Filter out un-needed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0503e2-2df3-4fcd-934a-561f4f4555ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of the dataframe\n",
    "df_crash1 = df_crash.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7133d74",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the data that doesn't have a latitude and longitude (case 1: field is empty)\n",
    "df_crash2 = df_crash1[(df_crash1['LATITUDE'].notnull()) & (df_crash1['LONGITUDE'].notnull())]\n",
    "print(\"number of dropped data =\", df_crash1.shape[0]-df_crash2.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06150d04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the data that doesn't have a latitude and longitude (case 2: both fields are zero)\n",
    "df_crash3 = df_crash2[(df_crash2['LATITUDE']!=0) & (df_crash2['LONGITUDE']!=0)]\n",
    "print(\"number of dropped data =\", df_crash2.shape[0]-df_crash3.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8df653d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the data that doesn't involve a cyclist accident\n",
    "df_crash4 = df_crash3[(df_crash3['NUMBER OF CYCLIST INJURED']>0) | (df_crash3['NUMBER OF CYCLIST KILLED']>0)]\n",
    "print(\"number of dropped data =\", df_crash3.shape[0]-df_crash4.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d063a080",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the data that was after 2022/01/01\n",
    "df_crash4['CRASH DATE'] = pd.to_datetime(df_crash4['CRASH DATE']).dt.date\n",
    "df_crash5 = df_crash4[(df_crash4['CRASH DATE']<datetime.date(2022,1,1))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13e337b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the data that was before 2017/01/01\n",
    "df_crash6 = df_crash5[(df_crash5['CRASH DATE']>=datetime.date(2017,1,1))]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7bf233d9",
   "metadata": {},
   "source": [
    "## Export data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a13ea6c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the filtered dataframe as csv file\n",
    "df_crash6.to_csv('Motor_Vehicle_Collisions_Crashes_data_clean.csv',index=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a04bbf0f",
   "metadata": {},
   "source": [
    "# Plot Cyclist Injury and Fatality Trend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b48e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "years = list(range(2015,2021+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0030b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute cyclist injury and fatality number for each year\n",
    "cyclist_injuries = []\n",
    "cyclist_fatalities = []\n",
    "\n",
    "for year in years:\n",
    "    cyclist_injuries.append(df_crash1[(pd.to_datetime(df_crash1['CRASH DATE']).dt.year == year) & (df_crash1['NUMBER OF CYCLIST INJURED']>0)].shape[0])\n",
    "    cyclist_fatalities.append(df_crash1[(pd.to_datetime(df_crash1['CRASH DATE']).dt.year == year) & (df_crash1['NUMBER OF CYCLIST KILLED']>0)].shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ba2996e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute all crash-related injury and fatality number for each year\n",
    "all_injuries = []\n",
    "all_fatalities = []\n",
    "\n",
    "for year in years:\n",
    "    all_injuries.append(df_crash1[(pd.to_datetime(df_crash1['CRASH DATE']).dt.year == year) & (df_crash1['NUMBER OF PERSONS INJURED']>0)].shape[0])\n",
    "    all_fatalities.append(df_crash1[(pd.to_datetime(df_crash1['CRASH DATE']).dt.year == year) & (df_crash1['NUMBER OF PERSONS KILLED']>0)].shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fff26290",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute cyclist-to-all-crash accident ratio\n",
    "injury_ratio = [x/y for x,y in zip(cyclist_injuries, all_injuries)]\n",
    "fatality_ratio = [x/y for x,y in zip(cyclist_fatalities, all_fatalities)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebb87869",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot cyclist injury trend\n",
    "x = years[1:]\n",
    "y = cyclist_injuries[1:]\n",
    "plt.plot(x, y, label = \"# of cyclist injuries in NYC\")\n",
    "\n",
    "# Labels\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Number of People')\n",
    "plt.title('Cyclist Injury Trend')\n",
    "plt.legend()\n",
    "\n",
    "# Set the ticks\n",
    "plt.xticks(np.arange(min(x), max(x)+1, 1))\n",
    "\n",
    "# Save the plot\n",
    "plt.savefig('cyclist_injury_trend.jpg', dpi=300, bbox_inches=\"tight\")\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b98d0a10",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot cyclist fatality trend\n",
    "x = years[1:]\n",
    "y = cyclist_fatalities[1:]\n",
    "plt.plot(x, y, label = \"# of cyclist fatalities in NYC\")\n",
    "\n",
    "# Label\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Number of People')\n",
    "plt.title('Cyclist Fatality Trend')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "# Set the ticks\n",
    "plt.xticks(np.arange(min(x), max(x)+1, 1))\n",
    "\n",
    "# Save the plot\n",
    "plt.savefig('cyclist_fatality_trend.jpg', dpi=300, bbox_inches=\"tight\")\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5bc5847",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot cyclist injury ratio trend\n",
    "x = years[1:]\n",
    "y = injury_ratio[1:]\n",
    "plt.plot(x, y, label = \"Ratio of cyclist injury to all crash injury\")\n",
    "\n",
    "# Label\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Percentage')\n",
    "plt.title('Cyclist Injury Ratio Trend')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "# Set the ticks\n",
    "plt.xticks(np.arange(min(x), max(x)+1, 1))\n",
    "\n",
    "# Save the plot\n",
    "plt.savefig('cyclist_injury_ratio_trend.jpg', dpi=300, bbox_inches=\"tight\")\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46be1b54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot cyclist fatality ratio trend\n",
    "x = years[1:]\n",
    "y = fatality_ratio[1:]\n",
    "plt.plot(x, y, label = \"Ratio of cyclist fatality to all crash fatality\")\n",
    "\n",
    "# Label\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Percentage')\n",
    "plt.title('Cyclist Fatality Ratio Trend')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "# Set the ticks\n",
    "plt.xticks(np.arange(min(x), max(x)+1, 1))\n",
    "\n",
    "# Save the plot\n",
    "plt.savefig('cyclist_fatality_ratio_trend.jpg', dpi=300, bbox_inches=\"tight\")\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3cd813ea",
   "metadata": {},
   "source": [
    "# Clean Merged ArcGIS Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c23d6e-fdc2-4ae1-9d19-3394f4d7967f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reading in data\n",
    "df = pd.read_excel(\"Merged_data_from_ArcGIS.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d88c4eb-08cb-496e-b4aa-357d490b430b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are the first five lines\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afe89874-059d-4fb4-9372-c84f9382a6e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What is the shape of the data? # rows and # cols\n",
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4df2f9b5-064c-4f99-b0a7-74be514f6bbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here are all the column names\n",
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0503e2-2df3-4fcd-934a-561f4f4555ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of the dataframe\n",
    "df1 = df.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e419d19a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the GEOID column from float64 (e.g. 3.604700e+10) to int64 (36005000100)\n",
    "df1['CT2020_GEOID'] = df1['CT2020_GEOID'].astype('Int64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "182f0085",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking -- note: there will still be dozens of records that don't have GEOID, those crashes usually are on bridges or turnpike\n",
    "df1[df1['CT2020_GEOID'].isnull()][['CRASH_DATE','CRASH_TIME','BOROUGH','ZIP_CODE','LATITUDE','LONGITUDE', 'ON_STREET_NAME','CROSS_STREET_NAME']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16f52f68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the columns created because of ArcGIS's \"spatial join\" operation (6 join * 3 col/join = 18 columns)\n",
    "ArcGIS_drop_columns=['OBJECTID', 'Join_Count', 'TARGET_FID']\n",
    "\n",
    "# Drop all columns in census tract dataset except 'CT2020_GEOID'\n",
    "census_tract_drop_columns = [\n",
    "    # 'CT2020_GEOID'\n",
    "    'CTLabel', \n",
    "    'BoroCode', \n",
    "    'BoroName',\n",
    "    'CT2020', \n",
    "    'BoroCT2020', \n",
    "    'CDEligibil', \n",
    "    'NTAName', \n",
    "    'NTA2020', \n",
    "    'CDTA2020',\n",
    "    'CDTANAME', \n",
    "    'Shape_Leng'\n",
    "]\n",
    "\n",
    "# Drop unuseful columns from crash datasets\n",
    "crash_drop_columns = [\n",
    "\t# 'CRASH_DATE',\n",
    "    'CRASH_TIME',\n",
    "    'BOROUGH',\n",
    "    # 'ZIP_CODE',\n",
    "    'LATITUDE',\n",
    "    'LONGITUDE',\n",
    "    'LOCATION',\n",
    "    'ON_STREET_NAME',\n",
    "    'CROSS_STREET_NAME',\n",
    "    'OFF_STREET_NAME',\n",
    "    # 'NUMBER_OF_PERSONS_INJURED',\n",
    "    # 'NUMBER_OF_PERSONS_KILLED',\n",
    "    'NUMBER_OF_PEDESTRIANS_INJURED',\n",
    "    'NUMBER_OF_PEDESTRIANS_KILLED',\n",
    "    # 'NUMBER_OF_CYCLIST_INJURED',\n",
    "    # 'NUMBER_OF_CYCLIST_KILLED',\n",
    "    'NUMBER_OF_MOTORIST_INJURED',\n",
    "    'NUMBER_OF_MOTORIST_KILLED',\n",
    "    'CONTRIBUTING_FACTOR_VEHICLE_1',\n",
    "    'CONTRIBUTING_FACTOR_VEHICLE_2',\n",
    "    'CONTRIBUTING_FACTOR_VEHICLE_3',\n",
    "    'CONTRIBUTING_FACTOR_VEHICLE_4',\n",
    "    'CONTRIBUTING_FACTOR_VEHICLE_5',\n",
    "    'COLLISION_ID',\n",
    "    'VEHICLE_TYPE_CODE_1',\n",
    "    'VEHICLE_TYPE_CODE_2',\n",
    "    'VEHICLE_TYPE_CODE_3',\n",
    "    'VEHICLE_TYPE_CODE_4',\n",
    "    'VEHICLE_TYPE_CODE_5',\n",
    "]\n",
    "\n",
    "df2 = df1.drop(columns=ArcGIS_drop_columns+census_tract_drop_columns+crash_drop_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcf265cc",
   "metadata": {},
   "source": [
    "### Create time-related variables from CRASH DATE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f074ccf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df3 = df2.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d76468f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a numerical variable column to indicate the year of crash accident\n",
    "df3['CRASH_YEAR'] = pd.to_datetime(df3['CRASH_DATE']).dt.year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a numerical variable column to indicate the month of crash accident\n",
    "df3['CRASH_MONTH'] = pd.to_datetime(df3['CRASH_DATE']).dt.month"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68db08fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a numerical variable column to indicate the year+month of crash accident\n",
    "df3['CRASH_YEAR-MONTH'] = pd.to_datetime(df3['CRASH_DATE']).dt.to_period('m')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "730307b9",
   "metadata": {},
   "source": [
    "# Merge All Census Data to the Merged ArcGIS Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43c0eac",
   "metadata": {},
   "source": [
    "### Merge all census data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "391b3dc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import census data\n",
    "df_acs2018 = pd.read_csv(\"ACS_data_2014-2018.txt\",sep='\\t')\n",
    "df_acs2019 = pd.read_csv(\"ACS_data_2015-2019.txt\",sep='\\t')\n",
    "df_acs2020 = pd.read_csv(\"ACS_data_2016-2020.txt\",sep='\\t')\n",
    "df_acs2021 = pd.read_csv(\"ACS_data_2017-2021.txt\",sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d919c99",
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_acs_data(dataframe, year):\n",
    "    '''\n",
    "    This function takes in one single ACS dataframe and its corresponding year, cleans it, and outputs the dataframe.\n",
    "    '''\n",
    "    df = dataframe.copy()\n",
    "\n",
    "    # Convert the data type for later join operation & create as a new column under the same name as census tract data\n",
    "    # Census tract uses CT2020_GEOID, ACS uses Geo_FIPS\n",
    "    df['CT2020_GEOID'] = df['Geo_FIPS'].astype('Int64')\n",
    "\n",
    "    # Keep the borough column by renaming it before dropping all the others\n",
    "    df['borough'] = df['Geo_COUNTY']\n",
    "\n",
    "    # Drop all columns that start with \"Geo_\" (ie. geo data, non-demographic data)\n",
    "    df = df.loc[:,~df.columns.str.startswith('Geo_')]\n",
    "\n",
    "    # Rename columns\n",
    "    # A00001_001:     Total Population\n",
    "    # A00002_002:     Population Density (Per Sq. Mile)\n",
    "    # B12001_001:     Population 25 Years and Over\n",
    "    # B12001_002:     Population 25 Years and Over: Less than High School\n",
    "    # B12001_003:     Population 25 Years and Over: High School Diploma\n",
    "    # B12001_004:     Population 25 Years and Over: Bachelor's Degree or Better\n",
    "    # A14006_001:     Median Household Income (In 2021 Inflation Adjusted Dollars) [Dollars adjusted for inflation to match value in 2021]\n",
    "    # A09005_001:     Workers 16 Years and Over:\n",
    "    # A09005_002:     Workers 16 Years and Over: Car, Truck, or Van\n",
    "    # A09005_003:     Workers 16 Years and Over: Public Transportation (Includes Taxicab)\n",
    "    # A09005_005:     Workers 16 Years and Over: Bicycle\n",
    "    # A09003_001:     Average Commute to Work (In Min)\n",
    "    df = df.rename({\n",
    "        'SE_A00001_001':'ttl_pop',\n",
    "        'SE_A00002_002':'pop_density_per_sq_mil',\n",
    "        'SE_B12001_001':'pop_25_yr_over',\n",
    "        'SE_B12001_002':'educ_less_hs',\n",
    "        'SE_B12001_003':'educ_hs',\n",
    "        'SE_B12001_004':'educ_bs_over',\n",
    "        'SE_A14006_001':'median_household_inc',\n",
    "        'SE_A09005_001':'workers_16_yr_over',\n",
    "        'SE_A09005_002':'tranport_mean_car',\n",
    "        'SE_A09005_003':'tranport_mean_public',\n",
    "        'SE_A09005_005':'tranport_mean_bike',\n",
    "        'SE_A09003_001':'avg_commmute_to_work_min'\n",
    "        }, axis='columns')\n",
    "\n",
    "    # Drop all remaining columns that start with \"SE\"\n",
    "    df = df.loc[:,~df.columns.str.startswith('SE_')]\n",
    "\n",
    "    # Compute the \"population over 25 years and over for education\"\n",
    "    df['educ_less_hs_pct'] = df['educ_less_hs']/df['pop_25_yr_over']\n",
    "    df['educ_hs_pct'] = df['educ_hs']/df['pop_25_yr_over']\n",
    "    df['educ_bs_over_pct'] = df['educ_bs_over']/df['pop_25_yr_over']\n",
    "\n",
    "    # Compute the \"workers over 16 years and over for tranportation mean\"\n",
    "    df['tranport_mean_car_pct'] = df['tranport_mean_car']/df['workers_16_yr_over']\n",
    "    df['tranport_mean_public_pct'] = df['tranport_mean_public']/df['workers_16_yr_over']\n",
    "    df['tranport_mean_bike_pct'] = df['tranport_mean_bike']/df['workers_16_yr_over']\n",
    "\n",
    "    # Drop the columns after we finished the computation\n",
    "    drop_columns = [\n",
    "        'pop_25_yr_over',\n",
    "        'educ_less_hs',\n",
    "        'educ_hs',\n",
    "        'educ_bs_over',\n",
    "        'workers_16_yr_over',\n",
    "        'tranport_mean_car',\n",
    "        'tranport_mean_public',\n",
    "        'tranport_mean_bike',\n",
    "        'avg_commmute_to_work_min'\n",
    "        ]\n",
    "    df = df.drop(columns=drop_columns)\n",
    "\n",
    "    # Add year to column name\n",
    "    for col in df.columns:\n",
    "        if col != 'CT2020_GEOID':\n",
    "            df.rename({col:'ACS'+str(year)+'_'+col}, axis='columns', inplace=True)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf093007",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Process and create new dataframe for ACS data\n",
    "df_acs2018_processed = process_acs_data(df_acs2018, 2018)\n",
    "df_acs2019_processed = process_acs_data(df_acs2019, 2019)\n",
    "df_acs2020_processed = process_acs_data(df_acs2020, 2020)\n",
    "df_acs2021_processed = process_acs_data(df_acs2021, 2021)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81e408b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge data by appending the ACS info to the main df based on GEOID\n",
    "dfs_to_merge = [df3, df_acs2018_processed, df_acs2019_processed, df_acs2020_processed, df_acs2021_processed]\n",
    "df4 = reduce(lambda left, right: pd.merge(left, right, how='inner', on='CT2020_GEOID'), dfs_to_merge)\n",
    "df4.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ab65882",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine ACS data depending on the year\n",
    "def combine_ACS_data(row, variable_name):\n",
    "    '''\n",
    "    This function takes in a default argument for the apply function.\n",
    "    The function will combine 4 years of ACS data into 1 and show the data that corresponds to the crash year. \n",
    "    For example, if the crash happens in year 2021, ACS_ttl_pop will be the data from ACS2021_ttl_pop.\n",
    "    '''\n",
    "    if row['CRASH_YEAR'] == 2018:\n",
    "        return row['ACS2018_'+variable_name]\n",
    "    elif row['CRASH_YEAR'] == 2019:\n",
    "        return row['ACS2019_'+variable_name]\n",
    "    elif row['CRASH_YEAR'] == 2020:\n",
    "        return row['ACS2020_'+variable_name]\n",
    "    elif row['CRASH_YEAR'] == 2021:\n",
    "        return row['ACS2021_'+variable_name]\n",
    "\n",
    "df4['ACS_ttl_pop'] = df4.apply(combine_ACS_data, variable_name='ttl_pop', axis=1)\n",
    "df4['ACS_pop_density_per_sq_mil'] = df4.apply(combine_ACS_data, variable_name='pop_density_per_sq_mil', axis=1)\n",
    "df4['ACS_pop_25_yr_over'] = df4.apply(combine_ACS_data, variable_name='pop_25_yr_over', axis=1)\n",
    "df4['ACS_workers_16_yr_over'] = df4.apply(combine_ACS_data, variable_name='workers_16_yr_over', axis=1)\n",
    "df4['ACS_educ_bs_over'] = df4.apply(combine_ACS_data, variable_name='educ_bs_over', axis=1)\n",
    "df4['ACS_median_household_inc'] = df4.apply(combine_ACS_data, variable_name='median_household_inc', axis=1)\n",
    "df4['ACS_tranport_mean_car'] = df4.apply(combine_ACS_data, variable_name='tranport_mean_car', axis=1)\n",
    "df4['ACS_tranport_mean_public'] = df4.apply(combine_ACS_data, variable_name='tranport_mean_public', axis=1)\n",
    "df4['ACS_tranport_mean_bike'] = df4.apply(combine_ACS_data, variable_name='tranport_mean_bike', axis=1)\n",
    "df4['ACS_borough'] = df4.apply(combine_ACS_data, variable_name='borough', axis=1)\n",
    "\n",
    "# Checking output\n",
    "df4[df4['CRASH_YEAR'].notnull()][['CRASH_YEAR','ACS_ttl_pop','ACS_pop_density_per_sq_mil','ACS_tranport_mean_car']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81432db8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop unneeded columns\n",
    "drop_columns = [\n",
    "    'ACS2018_ttl_pop', \n",
    "    'ACS2018_pop_density_per_sq_mil', \n",
    "    'ACS2018_pop_25_yr_over',\n",
    "    'ACS2018_workers_16_yr_over',\n",
    "    'ACS2018_educ_bs_over', \n",
    "    'ACS2018_median_household_inc', \n",
    "    'ACS2018_tranport_mean_car', \n",
    "    'ACS2018_tranport_mean_public', \n",
    "    'ACS2018_tranport_mean_bike', \n",
    "    'ACS2018_borough', \n",
    "    'ACS2019_ttl_pop', \n",
    "    'ACS2019_pop_density_per_sq_mil', \n",
    "    'ACS2019_pop_25_yr_over',\n",
    "    'ACS2019_workers_16_yr_over',\n",
    "    'ACS2019_educ_bs_over', \n",
    "    'ACS2019_median_household_inc', \n",
    "    'ACS2019_tranport_mean_car', \n",
    "    'ACS2019_tranport_mean_public', \n",
    "    'ACS2019_tranport_mean_bike', \n",
    "    'ACS2019_borough', \n",
    "    'ACS2020_ttl_pop', \n",
    "    'ACS2020_pop_density_per_sq_mil', \n",
    "    'ACS2020_pop_25_yr_over',\n",
    "    'ACS2020_workers_16_yr_over',\n",
    "    'ACS2020_educ_bs_over', \n",
    "    'ACS2020_median_household_inc', \n",
    "    'ACS2020_tranport_mean_car', \n",
    "    'ACS2020_tranport_mean_public', \n",
    "    'ACS2020_tranport_mean_bike', \n",
    "    'ACS2020_borough', \n",
    "    'ACS2021_ttl_pop', \n",
    "    'ACS2021_pop_density_per_sq_mil', \n",
    "    'ACS2021_pop_25_yr_over',\n",
    "    'ACS2021_workers_16_yr_over',\n",
    "    'ACS2021_educ_bs_over', \n",
    "    'ACS2021_median_household_inc', \n",
    "    'ACS2021_tranport_mean_car', \n",
    "    'ACS2021_tranport_mean_public', \n",
    "    'ACS2021_tranport_mean_bike', \n",
    "    'ACS2021_borough']\n",
    "df5 = df4.drop(columns=drop_columns)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "608f3a30",
   "metadata": {},
   "source": [
    "# Group by Census Tract and Crash Date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d01f9ba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df6 = df5.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7b18899",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Group by (year-month, census tract), then sum by death and injury number, use max for the ACS numbers\n",
    "\n",
    "agg_func_math = {\n",
    "    'NUMBER_OF_PERSONS_INJURED': ['sum'],\n",
    "    'NUMBER_OF_PERSONS_KILLED': ['sum'],\n",
    "    'NUMBER_OF_CYCLIST_INJURED': ['sum'],\n",
    "    'NUMBER_OF_CYCLIST_KILLED': ['sum'],\n",
    "    'ACS_ttl_pop': ['median'],\n",
    "    'ACS_pop_density_per_sq_mil': ['median'],\n",
    "    'ACS_pop_25_yr_over': ['median'],\n",
    "    'ACS_workers_16_yr_over': ['median'],\n",
    "    'ACS_educ_bs_over': ['median'],\n",
    "    'ACS_median_household_inc': ['median'],\n",
    "    'ACS_tranport_mean_car': ['median'],\n",
    "    'ACS_tranport_mean_public': ['median'],\n",
    "    'ACS_tranport_mean_bike': ['median'],\n",
    "    'ACS_borough': ['median']\n",
    "}\n",
    "df7 = df6.groupby(['CRASH_YEAR-MONTH','CT2020_GEOID'], as_index=False).agg(agg_func_math)\n",
    "df7.columns = df7.columns.droplevel(-1)\n",
    "df7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a056362f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a numerical variable column to indicate the year of crash accident\n",
    "df7['CRASH_YEAR'] = pd.to_datetime(df7['CRASH_YEAR-MONTH'].astype('datetime64[ns]')).dt.year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0a4be84",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a numerical variable column to indicate the month of crash accident\n",
    "df7['CRASH_MONTH'] = pd.to_datetime(df7['CRASH_YEAR-MONTH'].astype('datetime64[ns]')).dt.month"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c4156e77",
   "metadata": {},
   "source": [
    "# Compute Ridability Score (Per Census Tract)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5e768f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import length data\n",
    "df_road_length = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='road_by_ct_Statistics')\n",
    "df_bikelane2021_class1 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2021_class1_Statistics')\n",
    "df_bikelane2021_class2 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2021_class2_Statistics')\n",
    "df_bikelane2021_class3 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2021_class3_Statistics')\n",
    "df_bikelane2020_class1 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2020_class1_Statistics')\n",
    "df_bikelane2020_class2 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2020_class2_Statistics')\n",
    "df_bikelane2020_class3 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2020_class3_Statistics')\n",
    "df_bikelane2019_class1 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2019_class1_Statistics')\n",
    "df_bikelane2019_class2 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2019_class2_Statistics')\n",
    "df_bikelane2019_class3 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2019_class3_Statistics')\n",
    "df_bikelane2018_class1 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2018_class1_Statistics')\n",
    "df_bikelane2018_class2 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2018_class2_Statistics')\n",
    "df_bikelane2018_class3 = pd.read_excel(\"Length_Data_for_Ridability_Score_Calculation.xlsx\", sheet_name='bikelane2018_class3_Statistics')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e2523a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_road_length.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a46ed0d2",
   "metadata": {},
   "source": [
    "## Clean, process, merge length data from ArcGIS output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dcf919c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop columns\n",
    "df_road_length = df_road_length.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "\n",
    "df_bikelane2021_class1 = df_bikelane2021_class1.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2021_class2 = df_bikelane2021_class2.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2021_class3 = df_bikelane2021_class3.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "\n",
    "df_bikelane2020_class1 = df_bikelane2020_class1.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2020_class2 = df_bikelane2020_class2.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2020_class3 = df_bikelane2020_class3.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "\n",
    "df_bikelane2019_class1 = df_bikelane2019_class1.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2019_class2 = df_bikelane2019_class2.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2019_class3 = df_bikelane2019_class3.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "\n",
    "df_bikelane2018_class1 = df_bikelane2018_class1.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2018_class2 = df_bikelane2018_class2.drop(columns=['OBJECTID', 'FREQUENCY'])\n",
    "df_bikelane2018_class3 = df_bikelane2018_class3.drop(columns=['OBJECTID', 'FREQUENCY'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3018c1f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename the length column (so that all have the same name)\n",
    "df_road_length = df_road_length.rename({'SUM_road_length':'road_length'}, axis='columns')\n",
    "\n",
    "df_bikelane2021_class1 = df_bikelane2021_class1.rename({'SUM_bikelane2021_class1_length':'class1_length'}, axis='columns')\n",
    "df_bikelane2021_class2 = df_bikelane2021_class2.rename({'SUM_bikelane2021_class2_length':'class2_length'}, axis='columns')\n",
    "df_bikelane2021_class3 = df_bikelane2021_class3.rename({'SUM_bikelane2021_class3_length':'class3_length'}, axis='columns')\n",
    "\n",
    "df_bikelane2020_class1 = df_bikelane2020_class1.rename({'SUM_bikelane2020_class1_length':'class1_length'}, axis='columns')\n",
    "df_bikelane2020_class2 = df_bikelane2020_class2.rename({'SUM_bikelane2020_class2_length':'class2_length'}, axis='columns')\n",
    "df_bikelane2020_class3 = df_bikelane2020_class3.rename({'SUM_bikelane2020_class3_length':'class3_length'}, axis='columns')\n",
    "\n",
    "df_bikelane2019_class1 = df_bikelane2019_class1.rename({'SUM_bikelane2019_class1_length':'class1_length'}, axis='columns')\n",
    "df_bikelane2019_class2 = df_bikelane2019_class2.rename({'SUM_bikelane2019_class2_length':'class2_length'}, axis='columns')\n",
    "df_bikelane2019_class3 = df_bikelane2019_class3.rename({'SUM_bikelane2019_class3_length':'class3_length'}, axis='columns')\n",
    "\n",
    "df_bikelane2018_class1 = df_bikelane2018_class1.rename({'SUM_bikelane2018_class1_length':'class1_length'}, axis='columns')\n",
    "df_bikelane2018_class2 = df_bikelane2018_class2.rename({'SUM_bikelane2018_class2_length':'class2_length'}, axis='columns')\n",
    "df_bikelane2018_class3 = df_bikelane2018_class3.rename({'SUM_bikelane2018_class3_length':'class3_length'}, axis='columns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05dc5f97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert from string to interger\n",
    "df_bikelane2021_class1['CT2020_GEOID'] = df_bikelane2021_class1['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2021_class2['CT2020_GEOID'] = df_bikelane2021_class2['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2021_class3['CT2020_GEOID'] = df_bikelane2021_class3['CT2020_GEOID'].astype('Int64')\n",
    "\n",
    "df_bikelane2020_class1['CT2020_GEOID'] = df_bikelane2020_class1['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2020_class2['CT2020_GEOID'] = df_bikelane2020_class2['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2020_class3['CT2020_GEOID'] = df_bikelane2020_class3['CT2020_GEOID'].astype('Int64')\n",
    "\n",
    "df_bikelane2019_class1['CT2020_GEOID'] = df_bikelane2019_class1['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2019_class2['CT2020_GEOID'] = df_bikelane2019_class2['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2019_class3['CT2020_GEOID'] = df_bikelane2019_class3['CT2020_GEOID'].astype('Int64')\n",
    "\n",
    "df_bikelane2018_class1['CT2020_GEOID'] = df_bikelane2018_class1['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2018_class2['CT2020_GEOID'] = df_bikelane2018_class2['CT2020_GEOID'].astype('Int64')\n",
    "df_bikelane2018_class3['CT2020_GEOID'] = df_bikelane2018_class3['CT2020_GEOID'].astype('Int64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32b1abf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make copies of the road_length column for later use (different years)\n",
    "df_road_length2021 = df_road_length.copy()\n",
    "df_road_length2020 = df_road_length.copy()\n",
    "df_road_length2019 = df_road_length.copy()\n",
    "df_road_length2018 = df_road_length.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9350ddca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a year column\n",
    "df_road_length2021['CRASH_YEAR'] = 2021\n",
    "df_road_length2020['CRASH_YEAR'] = 2020\n",
    "df_road_length2019['CRASH_YEAR'] = 2019\n",
    "df_road_length2018['CRASH_YEAR'] = 2018\n",
    "\n",
    "df_bikelane2021_class1['CRASH_YEAR'] = 2021\n",
    "df_bikelane2021_class2['CRASH_YEAR'] = 2021\n",
    "df_bikelane2021_class3['CRASH_YEAR'] = 2021\n",
    "\n",
    "df_bikelane2020_class1['CRASH_YEAR'] = 2020\n",
    "df_bikelane2020_class2['CRASH_YEAR'] = 2020\n",
    "df_bikelane2020_class3['CRASH_YEAR'] = 2020\n",
    "\n",
    "df_bikelane2019_class1['CRASH_YEAR'] = 2019\n",
    "df_bikelane2019_class2['CRASH_YEAR'] = 2019\n",
    "df_bikelane2019_class3['CRASH_YEAR'] = 2019\n",
    "\n",
    "df_bikelane2018_class1['CRASH_YEAR'] = 2018\n",
    "df_bikelane2018_class2['CRASH_YEAR'] = 2018\n",
    "df_bikelane2018_class3['CRASH_YEAR'] = 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99879603",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each year, merge length data of road and bike lane\n",
    "dfs_length2021 = [df_road_length2021,df_bikelane2021_class1,df_bikelane2021_class2,df_bikelane2021_class3]\n",
    "df_length2021_merge = reduce(lambda left, right: pd.merge(left, right, how='outer', on=['CT2020_GEOID','CRASH_YEAR']), dfs_length2021)\n",
    "\n",
    "dfs_length2020 = [df_road_length2020,df_bikelane2020_class1,df_bikelane2020_class2,df_bikelane2020_class3]\n",
    "df_length2020_merge = reduce(lambda left, right: pd.merge(left, right, how='outer', on=['CT2020_GEOID','CRASH_YEAR']), dfs_length2020)\n",
    "\n",
    "dfs_length2019 = [df_road_length2019,df_bikelane2019_class1,df_bikelane2019_class2,df_bikelane2019_class3]\n",
    "df_length2019_merge = reduce(lambda left, right: pd.merge(left, right, how='outer', on=['CT2020_GEOID','CRASH_YEAR']), dfs_length2019)\n",
    "\n",
    "dfs_length2018 = [df_road_length2018,df_bikelane2018_class1,df_bikelane2018_class2,df_bikelane2018_class3]\n",
    "df_length2018_merge = reduce(lambda left, right: pd.merge(left, right, how='outer', on=['CT2020_GEOID','CRASH_YEAR']), dfs_length2018)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a161774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Concatenate the data from different years in the vertical direction\n",
    "dfs_to_concat = [df_length2021_merge,\n",
    "    df_length2020_merge,\n",
    "    df_length2019_merge,\n",
    "    df_length2018_merge\n",
    "]\n",
    "dfs_length_concat = pd.concat(dfs_to_concat, axis=0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "df140272",
   "metadata": {},
   "source": [
    "## Merge bike lane & road length data with the crash & census data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2339d61b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df8 = df7.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7153a3df",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_non2017 = df8[(df8['CRASH_YEAR']!=2017)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06fb20e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merge the road & bike lane length data to the crash data by CT2020_ID and year\n",
    "df9 = pd.merge(df_non2017,dfs_length_concat, on=['CT2020_GEOID','CRASH_YEAR'], how ='left')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "82fc1e8b",
   "metadata": {},
   "source": [
    "## Compute score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7747ffba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replace NaN by 0\n",
    "df9['class1_length'] = df9['class1_length'].fillna(0)\n",
    "df9['class2_length'] = df9['class2_length'].fillna(0)\n",
    "df9['class3_length'] = df9['class3_length'].fillna(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1640f567",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the percent of bike lane and non-bike lane out of all roads\n",
    "df9['class1_percent'] = df9['class1_length']/df9['road_length']\n",
    "df9['class2_percent'] = df9['class2_length']/df9['road_length']\n",
    "df9['class3_percent'] = df9['class3_length']/df9['road_length']\n",
    "df9['no_bikelane_percent'] = (df9['road_length']-df9['class1_length']-df9['class2_length']-df9['class3_length'])/df9['road_length']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd6aa485",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute ridability score\n",
    "df9['score1'] = df9['no_bikelane_percent']*(-1) + df9['class1_percent']*2 + df9['class2_percent']*1 + df9['class3_percent']*0\n",
    "df9['score2'] = df9['no_bikelane_percent']*0 + df9['class1_percent']*3 + df9['class2_percent']*2 + df9['class3_percent']*1\n",
    "df9['score3'] = df9['no_bikelane_percent']*1 + df9['class1_percent']*1000 + df9['class2_percent']*100 + df9['class3_percent']*10"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "59abd89b",
   "metadata": {},
   "source": [
    "# Data for Thesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52fa011f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df10 = df9.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc383fee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop columns\n",
    "drop_columns = [\n",
    "   'CRASH_YEAR-MONTH',\n",
    "   # 'CT2020_GEOID',\n",
    "   'NUMBER_OF_PERSONS_INJURED',\n",
    "   'NUMBER_OF_PERSONS_KILLED',\n",
    "   # 'NUMBER_OF_CYCLIST_INJURED',\n",
    "   # 'NUMBER_OF_CYCLIST_KILLED',\n",
    "   # 'ACS_ttl_pop',\n",
    "   # 'ACS_pop_density_per_sq_mil',\n",
    "   # 'ACS_pop_25_yr_over',\n",
    "   # 'ACS_workers_16_yr_over',\n",
    "   # 'ACS_educ_bs_over',\n",
    "   # 'ACS_median_household_inc',\n",
    "   # 'ACS_tranport_mean_car',\n",
    "   # 'ACS_tranport_mean_public',\n",
    "   # 'ACS_tranport_mean_bike',\n",
    "   # 'ACS_borough',\n",
    "   # 'CRASH_YEAR',\n",
    "   # 'CRASH_MONTH',\n",
    "   # 'road_length',\n",
    "   # 'class1_length',\n",
    "   # 'class2_length',\n",
    "   # 'class3_length',\n",
    "   'class1_percent',\n",
    "   'class2_percent',\n",
    "   'class3_percent',\n",
    "   'no_bikelane_percent',\n",
    "   # 'score1',\n",
    "   'score2',\n",
    "   'score3'\n",
    "]\n",
    "df10 = df10.drop(columns=drop_columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f145a5d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename columns\n",
    "df10 = df10.rename({\n",
    "    'CT2020_GEOID':'census_tract_id',\n",
    "    'score1':'ridability_score',\n",
    "    'NUMBER_OF_CYCLIST_INJURED':'cyclist_injuries',\n",
    "    'NUMBER_OF_CYCLIST_KILLED':'cyclist_death',\n",
    "    'ACS_ttl_pop':'ttl_pop',\n",
    "    'ACS_pop_25_yr_over': 'pop_25yr_educ',\n",
    "    'ACS_workers_16_yr_over': 'worker_16yr_transport',\n",
    "    'ACS_pop_density_per_sq_mil':'pop_density',\n",
    "    'ACS_median_household_inc':'income',\n",
    "    'ACS_educ_bs_over':'educ',\n",
    "    'ACS_tranport_mean_car':'car',\n",
    "    'ACS_tranport_mean_public':'public_transportation',\n",
    "    'ACS_tranport_mean_bike':'bike',\n",
    "    'ACS_borough':'borough',\n",
    "    'CRASH_YEAR':'crash_year',\n",
    "    'CRASH_MONTH':'crash_month',\n",
    "    }, axis='columns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61b89b37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop data entry that has missing values (which cannot be used in regression in STATA)\n",
    "df10 = df10.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf233d9",
   "metadata": {},
   "source": [
    "# Export data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e4b565b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export thesis data\n",
    "df10.to_csv('Final Data/final_data.csv',index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.4 64-bit",
   "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.10.4 (tags/v3.10.4:9d38120, Mar 23 2022, 23:13:41) [MSC v.1929 64 bit (AMD64)]"
  },
  "vscode": {
   "interpreter": {
    "hash": "26de051ba29f2982a8de78e945f0abaf191376122a1563185a90213a26c5da77"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}