CODE: Additional Feature Engineering

## Module
from pathlib import Path
import geopandas as gpd
import rasterio
from rasterio.mask import mask
import numpy as np
import pandas as pd

## File Paths
tracts_path = Path("data/nyc_tracts_2020/nyc_tracts_2020.shp")
NDVI_dir = Path("data/raster/NDVI.tif")

# 1. Load NYC census tracts shapefile
tracts = gpd.read_file(tracts_path)

# 2. Open NDVI raster and ensure CRS matches the vector layer
with rasterio.open(NDVI_dir) as src:
    raster_crs = src.crs
    nodata = src.nodata

    # Reproject tracts to match raster CRS if necessary
    if tracts.crs != raster_crs:
        tracts = tracts.to_crs(raster_crs)

    ndvi_means = []

    # 3. Compute zonal mean NDVI for each tract
    for _, row in tracts.iterrows():
        geom = [row.geometry]

        # Mask raster using the tract polygon
        out_image, out_transform = mask(src, geom, crop=True)

        data = out_image[0].astype("float32")

        # Replace NoData values with NaN to avoid affecting mean
        if nodata is not None:
            data[data == nodata] = np.nan

        # Compute mean NDVI for this tract
        mean_ndvi = float(np.nanmean(data))
        ndvi_means.append(mean_ndvi)

# 4. Add NDVI results to GeoDataFrame
tracts["NDVI"] = ndvi_means

tracts.head()

## File Paths
water_path = Path("data/Water shp/NYC_water.shp")

water = gpd.read_file(water_path)

tracts = tracts.to_crs(water.crs)

# 1. Compute tract area
tracts["tract_area"] = tracts.geometry.area

# 2. Intersect tracts with water polygons
#    This will create pieces of water polygons clipped by tract boundaries
water_in_tracts = gpd.overlay(
    tracts[["geoid", "geometry"]],
    water[["geometry"]],
    how="intersection"
)

# 3. Compute water area inside each tract
water_in_tracts["water_area"] = water_in_tracts.geometry.area

# 4. Aggregate water area by tract (GEOID)
water_area_by_tract = (
    water_in_tracts
    .groupby("geoid")["water_area"]
    .sum()
)

# 5. Map aggregated water area back to the tracts GeoDataFrame
tracts["water_area"] = tracts["geoid"].map(water_area_by_tract).fillna(0)

# 6. Compute Water Coverage Ratio (WCR = water area / tract area)
tracts["WCR"] = tracts["water_area"] / tracts["tract_area"]

# Check the result
tracts[["geoid", "NDVI", "WCR"]].head()

building_path = Path("data/Buildings/geo_export_10da9e2c-833d-4ba4-9fe2-1f999ac16759.shp")
buildings = gpd.read_file(building_path)

buildings = buildings.to_crs(water.crs)
buildings.crs

# 1. Use height_roo as building height (it is in feet), convert to meters
buildings["bldg_height"] = buildings["height_roo"] * 0.3048

# 2. Compute building footprint area (in m² because CRS is UTM)
buildings["bldg_area"] = buildings.geometry.area

# 3. Compute tract area
tracts["tract_area"] = tracts.geometry.area

# 4. Spatial join — assign buildings to tracts
bldg_in_tracts = gpd.sjoin(
    buildings[["bldg_height", "bldg_area", "geometry"]],
    tracts[["geoid", "geometry"]],
    how="inner",
    predicate="intersects"
)

# 5. Weighted height sum for AH
bldg_in_tracts["height_area"] = (
    bldg_in_tracts["bldg_height"] * bldg_in_tracts["bldg_area"]
)

# 6. Aggregate by tract
bldg_stats = (
    bldg_in_tracts
    .groupby("geoid")
    .agg(
        total_bldg_area=("bldg_area", "sum"),
        total_height_area=("height_area", "sum")
    )
)

# 7. Join results back
tracts = tracts.join(bldg_stats, on="geoid")

tracts["total_bldg_area"] = tracts["total_bldg_area"].fillna(0)

# 8. BD (building density: footprint area ratio)
tracts["BD"] = tracts["total_bldg_area"] / tracts["tract_area"]

# 9. AH (area-weighted building height)
tracts["AH"] = tracts["total_height_area"] / tracts["total_bldg_area"]
tracts.loc[tracts["total_bldg_area"] == 0, "AH"] = np.nan

# Preview
tracts[["geoid", "BD", "AH"]].head()

## Save the reusults
# Select only the variables we need (drop geometry implicitly)
tracts_vars = tracts[["geoid", "BD", "AH", "NDVI", "WCR"]].copy()

tracts_vars.columns = tracts_vars.columns.str.upper()

# Save to CSV (no index column)
tracts_vars.to_csv("data/additional_features/nyc_tracts_new_variables.csv", index=False)

tracts_vars

# Predictor lists.
env_variables = ["TREE_CANOPY_PCT", "IMPERVIOUS_RATIO", "WCR","NDVI"]
acs_variables = ["PCT_BACHELORS_PLUS", "PCT_RENTERS", "PCT_LIMITED_ENGLISH", "MEDIAN_INCOME"]
urban_variables = ["KNN_SUBWAY", "POI_DENSITY", "KNN_PARKS","BD","AH"]