Assignment 7: Predictive Modeling of Housing Prices in Philadelphia¶

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import geopandas as gpd
import hvplot.pandas
import esri2gpd
import carto2gpd
import seaborn as sns
import requests
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
from sklearn.neighbors import NearestNeighbors
import osmnx as ox
from mpl_toolkits.axes_grid1 import make_axes_locatable
import hvplot.pandas
import holoviews as hv
hv.extension("bokeh")

np.random.seed(42)

%matplotlib inline

pd.options.display.max_columns = 999
pd.options.display.max_rows = 999

carto_url = "https://phl.carto.com/api/v2/sql"

table_name = "opa_properties_public"

sale = "sale_date >= '2021-01-01' AND sale_date <= '2021-12-31'" 
single = " and category_code_description IN ('SINGLE FAMILY','Single Family')"
where = sale + single

salesRaw = carto2gpd.get(carto_url, table_name, where=where)

salesRaw

join_hoods = gpd.sjoin(salesRaw, neighborhoods, how='left', predicate='intersects')

join_hoods = join_hoods.drop(labels=['index_right'], axis=1)

joined = gpd.sjoin(join_hoods, tracts, how='left', predicate='intersects')

joined

sales

graffiti_removal = req_date + " AND service_name = 'Graffiti Removal'"

graffiti = carto2gpd.get(carto_url, table_311, where=graffiti_removal)

graffiti = graffiti.loc[graffiti.geometry.notnull()]

graffiti_3857 = graffiti.to_crs(epsg=3857)

graffitiXY = xy_from_geom(graffiti_3857)

graffiti_nbrs = NearestNeighbors(n_neighbors=5)

graffiti_nbrs.fit(graffitiXY)

grafDists, grafIndices = graffiti_nbrs.kneighbors(salesXY) 

avgGrafDist = grafDists.mean(axis=1)

avgGrafDist[avgGrafDist==0] = 1e-5

sales['logDistGraffiti'] = np.log10(avgGrafDist)

#Fire Safety Complaints
fire_safety = req_date + " AND service_name = 'Fire Safety Complaint'"

fire = carto2gpd.get(carto_url, table_311, where=fire_safety)

fire = fire.loc[fire.geometry.notnull()]

fire_3857 = fire.to_crs(epsg=3857)

fireXY = xy_from_geom(fire_3857)

fire_nbrs = NearestNeighbors(n_neighbors=3)

fire_nbrs.fit(fireXY)

fireDists, fireIndices = fire_nbrs.kneighbors(salesXY) 

avgFireDist = fireDists.mean(axis=1)

avgFireDist[avgFireDist==0] = 1e-5

sales['logDistFire'] = np.log10(avgFireDist)

incidents_table = "incidents_part1_part2"

dispatch = "dispatch_date >= '2021-01-01' AND dispatch_date < '2022-01-01'"
agg = " AND Text_General_Code IN ('Aggravated Assault No Firearm', 'Aggravated Assault Firearm')"
where_assaults = dispatch + agg

assaults = carto2gpd.get(carto_url, incidents_table, where=where_assaults)

assaults = assaults.loc[assaults.geometry.notnull()]
    
assaultsXY = xy_from_geom(assaults.to_crs(epsg=3857))

assaults_nbrs = NearestNeighbors(n_neighbors=5)

assaults_nbrs.fit(assaultsXY)

assaultDists, _ = assaults_nbrs.kneighbors(salesXY)

# Add the new feature
sales['logDistAssaults'] = np.log10(assaultDists.mean(axis=1))

#Three Nearest Schools
schools_url = "https://services.arcgis.com/fLeGjb7u4uXqeF9q/arcgis/rest/services/Schools/FeatureServer/0"
schools = esri2gpd.get(schools_url)

schoolsXY = xy_from_geom(schools.to_crs(epsg=3857))

schools_nbrs = NearestNeighbors(n_neighbors=3)

schools_nbrs.fit(schoolsXY)

schoolsDists, _ = schools_nbrs.kneighbors(salesXY)

