Urban Heat Inequality in New York City: Tree Canopy, Neighborhood Income, and Land Surface Temperature

Author

Sarah Abdelrahman

Published

November 9, 2026

1 Data Preparation

This proposal examines neighborhood-level variation in land surface temperature (LST) across New York City (NYC). The analysis uses a publicly accessible ArcGIS feature layer containing census tract identifiers, land surface temperature, tree-canopy coverage, and median household income.

Using a single compiled geographic layer avoids an unverified spatial join between incompatible geographic units.

Reproducibility: The R code downloads the records from the cited ArcGIS FeatureServer using object-ID batches, avoiding the service’s per-request record limit. A local CSV snapshot is saved to support reproducible analysis.

required_packages <- c(
  "jsonlite",
  "dplyr",
  "ggplot2",
  "knitr"
)

not_installed <- required_packages[
  !vapply(
    required_packages,
    requireNamespace,
    logical(1),
    quietly = TRUE
  )
]

if (length(not_installed) > 0) {
  stop(
    "Install these packages first: ",
    paste(not_installed, collapse = ", ")
  )
}

library(dplyr)
library(ggplot2)

options(timeout = 180)

1.1 Download the Data

The primary data source is the NYC_HVI_Sampat ArcGIS FeatureServer, layer 1.

service_url <- paste0(
  "https://services2.arcgis.com/",
  "IsDCghZ73NgoYoz5/ArcGIS/rest/services/",
  "NYC_HVI_Sampat/FeatureServer/1/query"
)

snapshot_file <- "nyc_heat_tracts_snapshot.csv"

if (file.exists(snapshot_file)) {

  raw_data <- read.csv(
    snapshot_file,
    colClasses = c(GEOID = "character"),
    check.names = FALSE
  )

} else {

  ids_url <- paste0(
    service_url,
    "?where=1%3D1&returnIdsOnly=true&f=json"
  )

  id_response <- jsonlite::fromJSON(ids_url)

  if (
    !is.null(id_response$error) ||
    length(id_response$objectIds) == 0
  ) {
    stop("Unable to retrieve ArcGIS object IDs.")
  }

  object_ids <- sort(unique(id_response$objectIds))

  batches <- split(
    object_ids,
    ceiling(seq_along(object_ids) / 150)
  )

  get_batch <- function(id_batch) {

    request_url <- paste0(
      service_url,
      "?objectIds=",
      paste(id_batch, collapse = ","),
      "&outFields=",
      "OBJECTID%2CGEOID%2CBoroName%2CMEAN%2C",
      "canopy2017%2Cmed_income",
      "&returnGeometry=false&f=json"
    )

    response <- jsonlite::fromJSON(
      request_url,
      flatten = TRUE
    )

    if (!is.null(response$error)) {
      stop("ArcGIS returned a query error.")
    }

    if (
      is.null(response$features) ||
      nrow(response$features) == 0
    ) {
      return(data.frame())
    }

    result <- response$features

    names(result) <- sub(
      "^attributes\\.",
      "",
      names(result)
    )

    result
  }

  raw_data <- dplyr::bind_rows(
    lapply(batches, get_batch)
  )

  if (nrow(raw_data) == 0) {
    stop("No data were downloaded.")
  }

  write.csv(
    raw_data,
    snapshot_file,
    row.names = FALSE,
    na = ""
  )
}

required_fields <- c(
  "GEOID",
  "BoroName",
  "MEAN",
  "canopy2017",
  "med_income"
)

if (!all(required_fields %in% names(raw_data))) {
  stop("Required fields are missing from the dataset.")
}

1.2 Clean and Prepare the Data

The analysis preserves census tract identifiers as text, checks tree-canopy values, and excludes records with missing or invalid measurements.

heat_data <- raw_data %>%
  transmute(
    tract_id = as.character(GEOID),
    borough = as.factor(BoroName),
    surface_temp_c = as.numeric(MEAN),
    canopy_original = as.numeric(canopy2017),
    median_income_usd = as.numeric(med_income)
  ) %>%
  filter(
    !is.na(tract_id),
    tract_id != ""
  ) %>%
  distinct(
    tract_id,
    .keep_all = TRUE
  )

canopy_max <- max(
  heat_data$canopy_original,
  na.rm = TRUE
)

if (!is.finite(canopy_max)) {
  stop("Tree-canopy values are unavailable.")
}

