%reload_ext watermark
%watermark -v -p numpy,pandas,geopandas,matplotlib,matplotlib,seaborn,contextily

# list items in data folder
%ls data/

# read data as a dataframe
path = 'data/street-flooding-complaints.csv'
df = pd.read_csv(path, low_memory=False)

# preview data
print(f'shape of data: {df.shape}')
df.head()

# column info
df.info()

# importing nta boundaries
url = 'https://data.cityofnewyork.us/resource/9nt8-h7nd.geojson'
nta_gdf = gpd.read_file(url).to_crs(epsg=2263)

# full name
nta_gdf['ntaname_full'] = nta_gdf['ntaname'] + ', ' + nta_gdf['boroname']

# previewing first five rows in data
print(f'shape of data: {nta_gdf.shape}')
nta_gdf.head()

# sanity check
nta_gdf.plot()

path = 'data/streets.gpkg'
streets = gpd.read_file(path)

# sanity checks
print(f'shape of data: {streets.shape}')
print(f"street id is unique: {streets['physicalid'].is_unique}")
print(streets.crs)

streets.head()

# column info
streets.info()

# sanity check, examine counts of geom types
streets.geom_type.value_counts()

# sanity check plot
streets.plot()

crs=2263
geometry = gpd.points_from_xy(
    df['x_coordinate_state_plane'],
    df['y_coordinate_state_plane']
)

# make geodataframe form x,y points
gdf = gpd.GeoDataFrame(df, geometry=geometry, crs=crs)

# spatial join with nta information
gdf = gpd.sjoin(
    gdf, 
    nta_gdf, 
    how="inner", 
    predicate='within'
)

print('shape of data: {}'.format(gdf.shape))
gdf.head()

(gdf
 .groupby(by='ntaname_full')['unique_key']
 .count()
 .rename('count')
 .reset_index()
 .sort_values('count', ascending=False)
 .head(10)
)

crs = 2263
geometry = gpd.points_from_xy(
    df['x_coordinate_state_plane'],
    df['y_coordinate_state_plane']
)

# make geodataframe from x,y points
gdf = gpd.GeoDataFrame(df, geometry=geometry, crs=crs)

# spatial join with nta information
gdf = gpd.sjoin(
    gdf, 
    nta_gdf, 
    how="inner", 
    predicate='within'
)

print('shape of data: {}'.format(gdf.shape))
gdf.head()

# locate midland beach nta shape
name = 'New Dorp-Midland Beach, Staten Island'
midland_beach_shape = (
    nta_gdf
    .loc[nta_gdf['ntaname_full'] == name]
    .reset_index(drop=True)
)

# preview data
midland_beach_shape.head()

# clipping complaints that fall within Midland Beach
midland_beach_gdf = gpd.clip(gdf, midland_beach_shape)

# sanity check plot
midland_beach_gdf.plot()

# clip streets within Midland Beach
streets_clipped = gpd.clip(streets, midland_beach_shape)

# recalculate street length
streets_clipped['shape_leng'] = streets_clipped['geometry'].length

# sanity check plot
print(f'shape of data: {streets_clipped.shape}')
streets_clipped.plot()

# examine geom_type values
streets_clipped.geom_type.value_counts()

# retrieve only linestring geom types
streets_clipped = (
    streets_clipped.loc[streets_clipped.geom_type == 'LineString']
    .reset_index(drop=True)
)

# sanity check geom types
streets_clipped.geom_type.value_counts()

# offset 80 ft.
offset = 80
bbox = midland_beach_gdf.bounds + [-offset, -offset, offset, offset]
hits = bbox.apply(lambda row: list(streets_clipped.sindex.intersection(row)), axis=1)

hits.head()

tmp = pd.DataFrame({
    "pt_idx": np.repeat(hits.index, hits.apply(len)),
    "line_i": np.concatenate(hits.values)
})

# join back to the lines on line_i; we use reset_index() to 
# give us the ordinal position of each line
tmp = tmp.join(
    streets_clipped.reset_index(drop=True),
    on="line_i"
)

