Spatial Autocorrelation Patterns and Inequality Analysis
Abdi-Basid ADAN
Spatial Analysis of Hotspots and Gini Index for Synthetic Data
Objective Description
This analysis examines spatial autocorrelation patterns in synthetic raster data using the Getis-Ord Gi* statistic to identify hotspots and coldspots, and measures spatial inequality with the Gini index. It involves generating a synthetic raster (representing NDVI), visualizing it as an NDVI map, creating a shapefile, extracting raster values, computing spatial neighbors, calculating Gi* statistics, classifying hotspots, visualizing hotspots with a thematic map, and exporting outputs as shapefiles and rasters. The analysis highlights spatial clustering and inequality in a simulated dataset.
Packages
To install these packages, use:
install.packages('package_name').
library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(raster)## Loading required package: splibrary(terra)## terra 1.7.23library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tmap)
library(ineq)
library(exactextractr)
library(ggplot2)Data Creation
Generate a synthetic raster and a regular grid shapefile, then extract mean raster values for each polygon.
set.seed(42)
# Create raster
Raster <- raster(nrows=100, ncols=100, xmn=0, xmx=10, ymn=0, ymx=10)
values(Raster) <- runif(ncell(Raster), 0, 100)
# Create shapefile (regular grid)
grid <- st_make_grid(st_bbox(c(xmin=0, xmax=10, ymin=0, ymax=10)), cellsize=1, square=TRUE)
shape <- st_sf(geometry=grid)
# Extract mean raster values for each polygon
if (!require(exactextractr)) stop("Package 'exactextractr' is required")
shape$value <- exactextractr::exact_extract(Raster, shape, 'mean')## Warning in .local(x, y, ...): No CRS specified for polygons; assuming they have
## the same CRS as the raster.## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%if (any(is.na(shape$value))) stop("NA values detected in extracted raster data")NDVI Visualization
Visualize the synthetic raster as an NDVI map using a viridis color palette.
tmap_mode("plot")## tmap mode set to plottingtm_shape(Raster) +
tm_raster(palette = "viridis", title = "NDVI") +
tm_legend(outside = TRUE, text.size = 0.8) +
tm_layout(frame = FALSE)
Data Conversion and Preparation
Convert the shapefile to a spatial object and remove polygons with NA values.
shape_sp <- as(shape, "Spatial")
valid_shape <- shape_sp[!is.na(shape_sp$value), ]Spatial Analysis
Calculate spatial neighbors and compute Getis-Ord Gi* statistics.
nb <- poly2nb(valid_shape, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
coords <- coordinates(valid_shape)
localg <- localG(valid_shape$value, lw)
# Create dataframe with results
valid_shape_df <- as.data.frame(valid_shape)
valid_shape_df$localg <- as.numeric(localg)Hotspot Classification
Classify polygons into hotspots and coldspots based on Gi* values.
valid_shape_df <- valid_shape_df %>%
mutate(hotspotsg = case_when(
localg <= -2.58 ~ "Cold spot 99%",
localg > -2.58 & localg <= -1.96 ~ "Cold spot 95%",
localg > -1.96 & localg <= -1.65 ~ "Cold spot 90%",
localg >= 1.65 & localg < 1.96 ~ "Hot spot 90%",
localg >= 1.96 & localg <= 2.58 ~ "Hot spot 95%",
localg >= 2.58 ~ "Hot spot 99%",
TRUE ~ "Not Significant"
)) %>%
mutate(hotspotsg = factor(hotspotsg,
levels = c("Cold spot 99%", "Cold spot 95%",
"Cold spot 90%", "Not Significant",
"Hot spot 90%", "Hot spot 95%",
"Hot spot 99%")))
valid_shape@data <- valid_shape_dfHotspot Visualization
Create a thematic map of hotspots using tmap.
tmap_mode("plot")## tmap mode set to plottingpal <- c("#0000FF", "#4169E1", "#87CEEB", "#D3D3D3", "#FFA07A", "#FF4500", "#FF0000")
names(pal) <- levels(valid_shape_df$hotspotsg)
tm_shape(valid_shape) +
tm_polygons("hotspotsg", palette=pal, title="Hotspots Gi*") +
tm_layout(legend.outside=TRUE, frame=FALSE) +
tm_compass(position=c("right", "top")) +
tm_scale_bar(position=c("left", "bottom"))## Warning: Currect projection of shape valid_shape unknown. Long-lat (WGS84) is
## assumed.
Gini Index Calculation
Convert the raster to polygons and calculate the Gini index to measure spatial inequality.
Raster_polygons <- rasterToPolygons(Raster, dissolve = FALSE)
valid_shape_gini <- Raster_polygons[!is.na(Raster_polygons$layer), ]
print(nrow(valid_shape_gini)) # Number of valid polygons## [1] 10000gini_index <- Gini(valid_shape_gini$layer)
print(paste("Indice de Gini:", round(gini_index, 4)))## [1] "Indice de Gini: 0.3364"# Interpretation
if (gini_index < 0.2) {
cat("Distribution très homogène (inégalité faible)\n")
} else if (gini_index < 0.4) {
cat("Distribution relativement équitable\n")
} else if (gini_index < 0.6) {
cat("Distribution inégale\n")
} else {
cat("Distribution très inégale (forte concentration spatiale)\n")
}## Distribution relativement équitableResults Export
Save the hotspot shapefile and raster.
st_write(st_as_sf(valid_shape), "hotspots_shapefile.shp", delete_layer=TRUE)## writing: substituting ENGCRS["Undefined Cartesian SRS with unknown unit"] for missing CRS## Deleting layer `hotspots_shapefile' using driver `ESRI Shapefile'
## Writing layer `hotspots_shapefile' to data source
## `hotspots_shapefile.shp' using driver `ESRI Shapefile'
## Writing 100 features with 3 fields and geometry type Polygon.valid_shape$hotspot_num <- as.numeric(valid_shape$hotspotsg)
raster_hotspots <- rasterize(valid_shape, Raster, field="hotspot_num")
writeRaster(raster_hotspots, "hotspots_raster.tif", format="GTiff", overwrite=TRUE)```