# Gentrification 2014 - 2024
# Michael Minn
# 11 June 2026

library(sf)                        

# Load tracts

tracts = st_read("2020-2024-acs-tracts.geojson")

tracts$Percent_New_Housing = tracts$Percent_2020_2024_Units + tracts$Percent_2010_2019_Units

tracts = tracts[, c("GEOIDFQ", "ST", "Name", "Total_Population",
	"Percent_Bachelors_Degree", "Percent_Jobs_Management", "Median_Household_Income", 
	"Median_Home_Value", "Median_Monthly_Rent", "Percent_New_Housing", "geometry")]

tract_centroids = st_centroid(tracts)

tract_centroids = tract_centroids[tract_centroids$Total_Population > 50,]


# Spatial join to get 2014 data

tracts_2014 = st_read("2010-2014-acs-tracts.geojson")

tracts_2014 = tracts_2014[, c("Percent_Bachelors_Degree", "Percent_Jobs_Management",
"Median_Household_Income", "Median_Home_Value", "Median_Monthly_Rent", "geometry")]

names(tracts_2014) = 
c("Percent_Bachelors_Degree_2014", "Percent_Jobs_Management_2014",
"Median_Household_Income_2014", "Median_Home_Value_2014", "Median_Monthly_Rent_2014", "geometry")

tract_centroids = st_join(tract_centroids, tracts_2014)

tracts = merge(tracts, st_drop_geometry(tract_centroids[, c("GEOIDFQ", "Percent_Bachelors_Degree_2014",
	"Percent_Jobs_Management_2014", "Median_Household_Income_2014", 
	"Median_Home_Value_2014", "Median_Monthly_Rent_2014")]), all.x=T)


# Aggregate to get county values

tracts$County_GEOIDFQ = sub("1400000US", "0500000US", substr(tracts$GEOIDFQ, 1, 14))

tract_centroids = st_centroid(tracts)

