Mapping LISA cluster in R

General
, , , ,
1

Hi,
I am trying to map the bivariate spatial correlation in R (bivariate LISA). When I tried to run -- nb <- poly2nb(map)-- It is showing -- Error in poly2nb(map) : Not a polygon object.

Maybe I am doing something wrong. Please, I appreciate it if you could assist me. My code is here--

rm(list=ls(all=TRUE))
library(stringr)
library(spdep)
library(rgdal)
library(magrittr)
library(ggplot2)
library(sf)
library(sp)
library(rgdal)
library(lwgeom)
library(dplyr)
library(maptools)
library(plyr)
#======================================================
# Load data
setwd("D:/P_T/Europe")
file <-read.csv("Europe_points.csv", header = T)
Europe <-data.frame(file)
head(Europe)

#### Convert point data into spatial object
#map<-st_as_sf(Europe, coords = c("Longitude", "Latitude"), crs = 4326) 
coords <- cbind(Europe$Longitude, Europe$Latitude)
eu_data <- SpatialPointsDataFrame(coords = coords, data = Europe, proj = CRS("+proj=longlat +datum=WGS84"))

## read shapefile
countries <- readOGR("D:/shape_files/ne_110m_admin_0_countries", layer="ne_110m_admin_0_countries") 
sub_Europe <- countries[countries$CONTINENT == "Europe",]
#EU_proj <- spTransform(sub_Europe, CRS("+proj=longlat +datum=WGS84"))
sub_Europe@data$id <- rownames(sub_Europe@data)

# concer shp into dataframe
df_sub_Europe <-fortify(sub_Europe)
df_sub_Europe <- join(df_sub_Europe, sub_Europe@data, by="id")

## Now merge points attributes with shapefile
map <- merge(df_sub_Europe, eu_data)


# Variables to use in the correlation
x <- map[map$Variable=="Temp",]
y <- map[map$Variable=="Precip",]

#======================================================
# Some functions

# Bivariate Moran's I
moran_I <- function(x, y = NULL, W){
  if(is.null(y)) y = x
  
  xp <- scale(x)[, 1]
  yp <- scale(y)[, 1]
  W[which(is.na(W))] <- 0
  n <- nrow(W)
  
  global <- (xp%*%W%*%yp)/(n - 1)
  local  <- (xp*W%*%yp)
  
  list(global = global, local  = as.numeric(local))
}


# Permutations for the Bivariate Moran's I
simula_moran <- function(x, y = NULL, W, nsims = 2000){
  
  if(is.null(y)) y = x
  
  n   = nrow(W)
  IDs = 1:n
  
  xp <- scale(x)[, 1]
  W[which(is.na(W))] <- 0
  
  global_sims = NULL
  local_sims  = matrix(NA, nrow = n, ncol=nsims)
  
  ID_sample = sample(IDs, size = n*nsims, replace = T)
  
  y_s = y[ID_sample]
  y_s = matrix(y_s, nrow = n, ncol = nsims)
  y_s <- (y_s - apply(y_s, 1, mean))/apply(y_s, 1, sd)
  
  global_sims  <- as.numeric( (xp%*%W%*%y_s)/(n - 1) )
  local_sims  <- (xp*W%*%y_s)
  
  list(global_sims = global_sims,
       local_sims  = local_sims)
}



#======================================================
# Adjacency Matrix (Queen)
nb <- poly2nb(map) 
#or #nb = st_overlaps(map,sparse=FALSE)
#nb<-st_intersects(map)
lw <-nb2listw(nb,style="B", zero.policy = TRUE)
W  <- as(lw, "symmetricMatrix")
W  <- as.matrix(W/rowSums(W))
W[which(is.na(W))] <- 0


#======================================================
# Calculating the index and its simulated distribution
# for global and local values

m <- moran_I(x, y, W)

# Global Moral
global_moran <- m[[1]][1]
#> 0.2218409

# Local values
m_i <- m[[2]] 

# local simulations
local_sims <- simula_moran(x, y, W)$local_sims


# global pseudo p-value  
# get all simulated global moran
global_sims <- simula_moran(x, y, W)$global_sims

# Proportion of simulated global values taht are higher (in absolute terms) than the actual index 
moran_pvalue <- sum(abs(global_sims) > abs( global_moran )) / length(global_sims)
#> 0


# Identifying the significant values 
alpha <- .05  # for a 95% confidence interval
probs <- c(alpha/2, 1-alpha/2)
intervals <- t( apply(local_sims, 1, function(x) quantile(x, probs=probs)))
sig       <- ( m_i < intervals[,1] )  | ( m_i > intervals[,2] )