# join back to the original points to get their geometry
# rename the point geometry as "point"
tmp = tmp.join(
    midland_beach_gdf['geometry'].rename("point"),
    on="pt_idx"
)

# convert back to a GeoDataFrame, so we can do spatial ops
tmp = gpd.GeoDataFrame(
    tmp,
    geometry="geometry",
    crs=midland_beach_gdf.crs
)

tmp.head()

# calculate distance from line to point
tmp["snap_dist"] = tmp['geometry'].distance(gpd.GeoSeries(tmp.point))

# discard any lines that are greater than tolerance from points
tmp = (
    tmp
    .loc[tmp.snap_dist <= offset]
    .sort_values(by=["snap_dist"])
)

tmp.head()

# group by the index of the points and take the first, which is the closest line 
closest = tmp.groupby("pt_idx").first()

# construct a GeoDataFrame of the closest lines
closest = gpd.GeoDataFrame(closest, geometry="geometry")

closest.head()

# position of nearest point from start of the line
pos = closest['geometry'].project(gpd.GeoSeries(closest['point']))

# Get new point location geometry
new_pts = closest.geometry.interpolate(pos)

# Create a new GeoDataFrame from the columns from the closest line and
# new point geometries (which will be called "geometries")
snapped = gpd.GeoDataFrame(closest, geometry=new_pts)

snapped.head()

# join back to the original points
updated_points = (
    midland_beach_gdf.drop(columns=['geometry', 'shape_leng'])
    .join(snapped.drop(columns=['borocode']))
    .dropna(subset=["geometry"])
)

updated_points.head()

updated_points = gpd.GeoDataFrame(
    updated_points,
    geometry='point'
)

updated_points.plot()

fig, ax = plt.subplots(figsize=(8, 8)) 

updated_points.plot(
    color='black', 
    edgecolor='black',
    alpha=.4,
    ax=ax,
    zorder=3
)

streets_clipped.plot(ax=ax, zorder=2)

midland_beach_shape.plot(
    ax=ax, 
    zorder=1, 
    color='None', 
    edgecolor='black'
)

# adding basemap
ctx.add_basemap(
    ax, 
    crs=2263, 
    source=ctx.providers.CartoDB.PositronNoLabels,
    zorder=0
)

label = 'NYC 311 Street Flooding Complaints\n\
in New Dorp-Midland Beach, Staten Island from 2010 to 2020'
ax.set_title(label, fontsize=12, pad=10)

plt.axis('off')
plt.tight_layout()

# buffer lines
streets_clipped = (
    streets_clipped
    .set_geometry('geometry')
    .assign(new_geom=lambda x: x.geometry.buffer(40, cap_style=2))
    .set_geometry('new_geom')
)

streets_clipped.head()

# group count of complaints by street
gdf_count = (
    updated_points
    .groupby(by='physicalid')['created_date']
    .count()
    .reset_index()
    .rename(columns={"created_date": "count"})
)
gdf_count.head()

# joining count of complaints to streets
streets_with_count = streets_clipped.merge(
    gdf_count, 
    left_on='physicalid', 
    right_on='physicalid', 
    how='left'
)

streets_with_count['count'] = streets_with_count['count'].fillna(0).astype(int)

# examine df
streets_with_count.info()

# normalize counts
streets_with_count['count_per_100ft'] = (streets_with_count['count'] / streets_with_count['shape_leng']) * 100

streets_with_count.head()

# exammine distribution of counts
streets_with_count.loc[:, ['count', 'count_per_100ft']].describe()

fig, ax = plt.subplots(figsize=(8, 8)) 

divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cmap = plt.cm.Blues

# get max of array and place as top bounds
norm = mpl.colors.BoundaryNorm([0, 1, 5, 10, 30, 60], cmap.N)

streets_with_count.plot(
    column='count', 
    ax=ax,
    cax=cax,
    zorder=2,                        
    cmap=cmap,
    norm=norm,
    edgecolor='black',
    linewidth=.3,
    legend=True
)

(midland_beach_shape
 .plot(
     ax=ax, 
     zorder=1, 
     color='None', 
     edgecolor='black')
)