if (canopy_max <= 1) {

  heat_data$tree_canopy_pct <-
    heat_data$canopy_original * 100

  canopy_note <- "Canopy values converted from proportions to percentages."

} else if (canopy_max <= 100) {

  heat_data$tree_canopy_pct <-
    heat_data$canopy_original

  canopy_note <- "Canopy values already use a percentage scale."

} else {

  stop("Tree-canopy units require verification.")
}

analysis_data <- heat_data %>%
  filter(
    is.finite(surface_temp_c),
    is.finite(tree_canopy_pct),
    tree_canopy_pct >= 0,
    tree_canopy_pct <= 100,
    is.finite(median_income_usd),
    median_income_usd > 0
  ) %>%
  mutate(
    median_income_10k = median_income_usd / 10000
  )

stopifnot(
  nrow(analysis_data) >= 20,
  dplyr::n_distinct(analysis_data$tract_id) ==
    nrow(analysis_data),
  !anyNA(analysis_data$tract_id)
)

1.3 Data Quality Checks

quality_table <- data.frame(
  Measure = c(
    "Downloaded records",
    "Distinct census tracts",
    "Complete analytical records",
    "Excluded records",
    "Duplicate tract IDs in analytical sample",
    "Missing tract IDs in analytical sample"
  ),
  Value = c(
    nrow(raw_data),
    nrow(heat_data),
    nrow(analysis_data),
    nrow(heat_data) - nrow(analysis_data),
    sum(duplicated(analysis_data$tract_id)),
    sum(is.na(analysis_data$tract_id))
  )
)

knitr::kable(
  quality_table,
  caption = "Data quality and analytical sample"
)
Data quality and analytical sample
Measure Value
Downloaded records 2325
Distinct census tracts 2325
Complete analytical records 2190
Excluded records 135
Duplicate tract IDs in analytical sample 0
Missing tract IDs in analytical sample 0

The downloaded dataset contains 2325 records. After cleaning, 2190 unique census tracts remain for analysis. 135 records are excluded because required measurements are missing or invalid.

Canopy values already use a percentage scale.

2 Research Question

Primary research question:

Across New York City census tracts, how are tree-canopy coverage and median household income associated with land surface temperature?

2.1 Motivation

Urban heat is not distributed equally across neighborhoods. Differences in vegetation, socioeconomic conditions, land cover, and the built environment may contribute to geographic patterns in surface temperature.

Tree canopy may be associated with cooler urban surfaces through shading and evapotranspiration. Median household income may also be related to temperature through neighborhood characteristics and historical patterns of urban development.

This study investigates these relationships using tract-level data.

The research examines associations, not causal effects.

2.2 Scope of Inference

Findings will describe the NYC census tracts represented in the selected dataset and the time periods covered by its measurements.

Because the study uses observational data, its results cannot establish that increasing tree cover or household income would directly cause a change in temperature.

3 Cases

The observational unit is a 2020 U.S. Census tract within New York City represented in the selected ArcGIS layer.

Each record represents one census tract rather than an individual resident, household, or satellite pixel.

  • Original records: 2325
  • Distinct tracts before complete-case filtering: 2325
  • Complete analytical sample: 2190
  • Geographic scope: New York City’s five boroughs, subject to source coverage and missingness.

The analysis is conducted at the census-tract level. Results must not be interpreted as relationships between individual people, because doing so could lead to an ecological fallacy.

4 Data Collection

This project uses secondary data analysis. No original survey, experiment, or field measurements were conducted by the researcher.

The hosted ArcGIS feature layer combines several types of information.

4.1 Land Surface Temperature

Land surface temperature is represented by the source field MEAN, identified in the compiled layer as a surface-temperature summary.

Related temperature metadata describe Landsat 9 imagery from July 2, 2024.

Land surface temperature measures the temperature of surfaces, not ambient air temperature.

The values are treated as approximately degrees Celsius, pending confirmation from the original raster metadata.

4.2 Tree-Canopy Coverage

Tree-canopy coverage comes from the field canopy2017.

The variable name indicates a 2017 canopy dataset. The R script checks the numerical scale and expresses canopy coverage as a percentage between 0 and 100.

4.3 Median Household Income

Median household income comes from the source field med_income.

It represents neighborhood socioeconomic conditions and is measured in U.S. dollars.

The exact American Community Survey reference period for this compiled field has not been independently verified.

4.4 Geographic Information

The source includes:

  • GEOID: Census tract identifier.
  • BoroName: New York City borough.

All primary measures are included in the same compiled tract-level layer.

5 Type of Study

This project is an observational, cross-sectional, ecological study.

