analysis part 3

library(sf)

## Warning: package 'sf' was built under R version 4.5.3

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

stockton <- st_read(
  "C:/Users/noahs/Downloads/ES252/stockton.shp"
)

## Reading layer `stockton' from data source 
##   `C:\Users\noahs\Downloads\ES252\stockton.shp' using driver `ESRI Shapefile'
## Simple feature collection with 86 features and 38 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -121.4197 ymin: 37.86748 xmax: -121.1839 ymax: 38.05876
## Geodetic CRS:  NAD83

names(stockton)

##  [1] "tl_2020_06" "tl_2020_07" "tl_2020_08" "tl_2020_09" "tl_2020_10"
##  [6] "tl_2020_11" "tl_2020_12" "tl_2020_13" "tl_2020_14" "tl_2020_15"
## [11] "tl_2020_16" "tl_2020_17" "Sheet1__Ma" "Sheet1__GE" "Sheet1__F2"
## [16] "Sheet1__Ob" "Shape_Leng" "Shape_Area" "LST_trend_" "LST_pvalue"
## [21] "NDVI_trend" "NDVI_pvalu" "NDBI_trend" "NDBI_pvalu" "Impervious"
## [26] "Impervio_1" "Mean_Albed" "GEOID"      "LST_trend1" "LST_pval_1"
## [31] "NDVI_tre_1" "NDVI_pva_1" "NDBI_tre_1" "NDBI_pva_1" "Impervio_2"
## [36] "Impervio_3" "Mean_Alb_1" "OBJECTID_1" "geometry"

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spdep)

## Warning: package 'spdep' was built under R version 4.5.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.5.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

stockton_lisa <- stockton %>%
  filter(!is.na(LST_trend1))

neighbors <- poly2nb(
  stockton_lisa,
  queen = TRUE,
  snap = 100
)

weights <- nb2listw(
  neighbors,
  style = "W",
  zero.policy = TRUE
)

moran.test(
  stockton_lisa$LST_trend1,
  weights,
  zero.policy = TRUE
)

## 
##  Moran I test under randomisation
## 
## data:  stockton_lisa$LST_trend1  
## weights: weights    
## 
## Moran I statistic standard deviate = 4.1389, p-value = 1.745e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.261351757      -0.011764706       0.004354402

local_lisa <- localmoran(
  stockton_lisa$LST_trend1,
  weights,
  zero.policy = TRUE
)

stockton_lisa$local_I <- local_lisa[,1]
stockton_lisa$local_p <- local_lisa[,5]
stockton_lisa$local_z <- local_lisa[,4]

#standardize
stockton_lisa$LST_std <- scale(
  stockton_lisa$LST_trend1
)

#spatial lag
lag_LST <- lag.listw(
  weights,
  stockton_lisa$LST_std,
  zero.policy = TRUE
)

stockton_lisa$cluster <- "Not Significant"

# High-High
stockton_lisa$cluster[
  stockton_lisa$LST_std > 0 &
  lag_LST > 0 &
  stockton_lisa$local_p <= 0.05
] <- "High-High"

# Low-Low
stockton_lisa$cluster[
  stockton_lisa$LST_std < 0 &
  lag_LST < 0 &
  stockton_lisa$local_p <= 0.05
] <- "Low-Low"

# High-Low
stockton_lisa$cluster[
  stockton_lisa$LST_std > 0 &
  lag_LST < 0 &
  stockton_lisa$local_p <= 0.05
] <- "High-Low"

# Low-High
stockton_lisa$cluster[
  stockton_lisa$LST_std < 0 &
  lag_LST > 0 &
  stockton_lisa$local_p <= 0.05
] <- "Low-High"

library(tmap)

tm_shape(stockton_lisa) +
  tm_fill(
    "cluster",
    palette = c(
      "red",
      "blue",
      "orange",
      "lightblue",
      "grey"
    ),
    title = "LISA Clusters"
  ) +
  tm_borders()

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

Now for NDBI

# Remove NA values
stockton_ndbi <- stockton %>%
  filter(!is.na(NDBI_tre_1))

# Create neighbors
neighbors_ndbi <- poly2nb(
  stockton_ndbi,
  queen = TRUE
)

# Spatial weights
weights_ndbi <- nb2listw(
  neighbors_ndbi,
  style = "W",
  zero.policy = TRUE
)

# Global Moran's I
moran.test(
  stockton_ndbi$NDBI_tre_1,
  weights_ndbi,
  zero.policy = TRUE
)

## 
##  Moran I test under randomisation
## 
## data:  stockton_ndbi$NDBI_tre_1  
## weights: weights_ndbi    
## 
## Moran I statistic standard deviate = 5.0379, p-value = 2.353e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.318928571      -0.011764706       0.004308748

# Local Moran's I
local_ndbi <- localmoran(
  stockton_ndbi$NDBI_tre_1,
  weights_ndbi,
  zero.policy = TRUE
)

