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 that combines NYC census tract boundaries with land surface temperature, tree-canopy coverage, and median household income. By analyzing all three measures at the same geographic unit (census tract), the project avoids an unverified geographic join between incompatible spatial units.

Reproducibility: The R code below retrieves the observations directly from the cited public ArcGIS feature service, using an ID-based request to avoid the server’s per-request feature limit. It also saves a local CSV copy so the analysis can be rerun on the same data snapshot. The code uses documented field names rather than hypothetical column names.

# Install once in the RStudio Console if needed:
# install.packages(c("jsonlite", "dplyr", "ggplot2", "knitr", "rmarkdown"))

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("Please install the following R packages first: ",
       paste(not_installed, collapse = ", "))
}

library(dplyr)
library(ggplot2)
knitr::opts_chunk$set(echo = TRUE, fig.width = 8, fig.height = 5,
                      warning = FALSE, message = FALSE)
options(timeout = 180)

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

# Prefer a saved snapshot after the first download. This makes the submission
# self-contained if nyc_heat_tracts_snapshot.csv is included with it.
snapshot_file <- "nyc_heat_tracts_snapshot.csv"

if (file.exists(snapshot_file)) {
  raw_data <- read.csv(snapshot_file, colClasses = c(GEOID = "character"))
} 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("Could not retrieve tract IDs from ArcGIS. Check the service and internet access.")
  }

  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%2Ccanopy2017%2Cmed_income",
      "&returnGeometry=false&f=json"
    )
    response <- jsonlite::fromJSON(request_url, flatten = TRUE)
    if (!is.null(response$error)) {
      stop("ArcGIS query error: ", paste(unlist(response$error), collapse = " "))
    }
    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))
  stopifnot(nrow(raw_data) > 0)
  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("Dataset does not have the expected fields: ",
       paste(setdiff(required_fields, names(raw_data)), collapse = ", "))
}

# Preserve geographic identifiers as text and confirm the interpretation of
# canopy coverage by inspecting its numeric range.
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("The tree canopy field has no usable data.")
if (canopy_max <= 1) {
  heat_data$tree_canopy_pct <- heat_data$canopy_original * 100
  canopy_note <- "The source stores canopy as a 0–1 fraction; multiplied by 100."
} else if (canopy_max <= 100) {
  heat_data$tree_canopy_pct <- heat_data$canopy_original
  canopy_note <- "The source stores canopy on a 0–100 percentage scale."
} else {
  stop("The canopy field exceeds 100; its unit must be checked before analysis.")
}

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)

The imported layer has 2325 records. After removing missing identifiers and duplicate tract identifiers, there are 2325 distinct tracts. There are 2190 tracts with complete, valid measurements for the three principal numerical variables. 135 tracts are excluded from the complete-case descriptive comparisons because at least one required measure is missing or invalid. The source stores canopy on a 0–100 percentage scale.

2 Research Question

Primary research question: Across New York City census tracts represented in the selected dataset, how are tree-canopy coverage and median household income associated with land surface temperature?

Motivation. Urban heat is not distributed evenly within cities. Vegetation may be associated with cooler surfaces, while neighborhood socioeconomic conditions may be related to differences in land cover, built environments, and heat exposure. Studying these relationships may identify geographic patterns relevant to urban planning and environmental equity. This proposal asks about associations, not causal effects.

Scope of inference. The findings will describe the included NYC census tracts and the time periods represented by the available measurements. Because the study uses observational, geographically aggregated data rather than randomly assigned treatments, the analysis cannot establish that planting trees or changing household income would cause a particular temperature change. The results should not be generalized automatically to other cities or seasons.

3 Cases

A case is one 2020 U.S. Census tract within New York City represented in the downloaded ArcGIS layer. Each row corresponds to a tract, not to an individual resident, a household, or a satellite pixel. Census tract boundaries were adapted to the NYC shoreline as described in the source service.

  • Number of records downloaded: 2325.
  • Number of distinct tracts after identifier checks: 2325.
  • Number of complete tracts analyzed for the three numerical measures: 2190.
  • Study area: New York City’s five boroughs, subject to coverage and missingness in the source layer.

The unit of analysis matters: results at the tract level must not be interpreted as individual-level relationships (ecological fallacy).

4 Data Collection

This is a secondary analysis of existing spatial and socioeconomic data. I did not collect the observations myself. The selected hosted layer combines several types of information at a census-tract geographic scale:

  1. Land surface temperature: a remotely sensed surface-temperature summary for geographic areas (the published layer labels its numerical field MEAN as lst_mean). The associated 2024 temperature-input service describes Landsat 9 data for July 2, 2024. The analysis treats the field as a surface-temperature measure, not a measurement of ambient air temperature.
  2. Tree canopy: a preexisting canopy-cover attribute (canopy2017) associated with the tracts. The variable name indicates a 2017 canopy layer; the script checks whether it is recorded as a percentage or proportion before use.
  3. Median household income: a socioeconomic attribute (med_income) included in the compiled layer. The exact ACS reference period for this specific compiled field has not been independently confirmed and will be documented before the final report.
  4. Geography: tract IDs (GEOID) and borough names (BoroName) included in the hosted layer. The service describes the tract geography as based on NYC modifications of 2020 Census boundaries.

The variables are already published in a single feature layer at a common geographic unit. The analysis downloads non-geometric attributes and therefore does not require an additional spatial join for the proposal.

5 Type of Study

Observational, cross-sectional, ecological study. No treatments were assigned and no individuals or neighborhoods were randomly allocated to different amounts of tree cover or income. The recorded associations may be affected by unmeasured factors, including building density, impervious surface, proximity to water, and geographic clustering.

