Merging NYC Housing Data with Geolocations

import pandas as pd
import geopandas as gpd
import matplotlib.pylab as plt
import fiona
import descartes
from shapely.geometry import Point, Polygon
%matplotlib inline

# Use your own path and file name here
data = pd.read_csv('data/rollingsales_brooklyn.csv')

# Take a look at the data
data.head()

# Eliminate outliers
n = 1 # defines how many stand deviation to look for. We find n=1 is the best practice.
# Establish acceptable upper bound and lower bound for geocoordinates
lon_upper = data['Longitude'].mean() + n * data['Longitude'].std()
lon_lower = data['Longitude'].mean() - n * data['Longitude'].std()
lat_upper = data['Latitude'].mean() + n * data['Latitude'].std()
lat_lower = data['Latitude'].mean() - n * data['Latitude'].std()
# Extract those records fall within the bounds
ind = (data['Longitude'] > lon_lower) &\
 (data['Longitude'] < lon_upper) &\
 (data['Latitude'] > lat_lower) &\
 (data['Latitude'] < lat_upper)
# Copy to a a new dataframe
data2 = data.loc[ind].copy()

# Read shape file and remember to use your own path
tzs=gpd.read_file('data/ZIP_CODE_040114/ZIP_CODE_040114.shp')

# Transfer tzs to 'epsg:4326' coordinates
# 'epsg:4326' is the type of geolocation encoding method we used in our data set
tzs = tzs.to_crs({'init':'epsg:4326'})

# Merge tzs with our housing data
geometry = [Point(xy) for xy in zip(data2['Longitude'],data2['Latitude'])] 
crs = {'init':'epsg:4326'}
geo_df = gpd.GeoDataFrame(data2, crs = crs, geometry = geometry)

# Take a look at what we have now
geo_df.head()

# Take a look at the housing sales' location.
fig, ax = plt.subplots(figsize = (12, 12))
tzs.plot(ax = ax, alpha = 1, color = '#cacfd9') # pick the color you want
ax = geo_df.plot(ax = ax, color = '#db7846', marker = 'o', markersize = 0.05)
ax.set_xlabel("Latitude")
ax.set_ylabel("Longitude")

# Just like what we had done before, we formulate a geo file for Tandon here.
# Geocoordinate can be easily find through Google map.
c_tandon = pd.DataFrame([[40.693579, -73.985701]],columns = ['Latitude','Longitude']) 
geometry_tandon = [Point(xy) for xy in zip(c_tandon['Longitude'],c_tandon['Latitude'])]
geo_tandon_df = gpd.GeoDataFrame(c_tandon, crs = crs, geometry = geometry_tandon)
ax.set_xlabel("Latitude")
ax.set_ylabel("Longitude")

lon_dis = 52.8 # 1 unit of longitude = 52.8 miles (around NYC's latitude)
lat_dis = 69.0 # 1 unit of latitude = 69.0 miles (almost a constant everywhere)

# calculate straight line distance
dis_lon = (geo_df['Longitude'] - c_tandon['Longitude'][0]) * lon_dis
dis_lat = (geo_df['Latitude'] - c_tandon['Latitude'][0]) * lat_dis
geo_df['dis_to_Tandon'] = (dis_lon ** 2 + dis_lat ** 2) ** 0.5

# seperate into 6 groups, you can do any other numbers
dis_group = [0, 0.5, 1, 2, 5, 10, 20]
legend = list()

fig, ax = plt.subplots(figsize = (12, 12))
tzs.plot(ax = ax, alpha = 1, color = '#cacfd9')

for i in range(0, len(dis_group) - 1):
    cond = (geo_df['dis_to_Tandon']>dis_group[i]) & (geo_df['dis_to_Tandon'] < dis_group[i + 1])
    geo_df[cond].plot(ax = ax, marker = 'o', markersize = 0.05)
    legend.append(str(dis_group[i]) + 'mi < d < ' + str(dis_group[i + 1]) + 'mi')
    
plt.legend(legend, markerscale = 20)