CODE: NLCD Raster Calculations

Hot City, Heated Calls:
Understanding Extreme Heat and Quality of Life
Using New York City's 311 and SHAP

Zonal statistics for tree canopy and impervious surface from NLCD data.

← 4. Additional Features 6. Data Merging →

NLCD RASTERS

Extract tree canopy and impervious percentages.
Code Cell 1 
# Modules.
import os
from pathlib import Path
import numpy as np
import rasterio
from rasterio.mask import mask
from rasterio.warp import reproject, Resampling
import geopandas as gpd
from tqdm import tqdm
from rasterstats import zonal_stats
Code Cell 2 
# Paths.
nlcd_tree_path = Path("data/raster/nlcd_raster/nlcd_tree_canopy_2023.tiff")
nlcd_impervious_path = Path("data/raster/nyc_impervious_2024.tif")

tracts_path = Path("data/nyc_tracts_2020/nyc_tracts_2020.shp")

output_dir = Path("data/raster/processed")
output_dir.mkdir(parents = True, exist_ok = True)

tracts = gpd.read_file(tracts_path).to_crs(32118)
Code Cell 3 
# Zonal statistics for NCLD.
def zonal_mean(rpath, gdf_or_geom):
    """Apply CRS zonal mean for tracts or city boundary."""
    with rasterio.open(rpath) as src:
        r_crs = src.crs

    if isinstance(gdf_or_geom, (gpd.GeoSeries, gpd.GeoDataFrame)):
        gdf = gdf_or_geom.to_crs(r_crs)
    else:
        gdf = gpd.GeoSeries([gdf_or_geom], crs = 32118).to_crs(r_crs)

    return zonal_stats(gdf, rpath, stats = ["mean"], nodata = np.nan)
Code Cell 4 
# Print checks for the calculations.
print("Tree canopy zonal stats:")
tree = zonal_mean(nlcd_tree_path, tracts)
print(tree)
print("Impervious zonal stats:")
impervious = zonal_mean(nlcd_impervious_path, tracts)
print(impervious)

