Exploring the NDVI in Philadelphia

Xiong Zheng

In this project, we'll explore the NDVI in Philadelphia, including two parts:

1.Compare the median NDVI within the city limits and the immediate suburbs
2.Calculate the NDVI around street trees in the city

Code

import matplotlib
import pandas as pd
import geopandas as gpd
import hvplot.pandas
import geoviews as gv
import geoviews.tile_sources as gvts
import holoviews as hv
import numpy as np
import rasterio as rio
import contextily as cx
hv.extension('bokeh', 'matplotlib')
from matplotlib import pyplot as plt
from holoviews import opts
from rasterio.mask import mask
from rasterstats import zonal_stats
from mpl_toolkits.axes_grid1 import make_axes_locatable
%matplotlib inline
pd.options.display.max_columns = 999
pd.options.display.max_rows = 999

2.1 Comparing the NDVI in the city and the suburbs

In this section, we will separate the city area from the suburbs, and calculate the NDVI for the city and the suburbs. After that, we will calculate the median NDVI within the city and within the suburbs

Code

# Separate
landsat = rio.open("./data/landsat8_philly.tif")

c_limit = gpd.read_file("./data/City_Limits.geojson")
c_limit = c_limit.to_crs(landsat.crs.data['init'])

env = c_limit.envelope
envgdf = gpd.GeoDataFrame(gpd.GeoSeries(env))
envgdf = envgdf.rename(columns={0:'geometry'}).set_geometry('geometry')
envgdf = envgdf.to_crs(landsat.crs.data['init'])

sub_limit = gpd.overlay(envgdf,c_limit, how='difference')

city, mask_transform = mask(
    dataset=landsat,
    shapes=c_limit.geometry,
    crop=True,  
    all_touched=True,  
    filled=False, 
)

suburb, mask_transform = mask(
    dataset=landsat,
    shapes=sub_limit.geometry,
    crop=True,  
    all_touched=True,  
    filled=False, 
)

c_red = city[3]
c_nir = city[4]
sub_red = suburb[3]
sub_nir = suburb[4]

# Calculate
def calculate_NDVI(nir, red):
    """
    Calculate the NDVI from the NIR and red landsat bands
    """
    nir = nir.astype(float)
    red = red.astype(float)
    
    check = np.logical_and( red.mask == False, nir.mask == False )
    
    ndvi = np.where(check,  (nir - red ) / ( nir + red ), np.nan )
    return ndvi

city_NDVI = calculate_NDVI(c_nir, c_red)
sub_NDVI = calculate_NDVI(sub_nir, sub_red)

fig, ax = plt.subplots(figsize=(10,10))
landsat_extent = [
    landsat.bounds.left,
    landsat.bounds.right,
    landsat.bounds.bottom,
    landsat.bounds.top,
]
img = ax.imshow(city_NDVI,extent=landsat_extent)
c_limit.boundary.plot(ax=ax, edgecolor='gray', facecolor='none', linewidth=4)

plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI in City of Philadelphia", fontsize=18);

Code

med_sub = np.nanmedian(sub_NDVI)
med_city = np.nanmedian(city_NDVI)
print(" Conclusion: suburbs has a higher NDVI","\n","Median NDVI within the city:",med_city,"\n","Median NDVI within the suburbs:",med_sub)

 Conclusion: suburbs has a higher NDVI 
 Median NDVI within the city: 0.20268593532493442 
 Median NDVI within the suburbs: 0.3746654463028859

2.2 Calculating the NDVI for Philadelphia's street treets

Code

st_trees = gpd.read_file("./data/ppr_tree_canopy_points_2015.geojson")
st_trees = st_trees.to_crs(landsat.crs.data['init'])

tree_NDVI = zonal_stats(
    st_trees, 
    city_NDVI, 
    affine=landsat.transform, 
    stats=['mean'])

st_trees['NDVI'] = [s['mean'] for s in tree_NDVI] 

# Plot
fig, ax = plt.subplots(figsize=(8,6))

plt.grid(color='#A8A8A8', lw=0.5,linestyle='dashed')

ax.hist(st_trees['NDVI'], bins=50,color='#FCBB85')

ax.axvline(x=0, c='k', lw=2,linestyle='dashed')
ax.set_xlabel("Tree NDVI", fontsize=18)
ax.set_ylabel("Number of Trees", fontsize=18);
ax.set_title("Tree NDVI in Philadelphia", fontsize=18);

Code

evPA = gpd.read_file("./data/PA-tracts.geojson")
evphilly = evPA.loc[(evPA['pl']=="Philadelphia County, Pennsylvania")]
needcolumns = ['GEOID','geometry']
value_vars = ['e-{:02d}'.format(x) for x in range(3, 17)]
needcolumns.extend(value_vars)
evphilly = evphilly[needcolumns]
mlt_evphilly = evphilly.melt(
    id_vars=['GEOID','geometry'],
    value_vars=value_vars,
    var_name='year', 
    value_name='evictions')

geo_evphilly = mlt_evphilly.loc[(mlt_evphilly['year']=="e-03")].drop(['year','evictions'],axis=1)
geo_philly = geo_evphilly.to_crs(c_limit.crs)

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

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

c_limit.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=4)

geo_philly.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=0.1)

st_trees.plot(
    column='NDVI', 
    legend=True, 
    ax=ax,
    cax=cax,
    cmap='OrRd',
    s = 55,
    alpha = 0.5)

cx.add_basemap(ax,
               zoom=12, 
               crs=st_trees.crs)

plt.title("Brooklyn",fontsize=10)
ax.set_title("Tree NDVI in Philadelphia EPSG:32618", fontsize=18);