# Attach results
stockton_ndbi$local_I <- local_ndbi[,1]
stockton_ndbi$local_p <- local_ndbi[,5]
stockton_ndbi$local_z <- local_ndbi[,4]

# Standardize NDBI
stockton_ndbi$NDBI_std <- scale(
  stockton_ndbi$NDBI_tre_1
)

# Spatial lag
lag_NDBI <- lag.listw(
  weights_ndbi,
  stockton_ndbi$NDBI_std,
  zero.policy = TRUE
)

# Create cluster categories
stockton_ndbi$cluster <- "Not Significant"

# High-High
stockton_ndbi$cluster[
  stockton_ndbi$NDBI_std > 0 &
  lag_NDBI > 0 &
  stockton_ndbi$local_p <= 0.05
] <- "High-High"

# Low-Low
stockton_ndbi$cluster[
  stockton_ndbi$NDBI_std < 0 &
  lag_NDBI < 0 &
  stockton_ndbi$local_p <= 0.05
] <- "Low-Low"

# High-Low
stockton_ndbi$cluster[
  stockton_ndbi$NDBI_std > 0 &
  lag_NDBI < 0 &
  stockton_ndbi$local_p <= 0.05
] <- "High-Low"

# Low-High
stockton_ndbi$cluster[
  stockton_ndbi$NDBI_std < 0 &
  lag_NDBI > 0 &
  stockton_ndbi$local_p <= 0.05
] <- "Low-High"

# Map
tm_shape(stockton_ndbi) +
  tm_fill(
    "cluster",
    palette = c(
      "red",
      "blue",
      "orange",
      "lightblue",
      "grey"
    ),
    title = "NDBI LISA Clusters"
  ) +
  tm_borders()

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

library(tmap)

cluster_colors <- c(
  "High-High" = "#d73027",
  "High-Low" = "#313695",
  "Low-High" = "#fdae61",
  "Low-Low" = "#abd9e9",
  "Not Significant" = "grey70"
)

# Remove Missing from displayed legend/classes
stockton_lisa$cluster <- factor(
  stockton_lisa$cluster,
  levels = c("High-High", "High-Low", "Low-High", "Low-Low", "Not Significant")
)

stockton_ndbi$cluster <- factor(
  stockton_ndbi$cluster,
  levels = c("High-High", "High-Low", "Low-High", "Low-Low", "Not Significant")
)

map_lst <- tm_shape(stockton_lisa) +
  tm_fill("cluster", palette = cluster_colors, title = "") +
  tm_borders(col = "grey25", lwd = 0.45) +
  tm_layout(
    title = "(A) LST warming trend clusters",
    title.size = 0.9,
    legend.show = FALSE,
    frame = FALSE,
    bg.color = "white",
    inner.margins = c(0.01, 0.01, 0.01, 0.01)
  )

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

map_ndbi <- tm_shape(stockton_ndbi) +
  tm_fill("cluster", palette = cluster_colors, title = "") +
  tm_borders(col = "grey25", lwd = 0.45) +
  tm_layout(
    title = "(B) NDBI trend clusters",
    title.size = 0.9,
    legend.show = FALSE,
    frame = FALSE,
    bg.color = "white",
    inner.margins = c(0.01, 0.01, 0.01, 0.01)
  )

## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

legend_only <- tm_shape(stockton_lisa) +
  tm_fill(
    "cluster",
    palette = cluster_colors,
    title = "LISA Cluster"
  ) +
  tm_layout(
    legend.only = TRUE,
    legend.title.size = 1.0,
    legend.text.size = 0.85,
    bg.color = "white",
    frame = FALSE
  )

## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'

tmap_save(
  map_lst,
  "C:/Users/noahs/Downloads/ES252/Figure3A_LST_LISA.png",
  width = 5.5,
  height = 5.5,
  dpi = 600
)

## Map saved to C:\Users\noahs\Downloads\ES252\Figure3A_LST_LISA.png
## Resolution: 3300 by 3300 pixels
## Size: 5.5 by 5.5 inches (600 dpi)

tmap_save(
  map_ndbi,
  "C:/Users/noahs/Downloads/ES252/Figure3B_NDBI_LISA.png",
  width = 5.5,
  height = 5.5,
  dpi = 600
)

## Map saved to C:\Users\noahs\Downloads\ES252\Figure3B_NDBI_LISA.png
## Resolution: 3300 by 3300 pixels
## Size: 5.5 by 5.5 inches (600 dpi)

tmap_save(
  legend_only,
  "C:/Users/noahs/Downloads/ES252/LISA_Legend.png",
  width = 3,
  height = 1.8,
  dpi = 600
)

## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
## Map saved to C:\Users\noahs\Downloads\ES252\LISA_Legend.png
## 
## Resolution: 1800 by 1080 pixels
## 
## Size: 3 by 1.8 inches (600 dpi)