tracts["pct_tree_canopy"] = [t["mean"] for t in tree]
tracts["pct_impervious"] = [i["mean"] for i in impervious]
Tree canopy zonal stats:
[{'mean': 10.126984126984127}, {'mean': 13.314285714285715}, {'mean': 0.592}, {'mean': 3.076305220883534}, {'mean': 6.089887640449438}, {'mean': 2.736842105263158}, {'mean': 1.7009803921568627}, {'mean': 1.2598870056497176}, {'mean': 6.25}, {'mean': 2.8217821782178216}, {'mean': 0.4077669902912621}, {'mean': 1.0520833333333333}, {'mean': 0.047619047619047616}, {'mean': 1.2023809523809523}, {'mean': 0.0}, {'mean': 7.0647482014388485}, {'mean': 0.07936507936507936}, {'mean': 0.05154639175257732}, {'mean': 1.3650793650793651}, {'mean': 0.15025906735751296}, {'mean': 6.22}, {'mean': 0.6020408163265306}, {'mean': 0.39896373056994816}, {'mean': 0.6956521739130435}, {'mean': 0.19270833333333334}, {'mean': 1.517766497461929}, {'mean': 4.198979591836735}, {'mean': 0.387434554973822}, {'mean': 0.10256410256410256}, {'mean': 0.05235602094240838}, {'mean': 0.8134715025906736}, {'mean': 4.455958549222798}, {'mean': 3.572972972972973}, {'mean': 0.0}, {'mean': 4.678391959798995}, {'mean': 0.0}, {'mean': 0.0}, {'mean': 0.0}, {'mean': 1.3679245283018868}, {'mean': 0.0}, {'mean': 0.0}, {'mean': 0.0}, {'mean': 0.16580310880829016}, {'mean': 0.6375}, {'mean': 0.15625}, {'mean': 0.0}, {'mean': 0.07142857142857142}, {'mean': 0.32954545454545453}, {'mean': 0.0}, {'mean': 0.0}, {'mean': 0.0}, {'mean': 0.6306306306306306}, {'mean': 2.5786802030456855}, {'mean': 0.0}, {'mean': 1.1614583333333333}, {'mean': 0.10204081632653061}, {'mean': 0.59}, {'mean': 0.953125}, {'mean': 0.05154639175257732}, {'mean': 1.5339805825242718}, {'mean': 1.606060606060606}, {'mean': 0.7680412371134021}, {'mean': 1.7868020304568528}, {'mean': 0.05699481865284974}, {'mean': 0.15294117647058825}, {'mean': 5.6}, {'mean': 0.7525773195876289}, {'mean': 0.65}, {'mean': 0.8571428571428571}, {'mean': 0.4482758620689655}, {'mean': 1.1649484536082475}, {'mean': 1.6551724137931034}, {'mean': 1.6}, {'mean': 8.90983606557377}, {'mean': 0.9375}, {'mean': 0.6632653061224489}, {'mean': 5.6875}, {'mean': 0.3442622950819672}, {'mean': 7.12280701754386}, {'mean': 11.072916666666666}, {'mean': 0.26732673267326734}, {'mean': 1.5602836879432624}, {'mean': 0.8392857142857143}, {'mean': 1.6906474820143884}, {'mean': 1.36}, {'mean': 9.269035532994923}, {'mean': 10.030837004405287}, {'mean': 2.9076923076923076}, {'mean': 0.465}, {'mean': 14.548936170212766}, {'mean': 1.8172588832487309}, {'mean': 1.803030303030303}, {'mean': 4.918918918918919}, {'mean': 1.9824561403508771}, {'mean': 13.605042016806722}, {'mean': 3.8969072164948453}, {'mean': 6.851694915254237}, {'mean': 1.367741935483871}, {'mean': 6.197969543147208}, {'mean': 12.355555555555556}, {'mean': 20.541666666666668}, {'mean': 1.6717948717948719}, {'mean': 1.8934010152284264}, {'mean': 1.537037037037037}, {'mean': 22.463414634146343}, {'mean': 1.8363636363636364}, {'mean': 0.46153846153846156}, {'mean': 0.8493975903614458}, {'mean': 1.4311377245508983}, {'mean': 1.1818181818181819}, {'mean': 22.27220630372493}, {'mean': 1.349753694581280
... [output truncated]
Code Cell 5 
tracts.columns = tracts.columns.str.upper()
tracts = tracts.rename(columns = {"GEOMETRY": "geometry"})

tracts.head()
CTLABEL BOROCODE BORONAME CT2020 BOROCT2020 CDELIGIBIL NTANAME NTA2020 CDTA2020 CDTANAME GEOID SHAPE_LENG SHAPE_AREA geometry PCT_TREE_CANOPY PCT_IMPERVIOUS
0 1 1 Manhattan 000100 1000100 I The Battery-Governors Island-Ellis Island-Libe... MN0191 MN01 MN01 Financial District-Tribeca (CD 1 Equivalent) 36061000100 10833.043929 1.843005e+06 MULTIPOLYGON (((296291.11 58135.714, 296322.49... 10.126984 46.973545
1 14.01 1 Manhattan 001401 1001401 I Lower East Side MN0302 MN03 MN03 Lower East Side-Chinatown (CD 3 Equivalent) 36061001401 5075.332000 1.006117e+06 POLYGON ((300982.972 61050.77, 301053.207 6102... 13.314286 69.590476
2 14.02 1 Manhattan 001402 1001402 E Lower East Side MN0302 MN03 MN03 Lower East Side-Chinatown (CD 3 Equivalent) 36061001402 4459.156019 1.226206e+06 POLYGON ((301261.149 61343.714, 301332.269 613... 0.592000 83.656000
3 18 1 Manhattan 001800 1001800 I Lower East Side MN0302 MN03 MN03 Lower East Side-Chinatown (CD 3 Equivalent) 36061001800 6391.921174 2.399277e+06 POLYGON ((300857.167 61503.238, 300868.535 614... 3.076305 81.887550
4 22.01 1 Manhattan 002201 1002201 E Lower East Side MN0302 MN03 MN03 Lower East Side-Chinatown (CD 3 Equivalent) 36061002201 5779.062607 1.740174e+06 POLYGON ((301795.202 61438.262, 301792.173 614... 6.089888 78.735955
Code Cell 6 
tracts.columns
Index(['CTLABEL', 'BOROCODE', 'BORONAME', 'CT2020', 'BOROCT2020', 'CDELIGIBIL',
       'NTANAME', 'NTA2020', 'CDTA2020', 'CDTANAME', 'GEOID', 'SHAPE_LENG',
       'SHAPE_AREA', 'geometry', 'PCT_TREE_CANOPY', 'PCT_IMPERVIOUS'],
      dtype='object')
Code Cell 7 
tracts = tracts.drop(columns = ['CTLABEL', 'BOROCODE', 'BORONAME', 'CT2020',
                       'BOROCT2020', 'CDELIGIBIL', 'NTANAME', 'NTA2020',
                       'CDTA2020', 'CDTANAME', 'SHAPE_LENG', 'SHAPE_AREA'])
Code Cell 8 
tracts.columns
Index(['GEOID', 'geometry', 'PCT_TREE_CANOPY', 'PCT_IMPERVIOUS'], dtype='object')
Code Cell 9 
# Save as geojson.
geojson_out = output_dir.parent / "nlcd_calc_tracts.geojson"
geojson_out.parent.mkdir(parents = True, exist_ok = True)
tracts.to_file(geojson_out)
print("Saved:", geojson_out)

# Save as csv.
csv_out = output_dir.parent / "nlcd_calc_tracts.csv"
tracts.drop(columns = ["geometry"]).to_csv(csv_out, index = False)
print("Saved CSV:", csv_out)
Saved: data\raster\nlcd_calc_tracts.geojson
Saved CSV: data\raster\nlcd_calc_tracts.csv
← 4. Additional Features 6. Data Merging →

Notebooks

This Notebook

Source: 05_nlcd_calculations.ipynb

Code Cells: 9

Figures: 0

Data Pipeline

Weather → 311 → Census → Features → Merge → EDA → Modeling