#======================================================
# Preparing for plotting
# Convert shape file into sf object
map_sf     <- st_as_sf(map)
map_sf$sig <- sig


# Identifying the LISA clusters
xp <- scale(x)[,1]
yp <- scale(y)[,1]


patterns <- as.character( interaction(xp > 0, W%*%yp > 0) )
patterns <- patterns %>% 
  str_replace_all("TRUE","High") %>% 
  str_replace_all("FALSE","Low")

patterns[map_sf$sig==0] <- "Not significant"
map_sf$patterns <- patterns


# Rename LISA clusters
map_sf$patterns2 <- factor(map_sf$patterns, levels=c("High.High", "High.Low", "Low.High", "Low.Low", "Not significant"),
                           labels=c("High Temp - High Precip", "High Temp - Low Precip",
                                    "Low Temp - High Precip", "Low Temp - Low Precip",
                                    "Not significant"))

### PLOT

ggplot() +
  geom_sf(data=map_sf, aes(fill=patterns2), color="NA") +
  scale_fill_manual(values = c("red", "pink", "light blue", "dark blue", "grey80")) + 
  guides(fill = guide_legend(title="LISA clusters")) +
  theme_minimal()

My data looks like this---

Century Variable Mean Region Latitude Longitude Continent
1190 Temp 0.35375689 Kongressvatnet 78.0217 13.9311 Europe
1090 Temp -0.345180401 Lake Redon 42.6425 0.770278 Europe
1190 Temp 0.097248967 Lake Redon 42.6425 0.770278 Europe
990 Precip -0.26756433 Island of Grimsey 66.5265 -18.1957 Europe
1370 Precip 0.326454087 Island of Grimsey 66.5265 -18.1957 Europe
990 Precip 1.096977506 Island of Grimsey 66.5265 -18.1957 Europe
1370 Precip -1.024937192 Island of Grimsey 66.5265 -18.1957 Europe
1090 Temp 1.085211722 Seebergsee 46.116667 7.466667 Europe
1190 Temp 1.347592222 Seebergsee 46.116667 7.466667 Europe
1370 Temp -2.217270397 Seebergsee 46.116667 7.466667 Europe
1090 Precip -1.400811763 Uamh an Tartair 58.15 -4.983333 Europe
1190 Precip -0.092810391 Uamh an Tartair 58.15 -4.983333 Europe
1370 Precip 0.638131552 Uamh an Tartair 58.15 -4.983333 Europe
810 Temp -0.825153678 Spannagel Cave 47.09 11.67 Europe
860 Temp 0.5622551 Spannagel Cave 47.09 11.67 Europe
990 Temp 0.684673522 Spannagel Cave 47.09 11.67 Europe
1090 Temp 0.154193695 Spannagel Cave 47.09 11.67 Europe
1190 Temp -0.498704554 Spannagel Cave 47.09 11.67 Europe
1370 Temp -0.743541397 Spannagel Cave 47.09 11.67 Europe
810 Temp -0.230349118 Northern inland Iberia 42.69405 -3.9446 Europe
860 Temp 0.044131306 Northern inland Iberia 42.69405 -3.9446 Europe
990 Temp 0.335766755 Northern inland Iberia 42.69405 -3.9446 Europe
1090 Temp 0.136768448 Northern inland Iberia 42.69405 -3.9446 Europe
2

Your code is somewhat complex, and gives the impression of being a combination of several pre-existing scripts. As consequence it is not easy to comment upon.

But it seems your map object comes, indirectly via the df_sub_Europe object, from the result of a ggplot2::fortify() call.

This is an oldish function, used for plotting (no wonder, as it lives in ggplot...) and currently mostly superseded by ggplot2::geom_sf(). It outputs a regular data.frame. This kind of object is not an acceptable input for a spdep::poly2nb() call.

You need to provide the queen weights call with an sf or sfc object from {sf} context, or SpatialPolygons from the older {sp} context.

3

Hi, thank you so much for the info. I have points data (CSV format), and I want to estimate LISA (with two variables== x-Temp and y-Precip). I tried to convert point data into spatial object, like this==data_sf <- st_as_sf(file, coords = c('Latitude', 'Longitude'), crs="+proj=longlat +datum=WGS84"). Still it is showing = Error in poly2nb(data_sf) : Polygon geometries required. I am new to R, and difficult to create neighborhoods that the poly2nb function could easily track. Or maybe I am doing something wrong. Plz, I would appreciate it if you could guide me :slight_smile:

4

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.