# adding basemap
ctx.add_basemap(
    ax, 
    crs=2263, 
    source=ctx.providers.CartoDB.PositronNoLabels,
    zorder=0
)


# setting title
label = 'Count of 311 Street Flooding Complaints by Street Segment\n\
in New Dorp-Midland Beach, Staten Island from 2010 to 2020'
ax.set_title(label, fontsize=12, pad=10)

ax.axis('off')
plt.tight_layout()

fig, ax = plt.subplots(figsize=(8, 8)) 

divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cmap = plt.cm.Blues

# get max of array and place as top bounds
norm = mpl.colors.BoundaryNorm([0, 1, 3, 5, 10, 30], cmap.N)

streets_with_count.plot(
    column='count_per_100ft', 
    ax=ax,
    cax=cax,
    zorder=2,                        
    cmap=cmap,
    norm=norm,
    edgecolor='black',
    linewidth=.3,
    legend=True
)

(midland_beach_shape
 .plot(
     ax=ax, 
     zorder=1, 
     color='None', 
     edgecolor='black')
)

# adding basemap
ctx.add_basemap(
    ax, 
    crs=2263, 
    source=ctx.providers.CartoDB.PositronNoLabels,
    zorder=0
)


# setting title
label = 'Count of 311 Street Flooding Complaints per 100 ft.\n\
by Street Segment in New Dorp-Midland Beach, Staten Island from 2010 to 2020'
ax.set_title(label, fontsize=12, pad=10)

ax.axis('off')
plt.savefig('figures/midland-beach.png', dpi=250, bbox_inches='tight')

# create full street name
streets_with_count['street_and_id'] = (
    streets_with_count['full_stree']
    + ', id:'
    + streets_with_count['physicalid']
)

streets_with_count.head()

fig, axs = plt.subplots(2, 1, sharey=False, figsize=(8, 8))

# first plot
data = streets_with_count.sort_values(by='count', ascending=False).head(10)
sns.barplot(
    data=data,
    y='street_and_id',
    x='count',
    color='#1f77b4',
    ax=axs[0]
)

label = 'Count of NYC 311 Street Flooding Complaints\n\
by Street Segment in New Dorp-Midland Beach, Staten Island from 2010 to 2020'
axs[0].set_title(label, x=.3)
axs[0].set_xlabel('Count')
axs[0].set_ylabel('Street Segment')

# second plot
data = streets_with_count.sort_values(by='count_per_100ft', ascending=False).head(10)
sns.barplot(
    data=data,
    y='street_and_id',
    x='count_per_100ft',
    color='#1f77b4',
    ax=axs[1]
)

label = 'Count of NYC 311 Street Flooding Complaints per 100 ft.\n\
by Street Segment in New Dorp-Midland Beach, Staten Island from 2010 to 2020'
axs[1].set_title(label, x=.3)
axs[1].set_xlabel('Count per 100 ft.')
axs[1].set_ylabel('Street Segment')

fig.tight_layout()

# transform complaints to geodataframe
crs = 2263
geometry = gpd.points_from_xy(
    df['x_coordinate_state_plane'],
    df['y_coordinate_state_plane']
)

gdf = gpd.GeoDataFrame(df, geometry=geometry, crs=crs)

# sanity check
gdf.head()

# streets
path = 'data/streets.gpkg'
streets = gpd.read_file(path)

streets.head()

# examine streets df
streets.info()

# return only Hammels-Arverne-Edgemere, Queens
name = 'Rockaway Beach-Arverne-Edgemere, Queens'
arverne_edgemere_shape = (
    nta_gdf
    .loc[nta_gdf['ntaname_full'] == name]
    .reset_index(drop=True)
)

arverne_edgemere_shape

# clipping complaints that fall within arverne edgemere
arverne_edgemere_gdf = gpd.clip(
    gdf,
    arverne_edgemere_shape
)

arverne_edgemere_gdf.plot()

# clipping the points that only fall within the nta shapes
streets_clipped = gpd.clip(
    streets,
    arverne_edgemere_shape
)

# sanity check plot
streets_clipped.plot()