x = aggregate(tract_centroids$Percent_Bachelors_Degree_2014, 
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Bachelors_2014")
county_change = x

x = aggregate(tract_centroids$Percent_Jobs_Management_2014, 
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Management_2014")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Household_Income_2014, 
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Income_2014")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Household_Income_2014, 
	list(tract_centroids$County_GEOIDFQ), function(x) quantile(x, probs=0.40, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Income_40pct_2014")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Household_Income_2014, 
	list(tract_centroids$County_GEOIDFQ), function(x) quantile(x, probs=0.75, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Income_75pct_2014")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Percent_Bachelors_Degree - tract_centroids$Percent_Bachelors_Degree_2014,
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Bachelor_Change")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Home_Value - tract_centroids$Median_Home_Value_2014,
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Value_Change")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Home_Value - tract_centroids$Median_Home_Value_2014,
	list(tract_centroids$County_GEOIDFQ), function(x) quantile(x, probs=0.75, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Value_Change_75pct")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Monthly_Rent - tract_centroids$Median_Monthly_Rent_2014,
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Rent_Change")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Median_Monthly_Rent - tract_centroids$Median_Monthly_Rent_2014,
	list(tract_centroids$County_GEOIDFQ), function(x) quantile(x, probs=0.75, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_Rent_Change_75pct")
county_change = merge(county_change, x, all=T)

x = aggregate(tract_centroids$Percent_New_Housing,
	list(tract_centroids$County_GEOIDFQ), function(x) median(x, na.rm=T))
names(x) = c("County_GEOIDFQ", "County_New_Housing")
county_change = merge(county_change, x, all=T)

tracts = merge(tracts, county_change, all.x=T)


# County polygons for display

counties = st_read("2020-2024-acs-counties.geojson")

counties = counties[,c("GEOIDFQ", "ST", "Percent_Jobs_Management",
	"Percent_Bachelors_Degree", "Median_Household_Income", "geometry")]

names(county_change)[1] = "GEOIDFQ"

counties = merge(counties, county_change)

counties$SSI = counties$Percent_Jobs_Management + counties$Percent_Bachelors_Degree + log(counties$Median_Household_Income)

counties$SSI_2014 = counties$County_Management_2014 + counties$County_Bachelors_2014 + 
	log(counties$County_Income_2014)

counties$SSI_Change = counties$SSI - counties$SSI_2014

st_write(counties, "2024-gentrification-counties.geojson", delete_dsn=T)

# Rural_Urban_Core: https://www.ers.usda.gov/data-products/rural-urban-commuting-area-codes
# Metro core, Micro core, Small town core, Rural, Water

urban_core = read.csv("2024-gentrification-urban-tracts.csv")

urban_core$GEOIDFQ = ifelse(nchar(urban_core$TractFIPS23) < 11,
		paste0("1400000US0", urban_core$TractFIPS23),
		paste0("1400000US", urban_core$TractFIPS23))

urban_core = urban_core[,c("GEOIDFQ", "UrbanCoreType")]

names(urban_core) = c("GEOIDFQ", "Rural_Urban_Core")

tracts = merge(tracts, urban_core, all.x=T)


# Gentrification classification

tracts$Rent_Change = tracts$Median_Monthly_Rent - tracts$Median_Monthly_Rent_2014

tracts$Home_Value_Change = tracts$Median_Home_Value - tracts$Median_Home_Value_2014

tracts$Bachelor_Change = tracts$Percent_Bachelors_Degree - tracts$Percent_Bachelors_Degree_2014

# Freeman, Lance. 2005. "Displacement or Succession? Residential Mobility in Gentrifying Neighborhoods."
# Urban Affairs Review 40(4): 463-491. https://doi.org/10.1177/1078087404273341.

tracts$Gentrified_Freeman =
	ifelse ((tracts$Rural_Urban_Core != "Metro core") |
		(tracts$Median_Household_Income_2014 >= quantile(tracts$County_Income_2014, 0.4, na.rm=T)) |
		(tracts$Percent_New_Housing >= quantile(tracts$County_New_Housing, 0.4, na.rm=T)), "Ineligible",
		    ifelse ((tracts$Bachelor_Change < quantile(tracts$County_Bachelor_Change, 0.5, na.rm=T)) |
			   (tracts$Rent_Change < quantile(tracts$County_Rent_Change, 0.5, na.rm=T)), 
			   "Non-Gentrifying", "Gentrifying"))


# NYU Furman Center. 2016. “State of New York City’s Housing and Neighborhoods in 2015.” 
# Updated 9 May. https://furmancenter.org/research/sonychan/2015-report.

tracts$Gentrified_Furman = ifelse(tracts$Median_Household_Income_2014 >= tracts$County_Income_40pct_2014, "Higher-Income",
			ifelse(tracts$Rent_Change >= tracts$County_Rent_Change, "Gentrifying", "Non-Gentrifying"))


# Ding, Lei, Jackelyn Hwang, and Eileen Divringi. 2016.
# "Gentrification and Residential Mobility in Philadelphia." 
# Regional Science and Urban Economics 61: 38-51.
# https://doi.org/10.1016/j.regsciurbeco.2016.09.004.

names(x) = c("County_GEOIDFQ", "County_Income_75pct_2014")

tracts$Gentrified_Ding = ifelse(tracts$Median_Household_Income_2014 >= tracts$County_Income_75pct_2014, "Ineligible",
			ifelse(tracts$Bachelor_Change < tracts$County_Bachelor_Change, "Non-Gentrifying",
			ifelse(tracts$Rent_Change > tracts$County_Rent_Change_75pct, "Intense",
			ifelse(tracts$Rent_Change > tracts$County_Rent_Change, "Gentrifying",
			ifelse(tracts$Home_Value_Change > tracts$County_Value_Change_75pct, "Intense",
			ifelse(tracts$Home_Value_Change > tracts$County_Value_Change, "Gentrifying", "Non-Gentrifying"))))))

# Yonto, Daniel, and Claire Schuch. 2020. 
# “Developing and Ground-Truthing Multi-Scalar Approaches to Mapping Gentrification.”
# Papers in Applied Geography 6 (4): 352–68. https://doi.org/10.1080/23754931.2020.1789499.

tracts$SSI = round(tracts$Percent_Jobs_Management + tracts$Percent_Bachelors_Degree 
	+ log(tracts$Median_Household_Income), 2)

tracts$SSI_2014 = round(tracts$Percent_Jobs_Management_2014 + tracts$Percent_Bachelors_Degree_2014 
	+ log(tracts$Median_Household_Income_2014), 2)

tracts$SSI_Change = tracts$SSI - tracts$SSI_2014

st_write(tracts, "2024-gentrification-tracts.geojson", delete_dsn=T)