import rasterio as rio
import matplotlib.colors as mcolors
import earthpy.plot as ep
from rasterio.mask import mask
import geopandas as gpd
from matplotlib import pyplot as plt
import numpy as np
from rasterstats import zonal_stats
import warnings

warnings.filterwarnings('ignore')

landsat = rio.open("Desktop/School/MUSA 550/assignment-3-main/data/landsat8_philly.tif")

band1 = landsat.read(1)

cityLimits = gpd.read_file(
    "Desktop/School/MUSA 550/assignment-3-main/data/City_Limits.geojson")

cityLimits = cityLimits.to_crs(landsat.crs.data['init'])

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

landsat_extent = [
    landsat.bounds.left,
    landsat.bounds.right,
    landsat.bounds.bottom,
    landsat.bounds.top,
]

img = ax.imshow(band1, 
                norm=mcolors.LogNorm(), 
                extent=landsat_extent)  

cityLimits.plot(ax=ax, facecolor="none", edgecolor="white", linewidth=2)

plt.colorbar(img)
ax.set_axis_off()

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

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

ax.imshow(masked[0], cmap="viridis", extent=landsat_extent)

cityLimits.boundary.plot(ax=ax, color="gray", linewidth=3)
ax.set_axis_off()

envelope = cityLimits.envelope

suburbs = envelope.difference(cityLimits)

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

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

ax.imshow(maskedSuburb[0], cmap="viridis", extent=landsat_extent)

cityLimits.boundary.plot(ax=ax, color="gray", linewidth=3)
ax.set_axis_off()

redCity = masked[3]
nirCity = masked[4]

def calculate_NDVI(nirCity, redCity):
  
    nirCity = nirCity.astype(float)
    redCity = redCity.astype(float)
    
    check = np.logical_and( redCity.mask == False, nirCity.mask == False )
    
    ndvi = np.where(check,  (nirCity - redCity ) / ( nirCity + redCity ), np.nan )
    return ndvi 

NDVI_city = calculate_NDVI(nirCity, redCity)

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

img = ax.imshow(NDVI_city, extent=landsat_extent)

cityLimits.plot(ax=ax, edgecolor='gray', facecolor='none', linewidth=3)
plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI in Philadelphia", fontsize=18);

redSuburbs = maskedSuburb[3]
nirSuburbs = maskedSuburb[4]

def calculate_NDVI(nirSuburbs, redSuburbs):
    
    nirSuburbs = nirSuburbs.astype(float)
    redSuburbs = redSuburbs.astype(float)
  
    check = np.logical_and( redSuburbs.mask == False, nirSuburbs.mask == False )
 
    ndvi = np.where(check,  (nirSuburbs - redSuburbs ) / ( nirSuburbs + redSuburbs ), np.nan )
    return ndvi 

NDVI_suburbs = calculate_NDVI(nirSuburbs, redSuburbs)

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

img = ax.imshow(NDVI_suburbs, extent=landsat_extent)

cityLimits.plot(ax=ax, edgecolor='gray', facecolor='none', linewidth=3)
plt.colorbar(img)
ax.set_axis_off()
ax.set_title("NDVI in Suburbs", fontsize=18);

statsCity = zonal_stats(cityLimits, NDVI_city, affine=landsat.transform, stats=['median'])

statsSuburbs = zonal_stats(cityLimits, NDVI_suburbs, affine=landsat.transform, stats=['median'])

print(statsCity, statsSuburbs)

[{'median': 0.20246571107834677}] [{'median': 0.2682237997099458}]

The median NDVI in Philadelphia suburbs is slightly higher than the median NDVI within the city limits.

trees = gpd.read_file(
    "Desktop/School/MUSA 550/assignment-3-main/data/ppr_tree_canopy_points_2015.geojson")

trees = trees.to_crs(landsat.crs.data['init'])

treeStats = zonal_stats(trees, NDVI_city, affine=landsat.transform, stats=['mean'])

trees['mean_NDVI'] = [s['mean'] for s in treeStats]

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

ax.hist(trees['mean_NDVI'], bins='auto')
ax.axvline(x=0, c='k', lw=2)

ax.set_xlabel("NDVI", fontsize=18)
ax.set_ylabel("Number of Trees", fontsize=18);

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

ax.imshow(NDVI_city, cmap="viridis", extent=landsat_extent, alpha = 0.25)
cityLimits.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=3)
trees.plot(column='mean_NDVI', legend=True, ax=ax, cmap="viridis")

ax.set_axis_off()
ax.set_title("NDVI of Select Tree Locations", fontsize=18)

Text(0.5, 1.0, 'NDVI of Select Tree Locations')