# sanity check geom type
streets_clipped.geom_type.value_counts()

# locating only linestring geom type
streets_clipped = (
    streets_clipped
    .loc[streets_clipped.geom_type == 'LineString']
    .reset_index(drop=True)
)

streets_clipped.geom_type.value_counts()

# sanity check plot
streets_clipped.plot()

# offset by 80 ft.
offset = 80
bbox = arverne_edgemere_gdf.bounds + [-offset, -offset, offset, offset]
hits = bbox.apply(lambda row: list(streets_clipped.sindex.intersection(row)), axis=1)

tmp = pd.DataFrame({
    "pt_idx": np.repeat(hits.index, hits.apply(len)),
    "line_i": np.concatenate(hits.values)
})

tmp.head()

# join back to the lines on line_i; we use reset_index() to 
# give us the ordinal position of each line
tmp = tmp.join(
    streets_clipped.reset_index(drop=True), 
    on="line_i"
)

tmp.head()

# join back to the original points to get their geometry
# rename the point geometry as "point"
tmp = tmp.join(
    arverne_edgemere_gdf.geometry.rename("point"), 
    on="pt_idx"
)

# convert back to a GeoDataFrame, so we can do spatial ops
tmp = gpd.GeoDataFrame(
    tmp, 
    geometry="geometry", 
    crs=arverne_edgemere_gdf.crs
)

tmp.head()

tmp["snap_dist"] = tmp['geometry'].distance(gpd.GeoSeries(tmp.point))

# discard any lines that are greater than tolerance from points
# sort on ascending snap distance, so that closest goes to top
tmp = (
    tmp
    .loc[tmp.snap_dist <= offset]
    .sort_values(by=["snap_dist"])
)

tmp.head()

# group by the index of the points and take the first, which is the closest line 
closest = tmp.groupby("pt_idx").first()
closest = gpd.GeoDataFrame(closest, geometry="geometry")

# position of nearest point from start of the line
pos = closest.geometry.project(gpd.GeoSeries(closest.point))

# get new point location geometry
new_pts = closest.geometry.interpolate(pos)

new_pts.head()

# create a new GeoDataFrame from the columns from the closest line and 
# new point geometries (which will be called "geometries")
snapped = gpd.GeoDataFrame(closest, geometry=new_pts)

snapped.head()

# join back to the original points:
updated_points = (
    arverne_edgemere_gdf
    .drop(columns=["geometry"])
    .join(snapped)
    .dropna(subset=["geometry"])
)

updated_points.head()

updated_points.loc[:, ['snap_dist']].describe()

plt.figure(figsize=(6, 4))

sns.histplot(updated_points.snap_dist)
plt.title('Histogram of snap_dist (ft.)')

updated_points = gpd.GeoDataFrame(
    updated_points,
    geometry='point'
)

updated_points.plot()

fig, ax = plt.subplots(figsize=(10, 10))

updated_points.plot(
    color='black',
    edgecolor='black',
    alpha=.4,
    ax=ax,
    zorder=3
)

streets_clipped.plot(ax=ax, zorder=2)

(arverne_edgemere_shape
 .plot(ax=ax, 
       zorder=1, 
       color='None', 
       edgecolor='black')
)

# adding basemap
ctx.add_basemap(ax, crs=2263, source=ctx.providers.CartoDB.PositronNoLabels, zorder=0)

label = 'NYC 311 Street Flooding Complaints\nin Rockaway Beach-Arverne-Edgemere, \
Queens, New York from 2010 to 2020'
ax.set_title(label, fontsize=12, pad=10)

plt.axis('off')
plt.tight_layout()

# count of complaints by street
gdf_count = (
    updated_points
     .groupby(by='physicalid')['created_date']
     .count()
     .reset_index()
     .rename(columns={"created_date": "count"})
)

gdf_count.head()

# buffer lines
streets_clipped = (
    streets_clipped
    .set_geometry('geometry')
    .assign(new_geom=lambda x: x.geometry.buffer(40, cap_style=2))
    .set_geometry('new_geom')
)

streets_clipped.head()

