# set path
setwd()


###############################################################################
# 1) Load data
###############################################################################
dturban <- read.csv("dturban.csv")


###############################################################################
# 2) Calculate geocomplexity (GC) for selected variables
#    - Build an sf object from lon/lat
#    - Create spatial weights (contiguity-based)
#    - Compute GC using Moran-based method
###############################################################################
install.packages('geocomplexity',
                 repos = c("https://ausgis.r-universe.dev",
                           "https://cloud.r-project.org"),
                 dep = TRUE)
library(sf)
library(geocomplexity)

d1 <- dturban[, c(2, 3, 5:9)]
sf.data <- st_as_sf(d1, coords = c("lon", "lat"), crs = 4326)

wt <- sdsfun::spdep_contiguity_swm(sf_data, k = 44)
vars.gc <- geocd_vector(sf.data, wt, method = "moran", normalize = TRUE)

# Bind GC results back to the original data and keep the first 14 columns
dturban.gc <- cbind(dturban, vars.gc[, 1:5])
dturban.gc <- dturban.gc[, 1:14]


###############################################################################
# 3) Plot relationship: variable vs its GC pattern
#    - Scatter plot + GAM smoother
###############################################################################
library(ggplot2)

temp_data <- data.frame(
  x = dturban.gc$agriculture,
  y = dturban.gc$GC_agriculture
)

ggplot(data = temp_data, aes(x = x, y = y)) +
  geom_point(color = "darkgray", size = 0.5) +
  geom_smooth(method = "gam", formula = y ~ s(x)) +
  theme_bw()


###############################################################################
# 4) GD (Geodetector) for factor analysis
#    - Select variables
#    - Discretize continuous variables with multiple methods and interval counts
#    - Run gdm and visualize
###############################################################################
library(GD)

dt.gd <- dturban.gc[, 4:14]
discmethod <- c("equal", "quantile", "natural")
discitv <- 4:5
contvar <- names(dt.gd)[-1]

gd.result <- gdm(
  formula = urbanindex ~ .,
  continuous_variable = contvar,
  data = dt.gd,
  discmethod = discmethod,
  discitv = discitv
)

gd.result
plot(gd.result)

# Extract interaction detector results
gd.result.interaction <- gd.result$Interaction.detector$Interaction


###############################################################################
# 5) GOZH analysis (decision-tree-based partition + GD)
#    - First: GOZH for variables only (dturban)
#    - Second: GOZH for variables + GC patterns (dturban.gc)
###############################################################################
library(GD)
library(rpart)
library(rpart.plot)

# ---- GOZH for five variables ----
dt.vars <- dturban[, 4:9]
tree <- rpart(urbanindex ~ ., data = dt.vars)

# Create partition label from tree node assignment
dt.vars$tree <- as.character(as.numeric(tree$where))

# Run GD on the partition
q.vars <- gd(urbanindex ~ tree, dt.vars)
q.vars

# Visualize the decision tree
rpart.plot(tree, digits = 3, box.palette = "GnYlRd")


# ---- GOZH for variables and their geocomplexity patterns ----
dt.vars.gc <- dturban.gc[, 4:14]
tree <- rpart(urbanindex ~ ., data = dt.vars.gc)

# Create partition label from tree node assignment
dt.vars.gc$tree <- as.character(as.numeric(tree$where))

# Run GD on the partition
q.vars.gc <- gd(urbanindex ~ tree, dt.vars.gc)
q.vars.gc

# Visualize the decision tree
rpart.plot(tree, digits = 3, box.palette = "GnYlRd")


###############################################################################
# 6) Improvement metric: relative gain in q-value after adding GC patterns
###############################################################################
(q.vars.gc$Factor$qv - q.vars$Factor$qv) / q.vars$Factor$qv


###############################################################################
# 7) Create maps of variables and GC patterns
#    - Read study area shapefile (with PID)
#    - Join with dturban.gc by PID
#    - Plot: agriculture and GC_agriculture using different color schemes
###############################################################################
library(sf)
library(dplyr)
library(ggplot2)
library(patchwork)
library(scales)

# Read shapefile and standardize PID type
sa <- st_read("studyarea.shp", quiet = TRUE) %>%
  mutate(PID = as.character(PID))

# Standardize PID type in table
dturban.gc <- dturban.gc %>%
  mutate(PID = as.character(PID))

# Join attributes to polygons
sa2 <- sa %>%
  left_join(dturban.gc %>% select(PID, agriculture, GC_agriculture), by = "PID")

# --- Map 1: agriculture (YlGn) ---
ggplot(sa2) +
  geom_sf(aes(fill = agriculture), color = "gray") +
  scale_fill_distiller(palette = "YlGn", direction = 1) +
  theme_bw()

# --- Map 2: GC_agriculture (turbo) ---
ggplot(sa2) +
  geom_sf(aes(fill = GC_agriculture), color = "gray") +
  scale_fill_viridis_c(option = "turbo") +
  theme_bw()