No treatments are assigned, and census tracts are not randomly allocated to different levels of tree canopy or income.

The research analyzes existing measurements to identify statistical relationships.

Potential confounding variables include building density, impervious surface coverage, proximity to water, neighborhood development patterns, and geographic clustering.

Therefore, the findings will be interpreted as associations rather than causal effects.

6 Data Source

6.1 Primary Dataset

ArcGIS FeatureServer: NYC_HVI_Sampat, Layer 1 (tracts_with_canopy)

https://services2.arcgis.com/IsDCghZ73NgoYoz5/ArcGIS/rest/services/NYC_HVI_Sampat/FeatureServer/1

The dataset provides tract identifiers, borough names, land surface temperature, tree-canopy coverage, and median household income.

Accessed November 8, 2026.

6.3 Background Reference

New York City Department of Health and Mental Hygiene — Interactive Heat Vulnerability Index

https://a816-dohbesp.nyc.gov/IndicatorPublic/data-features/hvi/

The NYC Health Department source provides background context about urban heat vulnerability. It is not the direct source of the compiled analysis records.

Data limitation: The ArcGIS layer is a compiled secondary dataset rather than the official NYC Health Department HVI rankings dataset. The exact income reference period and temperature-processing metadata require further verification before the final project.

7 Describe Your Variables

7.1 Response Variable

Land Surface Temperature (surface_temp_c)

  • Source field: MEAN
  • Role: Response/dependent variable
  • Type: Numerical, continuous
  • Unit: Approximately degrees Celsius, pending source verification

The response variable represents tract-level surface temperature.

7.2 Explanatory Variable 1

Tree-Canopy Coverage (tree_canopy_pct)

  • Source field: canopy2017
  • Role: Explanatory/independent variable
  • Type: Numerical, continuous
  • Unit: Percentage (%)

This variable measures the percentage of a census tract covered by tree canopy.

7.3 Explanatory Variable 2

Median Household Income (median_income_usd)

  • Source field: med_income
  • Role: Explanatory/independent variable
  • Type: Numerical
  • Unit: U.S. dollars

For regression interpretation, income will also be expressed in $10,000 units.

7.4 Additional Variables

Borough (borough)

  • Source field: BoroName
  • Type: Categorical, nominal
  • Purpose: Geographic comparison and possible regression control

Census Tract ID (tract_id)

  • Source field: GEOID
  • Type: Categorical identifier
  • Purpose: Uniquely identifies each observation

7.5 Expected Relationships

The study will investigate whether greater tree-canopy coverage is associated with lower land surface temperature.

It will also evaluate whether median household income is related to temperature.

The direction and magnitude of the relationships will be determined from the data rather than assumed.

8 Relevant Summary Statistics

Descriptive statistics are calculated from the downloaded records.

For each primary numerical variable, the analysis reports the number of valid observations, missing values, mean, median, standard deviation, minimum, maximum, and quartiles.

8.1 Descriptive Statistics

variable_labels <- c(
  surface_temp_c = "Land surface temperature (C)",
  tree_canopy_pct = "Tree-canopy coverage (%)",
  median_income_usd = "Median household income (USD)"
)

num_summary <- dplyr::bind_rows(
  lapply(names(variable_labels), function(v) {

    x <- heat_data[[v]]

    ok <- is.finite(x)
    y <- x[ok]

    data.frame(
      Variable = unname(variable_labels[v]),
      N = sum(ok),
      Missing = sum(!ok),
      Mean = round(mean(y), 2),
      Median = round(median(y), 2),
      SD = round(sd(y), 2),
      Min = round(min(y), 2),
      Q1 = round(
        unname(quantile(y, 0.25)), 2
      ),
      Q3 = round(
        unname(quantile(y, 0.75)), 2
      ),
      Max = round(max(y), 2),
      check.names = FALSE
    )
  })
)

knitr::kable(
  num_summary,
  caption = "Descriptive statistics for primary study variables"
)
Descriptive statistics for primary study variables
Variable N Missing Mean Median SD Min Q1 Q3 Max
Land surface temperature (C) 2325 0 40.77 40.99 2.26 23.28 39.81 42.20 46.93
Tree-canopy coverage (%) 2317 8 18.89 18.37 7.08 2.43 14.67 21.81 50.47
Median household income (USD) 2192 133 88549.24 81981.50 42443.02 11612.00 60240.75 107326.75 250001.00

The descriptive statistics show variation in temperature, tree-canopy coverage, and household income across the included NYC census tracts.