# joining our nta population data to our nta shapes data
streets_with_count = streets_clipped.merge(
    gdf_count, 
    left_on='physicalid', 
    right_on='physicalid', 
    how='left'
)

streets_with_count['shape_leng'] = streets_with_count['geometry'].length
streets_with_count['count'] = streets_with_count['count'].fillna(0).astype(int)

streets_with_count.info()

# normalize counts
streets_with_count['count_per_100ft'] = (streets_with_count['count'] / streets_with_count['shape_leng']) * 100

streets_with_count.head()

streets_with_count.loc[:, ['count', 'count_per_100ft']].describe()

fig, ax = plt.subplots(figsize=(10, 10))

divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cmap = plt.cm.Blues

# get max of array and place as top bounds
norm = mpl.colors.BoundaryNorm([0, 1, 5, 10, 30, 50], cmap.N)

streets_with_count.plot(
    column='count', 
    ax=ax,
    cax=cax, 
    norm=norm,
    zorder=2,                        
    cmap='Blues',
    edgecolor='black',
    linewidth=.3,
    legend=True
)

(arverne_edgemere_shape
 .plot(ax=ax, 
       zorder=1, 
       color='None', 
       edgecolor='black')
)

# adding basemap
ctx.add_basemap(ax, crs=2263, source=ctx.providers.CartoDB.PositronNoLabels, zorder=0)

label = 'Count of NYC 311 Street Flooding Complaints \n\
by Street Segment in Rockaway Beach-Arverne-Edgemere, Queens from 2010 to 2020'
ax.set_title(label=label, fontsize=12, pad=10)

ax.axis('off')
plt.tight_layout()

fig, ax = plt.subplots(figsize=(10, 10))

divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cmap = plt.cm.Blues

# get max of array and place as top bounds
norm = mpl.colors.BoundaryNorm([0, 1, 3, 5, 10, 30], cmap.N)

streets_with_count.plot(
    column='count', 
    ax=ax,
    cax=cax, 
    norm=norm,
    zorder=2,                        
    cmap='Blues',
    edgecolor='black',
    linewidth=.3,
    legend=True
)

(arverne_edgemere_shape
 .plot(ax=ax, 
       zorder=1, 
       color='None', 
       edgecolor='black')
)

# adding basemap
ctx.add_basemap(ax, crs=2263, source=ctx.providers.CartoDB.PositronNoLabels, zorder=0)

label = 'Count of NYC 311 Street Flooding Complaints per 100 ft.\n\
by Street Segment in Rockaway Beach-Arverne-Edgemere, Queens from 2010 to 2020'
ax.set_title(label=label, fontsize=12, pad=10)

ax.axis('off')
plt.savefig('figures/arverne.png', dpi=250, bbox_inches='tight')

streets_with_count['street_and_id'] = (
    streets_with_count['full_stree']
    + ', id:'
    + streets_with_count['physicalid']
)

streets_with_count.head()

fig, axs = plt.subplots(2, 1, sharey=False, figsize=(8, 8))

# first plot
data = streets_with_count.sort_values(by='count', ascending=False).head(10)
sns.barplot(
    data=data,
    y='street_and_id',
    x='count',
    color='#1f77b4',
    ax=axs[0]
)

label = 'Count of NYC 311 Street Flooding Complaints\n\
by Street Segment in Hammels-Arverne-Edgemere, Queens from 2010 to 2020'
axs[0].set_title(label, x=.3)
axs[0].set_xlabel('Count')
axs[0].set_ylabel('Street Segment')

# second plot
data = streets_with_count.sort_values(by='count_per_100ft', ascending=False).head(10)
sns.barplot(
    data=data,
    y='street_and_id',
    x='count_per_100ft',
    color='#1f77b4',
    ax=axs[1]
)

label = 'Count of NYC 311 Street Flooding Complaints per 100 ft.\n\
by Street Segment in Hammels-Arverne-Edgemere, Queens from 2010 to 2020'
axs[1].set_title(label, x=.3)
axs[1].set_xlabel('Count per 100 ft.')
axs[1].set_ylabel('Street Segment')

plt.tight_layout()