# Moran's I examples
# Michael Minn
# 15 December 2020

library(sf)

library(spdep)

# Create sample polygon grid

polygons = list()

for (x in seq(-5, 4)) {
	for (y in seq(-5, 4)) {
		polygons = c(polygons, list(matrix(c(x, y, x, y+1, x+1, y+1, x+1, y, x, y), ncol=2, byrow=T)))
	} # for y
} # for x

polygons = lapply(polygons, function(x) { st_polygon(list(x)) })

polygons = st_sfc(polygons, crs=4326)

polygons = st_sf(polygons)

polygons = st_transform(polygons, 3857)


# https://mgimond.github.io/Spatial/spatial-autocorrelation-in-r.html

png("2020-morans-i-theoretical.png", width=700, height=500, pointsize=16)

layout(matrix(c(1,1,3,3,5,5,0,2,2,4,4,0), nrow=2, byrow=T))

par(mar=c(1,1,1,1))

polygons$value = rep(c(rep(c(0,1), 5), rep(c(1,0), 5)), 5)

neighbors = poly2nb(as_Spatial(polygons), queen=F)
weights = nb2listw(neighbors, style="W", zero.policy=T)
m = moran.test(polygons$value, weights)

title = paste0("Perfect Dispersion (I = ", round(m$estimate[1], 2), ")")
plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F)


polygons$value = rep(c(0,0,1,1,1,0,0), 15)[1:100]

neighbors = poly2nb(as_Spatial(polygons), queen=T)
weights = nb2listw(neighbors, style="W", zero.policy=T)
m = moran.test(polygons$value, weights)

title = paste0("Dispersed (I = ", round(m$estimate[1], 2), ")")
plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F)


set.seed(100)
polygons$value = round(runif(100))

neighbors = poly2nb(as_Spatial(polygons), queen=T)
weights = nb2listw(neighbors, style="W", zero.policy=T)
m = moran.test(polygons$value, weights)

title = paste0("Spatial Independence (I = ", round(m$estimate[1], 2), ")")
plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F)


polygons$value = c(rep(c(1,1,1,1,0,0,0,0,0,0), 3),
		   rep(c(0,0,0,0,0,1,1,1,1,0), 3),
		   rep(c(1,1,1,1,0,0,0,0,0,0), 4))

neighbors = poly2nb(as_Spatial(polygons), queen=T)
weights = nb2listw(neighbors, style="W", zero.policy=T)
m = moran.test(polygons$value, weights)

title = paste0("Clustering (I = ", round(m$estimate[1], 2), ")")
plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F)


polygons$value = c(rep(1, 50), rep(0,50))

neighbors = poly2nb(as_Spatial(polygons), queen=F)
weights = nb2listw(neighbors, style="W", zero.policy=T)
m = moran.test(polygons$value, weights)

title = paste0("Perfect Bifurcation (I = ", round(m$estimate[1], 2), ")")
plot(polygons, col=ifelse(polygons$value == 0, "#f0f0f0", "#a10022"), border="#a0a0a0", main=title, reset=F)

dev.off()



counties = st_read("2014-2018-acs-counties.geojson", stringsAsFactors = F)

counties = counties[(counties$FIPS >= "17000") & (counties$FIPS <= "17999"),]

neighbors = poly2nb(as_Spatial(counties), queen=T)

weights = nb2listw(neighbors, style="W", zero.policy=T)

#m = moran.test(counties$Median.Household.Income, weights)
#m

png("2020-morans-i-mortgage.png", width=400, height=600, pointsize=16)
plot(counties[,c("Median.Monthly.Mortgage", "geometry")],
	pal=colorRampPalette(c("#004d3c", "#f0f0f0", "#a10022"))(5),
	nbreaks=5, breaks="jenks", main="", key.pos=1)
m = moran.test(counties$Median.Monthly.Mortgage, weights)
print(m)
dev.off()

png("2020-morans-i-work-home.png", width=400, height=600, pointsize=16)
plot(counties[,c("X..Work.at.Home", "geometry")],
	pal=colorRampPalette(c("#004d3c", "#f0f0f0", "#a10022"))(5),
	nbreaks=5, breaks="jenks", main="", key.pos=1)
m = moran.test(counties$X..Work.at.Home, weights)
print(m)
dev.off()

png("2020-morans-i-income.png", width=400, height=600, pointsize=16)
plot(counties[,c("Median.Household.Income", "geometry")],
	pal=colorRampPalette(c("#004d3c", "#f0f0f0", "#a10022"))(5),
	nbreaks=5, breaks="jenks", main="", key.pos=1)
m = moran.test(counties$Median.Household.Income, weights)
print(m)
dev.off()

stop("Done")





# https://stats.idre.ucla.edu/r/faq/how-can-i-calculate-morans-i-in-r/

points = list()

for (x in seq(-5, 4)) {
	for (y in seq(-5, 4)) {
		polygons = c(polygons, list(matrix(c(x, y, x, y+1, x+1, y+1, x+1, y, x, y), ncol=2, byrow=T)))
	} # for y
} # for x

polygons = lapply(polygons, function(x) { st_polygon(list(x)) })

polygons = st_sfc(polygons, crs=4326)

polygons = st_sf(polygons)

polygons = st_transform(polygons, 3857)
library(ape)

#distances = st_distance(polygons, by_element=F)
#distances = 1 / distances
#diag(distances) = 0
#distances = drop_units(distances)
#Moran.I(polygons$value, distances)