Comparing the mean and median helps assess possible skewness, while quartiles and ranges provide information about dispersion and potential outliers.

8.2 Distribution of Land Surface Temperature

ggplot(
  analysis_data,
  aes(x = surface_temp_c)
) +
  geom_histogram(
    bins = 30,
    fill = "#2C5F78",
    color = "white"
  ) +
  labs(
    title = "Distribution of land surface temperature",
    subtitle = paste(
      "Complete-case census tracts:",
      nrow(analysis_data)
    ),
    x = "Land surface temperature (approximately C)",
    y = "Number of census tracts"
  ) +
  theme_minimal(base_size = 12)
Figure 1: Distribution of land surface temperature across NYC census tracts

The histogram illustrates the distribution of tract-level land surface temperature and helps identify its central concentration, spread, and potential extreme observations.

8.3 Tree Canopy and Land Surface Temperature

ggplot(
  analysis_data,
  aes(
    x = tree_canopy_pct,
    y = surface_temp_c
  )
) +
  geom_point(
    alpha = 0.45,
    color = "#23744D"
  ) +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = TRUE,
    color = "#153F63"
  ) +
  labs(
    title = "Tree-canopy coverage and land surface temperature",
    subtitle = "Each point represents one NYC census tract",
    x = "Tree-canopy coverage (%)",
    y = "Land surface temperature (approximately C)"
  ) +
  theme_minimal(base_size = 12)
Figure 2: Tree-canopy coverage and land surface temperature

The scatterplot examines the direction and approximate strength of the relationship between tree-canopy coverage and surface temperature.

The fitted line summarizes the observed bivariate association and should not be interpreted as causal evidence.

8.4 Median Household Income and Temperature

ggplot(
  analysis_data,
  aes(
    x = median_income_usd,
    y = surface_temp_c
  )
) +
  geom_point(
    alpha = 0.45,
    color = "#B46731"
  ) +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = TRUE,
    color = "#153F63"
  ) +
  scale_x_continuous(
    labels = function(x) {
      paste0(
        "$",
        format(
          x,
          big.mark = ",",
          trim = TRUE
        )
      )
    }
  ) +
  labs(
    title = "Median household income and land surface temperature",
    subtitle = "Each point represents one NYC census tract",
    x = "Median household income (USD)",
    y = "Land surface temperature (approximately C)"
  ) +
  theme_minimal(base_size = 12)
Figure 3: Median household income and land surface temperature

This scatterplot explores whether census tracts with different median household incomes tend to exhibit different surface temperatures.

8.5 Temperature Variation by Borough

ggplot(
  analysis_data,
  aes(
    x = borough,
    y = surface_temp_c
  )
) +
  geom_boxplot(
    fill = "#A9C3D1",
    outlier.alpha = 0.4
  ) +
  labs(
    title = "Variation in land surface temperature by borough",
    x = "Borough",
    y = "Land surface temperature (approximately C)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(
      angle = 20,
      hjust = 1
    )
  )
Figure 4: Land surface temperature by NYC borough

The boxplot compares temperature distributions across the five boroughs and helps identify geographic differences in medians, variability, and outliers.

8.6 Correlation Analysis

Pearson correlation coefficients are calculated to summarize the linear associations between land surface temperature and the two numerical explanatory variables.

correlation_table <- data.frame(
  Explanatory_variable = c(
    "Tree-canopy coverage (%)",
    "Median household income (USD)"
  ),
  Pearson_r = c(
    cor(
      analysis_data$surface_temp_c,
      analysis_data$tree_canopy_pct
    ),
    cor(
      analysis_data$surface_temp_c,
      analysis_data$median_income_usd
    )
  )
)

correlation_table$Pearson_r <- round(
  correlation_table$Pearson_r,
  3
)

knitr::kable(
  correlation_table,
  caption = "Exploratory Pearson correlations with land surface temperature"
)
Exploratory Pearson correlations with land surface temperature
Explanatory_variable Pearson_r
Tree-canopy coverage (%) -0.239
Median household income (USD) -0.132

8.7 Interpretation of Preliminary Findings

The exploratory analysis identifies two negative linear associations:

  • Tree-canopy coverage: Pearson’s r = -0.239.
  • Median household income: Pearson’s r = -0.132.

The negative correlation between tree-canopy coverage and land surface temperature suggests that census tracts with more tree canopy tend to have lower surface temperatures. The relationship is relatively weak in magnitude.

The negative correlation between median household income and land surface temperature suggests that higher-income tracts also tend to have somewhat lower surface temperatures. This association is weaker than the tree-canopy association in the observed data.

