Based on the sample data and code below, I am getting an error while running GWmodel::gwss function to compute the spatial Pearson autocorrelation values and plot it on a map.
Purpose:
To see if there is correlation between AvgTMin and AvgPpt and plot correlation values spatially as a map.
Error from the gwss function:
Error in gw_weight_mat(vdist, bw, kernel, adaptive) :
Expecting a single value: [extent=0].
In addition: Warning message:
In proj4string(data) :
CRS object has comment, which is lost in output; in tests, see
https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
What could be causing the issue and how can I fix this?
Please note that the dataset is gridded, what that means is that both AvgTMin and AvgPpt have the same coordinates. So, if I am using the wrong function to compute the spatial autocorrelation values between these two variables, kindly suggest an appropriate approach.
Desired Output (screenshot provided from another similar analysis with different variables and non-grid point data):
Sample Data (dummy_sp):
new("SpatialMultiPointsDataFrame", data = structure(list(AvgTMin = c(66.8129691710839,
65.7261495113373, 69.0089753627777, 67.0250748157501, 66.838925075531
), AvgPpt = c(186.947467290438, 111.889091470025, 114.544090270996,
154.866916656494, 122.16354282697), year = 1980:1984, CITYNAME = c("a",
"a", "a", "a", "a")), class = "data.frame", row.names = c(NA,
5L)), coords = list(structure(c(-80.2499999994118, -80.1666666660779,
26.3333333333333, 26.25), dim = c(2L, 2L)), structure(c(-80.2499999994118,
-80.1666666660779, 26.3333333333333, 26.25), dim = c(2L, 2L)),
structure(c(-80.2499999994118, -80.1666666660779, 26.3333333333333,
26.25), dim = c(2L, 2L)), structure(c(-80.2499999994118,
-80.1666666660779, 26.3333333333333, 26.25), dim = c(2L,
2L)), structure(c(-80.2499999994118, -80.1666666660779, 26.3333333333333,
26.25), dim = c(2L, 2L))), bbox = structure(c(-80.2499999994118,
26.25, -80.1666666660779, 26.3333333333333), dim = c(2L, 2L), dimnames = list(
NULL, c("min", "max"))), proj4string = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs"))
Code:
library(GWmodel)
##=============================================================
#quick.map function
##=============================================================
quick.map_function = function(spdf,var,legend.title,main.title)
{
x = spdf@data[,var]
cut.vals = pretty(x)
x.cut = cut(x,cut.vals)
cut.levels = levels(x.cut)
8
cut.band = match(x.cut,cut.levels)
colors = brewer.pal(length(cut.levels), "OrRd")
par(mar=c(1,1,1,1))
plot(Columbus,col="gray10",border="gray50", bg="gray50")
title(main.title)
plot(spdf,add=TRUE,col=colors[cut.band],pch=16, cex=0.7)
legend("topleft",cut.levels,col=colors,pch=16,bty="n",title=legend.title, cex=1)
}
# Using the gwss function
gwss.col = gwss(dummy_sp,
vars = c("AvgTMin", "AvgPpt"),
adaptive = T,
bw = 1,
kernel = "Gaussian")
# Plot the AvgTMin and AvgPpt coefficients results
quick.map_exam2(gwss.col$SDF,"Corr_AcgTMin.AvgPpt","Coefficients","Geographically Weighted Pearson Correlation")