6 Data Source

Primary data layer used for the analysis:

ArcGIS FeatureServer, NYC_HVI_Sampat, layer 1 (tracts_with_canopy), a publicly accessible compiled NYC tract layer with land surface temperature (MEAN), tree canopy (canopy2017), and median household income (med_income). Accessed November 8, 2026.

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

Related temperature metadata:

NYC_HVI_LST_2024_Sampat, layer 13 (HVI_Input_LST), described by its provider as Landsat 9 surface-temperature input for July 2, 2024.

https://services2.arcgis.com/IsDCghZ73NgoYoz5/ArcGIS/rest/services/NYC_HVI_LST_2024_Sampat/FeatureServer/13

Authoritative background and context:

New York City Department of Health and Mental Hygiene. Interactive Heat Vulnerability Index, NYC Environment & Health Data Portal.

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

Data caveat: The ArcGIS feature service is a compiled secondary layer, not the city’s official HVI rankings dataset. The NYC Health Department webpage is cited for background context rather than as the direct host of these analysis records. The code saves a local CSV snapshot for reproducibility. Because the feature service does not display explicit redistribution terms on the metadata page reviewed, verify reuse terms before publicly posting or redistributing that snapshot; the cited live source and fully reproducible code are provided regardless. These data were not taken from the DATA 606 lectures or textbook.

7 Describe Your Variables

Variable Field in source Role Type Unit / Interpretation
Land surface temperature MEAN Response (dependent) Numerical, continuous Surface temperature, approximately degrees Celsius; verify against the original raster metadata before final publication.
Tree-canopy coverage canopy2017 Explanatory 1 Numerical, continuous Percentage of each tract covered by tree canopy (normalized to 0–100 after checking the source scale).
Median household income med_income Explanatory 2 Numerical U.S. dollars per household; used as $10,000 units for regression interpretation.
Borough BoroName Context / optional control Categorical, nominal One of NYC’s five boroughs.
Census tract ID GEOID Identifier Categorical, nominal Unique geographic identifier; not a numeric predictor.

Expected patterns, without assuming the outcome: Tracts with greater tree-canopy coverage may have lower measured surface temperatures. Median income may also be associated with temperature, but the magnitude and direction must be established from the data rather than presumed.

8 Relevant Summary Statistics

The following summaries are calculated from the actual downloaded records, not hypothetical values. Each numerical variable is described using the nonmissing count, missing count, mean, median, standard deviation, quartiles, and observed range. Descriptive figures then show the shape of the temperature distribution and its association with the two explanatory measures.

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]]
  # For canopy, use the normalized percentage variable in heat_data.
  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 key concern in this section is whether the distributions show substantial variation and potential outliers. The median, quartiles, and histogram help identify skewness, while the scatterplots show the direction, form, and strength of the two bivariate relationships.

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

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

ggplot(analysis_data, aes(x = median_income_usd, y = surface_temp_c)) +
  geom_point(alpha = 0.45, color = "#B46731") +
  geom_smooth(method = "lm", 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; fitted line is descriptive",
    x = "Median household income (USD)",
    y = "Land surface temperature (approximately degrees C)"
  ) +
  theme_minimal(base_size = 12)

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 degrees C)"
  ) +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

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 surface temperature")
Exploratory Pearson correlations with surface temperature
Explanatory_variable Pearson_r
Tree-canopy coverage (%) -0.239
Median household income (USD) -0.132

Initial interpretation (to be finalized from the knitted output): The table and plots above provide a direct descriptive assessment of central tendency, spread, and potential relationships. A negative correlation would indicate that higher levels of the explanatory variable tend to occur in cooler tracts; a positive correlation would indicate the opposite. Correlation and a fitted trendline alone do not prove causation. The specific direction and magnitude of each relationship should be reported after the code has run successfully.

9 Planned Analysis for the Final Project

I plan to use multiple linear regression to examine how canopy coverage and median household income are jointly associated with tract-level land surface temperature:

\[ \text{LST}_i = \beta_0 + \beta_1\text{CanopyPct}_i + \beta_2\text{MedianIncome10k}_i + \varepsilon_i. \]

The estimated canopy coefficient would express the expected temperature difference associated with a one-percentage-point increase in canopy coverage, holding median income constant. The income coefficient would express the corresponding expected temperature difference per $10,000 of median household income, holding canopy coverage constant. I will report the estimated coefficients, their uncertainty, p-values where appropriate, and the model’s \(R^2\), interpreted as the proportion of observed variation in tract-level LST explained by the model.

I will inspect scatterplots and residuals for nonlinearity, unequal variance, influential observations, and possible spatial autocorrelation. If clustering is substantial, I will caution that ordinary least-squares standard errors may be too optimistic and consider appropriate robustness checks. Borough can be added as a categorical control in a secondary analysis.

10 Significance and Limitations

The study may contribute to understanding geographic inequalities in heat exposure and how environmental and socioeconomic characteristics vary alongside surface temperature. Such findings could be useful in discussions of urban greening, climate adaptation, and environmental justice.

However, the study has important limitations: (1) it is observational and ecological, not causal; (2) surface temperature is different from air temperature and human heat exposure; (3) canopy, temperature, and income observations may refer to different years, so the analysis is not a synchronized experiment; (4) missing measurements may affect which tracts are analyzed; (5) geographic proximity can produce spatial autocorrelation; and (6) the compiled data source’s income vintage and temperature processing metadata need further documentation before the final report. The proposal therefore focuses on defensible descriptive and associational findings, not causal claims.

11 References