These are descriptive, unadjusted correlations. They do not establish causation and do not show whether either relationship remains after controlling for the other explanatory variable.

The final project will use multiple regression to investigate these relationships jointly.

9 Planned Analysis for the Final Project

9.1 Multiple Linear Regression

The final project will use multiple linear regression to examine how tree-canopy coverage and median household income are jointly associated with census tract land surface temperature.

The proposed model is:

\[ LST_i = \beta_0 + \beta_1 CanopyPct_i + \beta_2 MedianIncome10k_i + \varepsilon_i \]

Where:

  • \(LST_i\) represents land surface temperature in census tract \(i\).
  • \(CanopyPct_i\) represents tree-canopy coverage as a percentage.
  • \(MedianIncome10k_i\) represents median household income in $10,000 units.
  • \(\beta_0\) is the intercept.
  • \(\beta_1\) is the tree-canopy regression coefficient.
  • \(\beta_2\) is the income regression coefficient.
  • \(\varepsilon_i\) is the model error term.

9.2 Planned Statistical Inference

The final analysis will test whether each explanatory variable has a nonzero regression coefficient after accounting for the other variable.

For tree-canopy coverage:

\[ H_0: \beta_1 = 0 \]

\[ H_A: \beta_1 \neq 0 \]

For median household income:

\[ H_0: \beta_2 = 0 \]

\[ H_A: \beta_2 \neq 0 \]

The final report will include coefficient estimates, standard errors, test statistics, p-values, and confidence intervals, where appropriate.

9.3 Model Interpretation

The estimated tree-canopy coefficient will represent the expected difference in land surface temperature associated with a one-percentage-point increase in tree canopy, holding median household income constant.

The income coefficient will represent the expected temperature difference associated with a $10,000 increase in median household income, holding tree canopy constant.

The coefficient of determination, \(R^2\), will be interpreted as the proportion of observed variation in tract-level temperature explained by the regression model.

9.4 Model Diagnostics

The final project will evaluate:

  • Linearity of the modeled relationships.
  • Residual patterns and unequal variance.
  • Influential observations and outliers.
  • Potential multicollinearity.
  • Geographic clustering and spatial autocorrelation.

If substantial spatial autocorrelation is present, ordinary least-squares standard errors may be misleading. This limitation will be discussed, and appropriate robustness checks will be considered.

Borough may also be included as a categorical control variable in a secondary model.

10 Significance and Limitations

10.1 Significance

This project may contribute to understanding geographic inequalities in urban heat exposure and their associations with environmental and socioeconomic conditions.

The findings may help inform discussions about urban greening, climate adaptation, neighborhood planning, and environmental equity.

Combining GIS-based geographic information with statistical analysis provides a framework for studying spatial patterns in environmental conditions.

10.2 Limitations

Several limitations must be considered.

  1. Observational design: The study cannot establish causal relationships.
  2. Ecological analysis: Census tract relationships do not necessarily describe individual residents.
  3. Surface versus air temperature: Land surface temperature is not equivalent to ambient air temperature or personal heat exposure.
  4. Different measurement periods: Temperature, tree canopy, and household income may represent different years.
  5. Missing observations: Excluding incomplete tracts may affect the analytical sample.
  6. Spatial dependence: Nearby census tracts may have correlated measurements.
  7. Compiled data source: The exact income reference period and original temperature-processing details require further verification.

These limitations will be considered when interpreting the final regression results.

11 References

  1. ArcGIS. NYC_HVI_Sampat, FeatureServer Layer 1 (tracts_with_canopy). Accessed November 8, 2026.
    https://services2.arcgis.com/IsDCghZ73NgoYoz5/ArcGIS/rest/services/NYC_HVI_Sampat/FeatureServer/1

  2. ArcGIS. NYC_HVI_LST_2024_Sampat, FeatureServer Layer 13 (HVI_Input_LST). Accessed November 8, 2026.
    https://services2.arcgis.com/IsDCghZ73NgoYoz5/ArcGIS/rest/services/NYC_HVI_LST_2024_Sampat/FeatureServer/13

  3. New York City Department of Health and Mental Hygiene. Interactive Heat Vulnerability Index. Accessed November 8, 2026.
    https://a816-dohbesp.nyc.gov/IndicatorPublic/data-features/hvi/

  4. Bryer, Jason. DATA 606 Fall 2026: Project Instructions.
    https://fall2026.data606.net/assignments/project.html