sales["logDistSchools"] = np.log10(schoolsDists.mean(axis=1))

#SEPTA Subway Lines
subway = ox.geometries_from_polygon(city_limits_outline, tags={"station": "subway"})

subway = subway.to_crs(epsg=3857)

subwayXY = xy_from_geom(subway.to_crs(epsg=3857))

subway_nbrs = NearestNeighbors(n_neighbors=1)

subway_nbrs.fit(subwayXY)

subwayDists, subwayIndices = subway_nbrs.kneighbors(salesXY)

sales["logDistSubway"] = np.log10(subwayDists.mean(axis=1))

#Three Nearest Restaurants, Pubs, and Restaurants
restaurant = ox.geometries_from_polygon(city_limits_outline, tags={"amenity": ["restaurant", "pub", "bar"]})

restaurant = restaurant.to_crs(epsg=3857)

restaurantXY = xy_from_geom(restaurant.to_crs(epsg=3857))

restaurant_nbrs = NearestNeighbors(n_neighbors=3)

restaurant_nbrs.fit(restaurantXY)

restaurantDists, restaurantIndices = restaurant_nbrs.kneighbors(salesXY)

sales["logDistRestaurant"] = np.log10(restaurantDists.mean(axis=1))

#Three Nearest Arts Centers and Theatres
arts = ox.geometries_from_polygon(city_limits_outline, tags={"amenity": ["arts_centre", "theatre"]})

arts = arts.to_crs(epsg=3857)

artsXY = xy_from_geom(arts.to_crs(epsg=3857))

arts_nbrs = NearestNeighbors(n_neighbors=3)

arts_nbrs.fit(artsXY)

artsDists, artsIndices = arts_nbrs.kneighbors(salesXY)

sales["logDistArts"] = np.log10(artsDists.mean(axis=1))

#Three Nearest Parks, Gardens, and Playgrounds
parks = ox.geometries_from_polygon(city_limits_outline, tags={"leisure":["park", "garden", "playground"]})

parks = parks.to_crs(epsg=3857)

parksXY = xy_from_geom(parks.to_crs(epsg=3857))

parks_nbrs = NearestNeighbors(n_neighbors=3)

parks_nbrs.fit(parksXY)

parksDists, parksIndices = parks_nbrs.kneighbors(salesXY)

sales["logDistParks"] = np.log10(parksDists.mean(axis=1))

#Nearest Golf Course
golf = ox.geometries_from_polygon(city_limits_outline, tags={"leisure": "golf_course"})

golf = golf.to_crs(epsg=3857)

golfXY = xy_from_geom(golf.to_crs(epsg=3857))

golf_nbrs = NearestNeighbors(n_neighbors=1)

golf_nbrs.fit(golfXY)

golfDists, golfIndices = golf_nbrs.kneighbors(salesXY)

sales["logDistGolf"] = np.log10(golfDists.mean(axis=1))

sales

nbrs = NearestNeighbors(n_neighbors=5)
nbrs.fit(salesXY)
saleDists, saleIndex = nbrs.kneighbors(salesXY)

price = sales.iloc[saleIndex[0]]
price['sale_price'].mean()

from pandas import DataFrame

price_df = DataFrame(price, columns=['lagprice'])

sales['lagprice'] = price

sales

#Regression stated I had NaN values that RandomForestRegressor can not read, new df ensures no features have NaN values
#Ultimately removed about 200 rows
sales_final = sales[sales.columns[sales.columns.isin([
    'market_value', 
    'number_of_bathrooms',
    'number_of_bedrooms',
    'number_stories',
    'logDistGraffiti',
    'logDistFire',
    'logDistSchools',
    'logDistSubway',
    'logDistRestaurant',
    'logDistArts',
    'logDistParks',
    'logDistGolf',
    'lagprice',
    'exterior_condition',
    'interior_condition',
    'mapname',
    'topography',
    'sale_price',
    'NAME10',
    'geometry'])]].dropna()
sales_final