[GIS] Spatial interpolation from categorical data in R

interpolationkernel densityrraster

I have a spatial point data set, where each point represents a person pronouncing a word in a specific fashion. Let's say I have six different pronunciations.

In an earlier iteration, I aggregated the point data to a hexagonal grid, the color represents the dominant pronunciation in the hexagon while the opacity represents its dominance, see: https://www.srf.ch/static/srf-data/data/2016/dialektkarten-karte/

Now I'd like to spatially interpolate the data to get a continuous surface, similar as Josh Katz did: http://www4.ncsu.edu/~jakatz2/files/dialectposter.png

He uses something he calls "k-nearest neighbor kernel smoothing" but he's not publishing his R code so I don't fully understand how that is supposed to work with categorical data. I perfectly understand on how to do it with continuous variables, but not with the latter.

The pseudo-algorithm I came up with so far:

Repeat the "final estimate" formula for each pronunciation, using yi = 1 for points where that pronunciation is used, yi = 0 if not. By doing so, you get a grid for each pronunciation. The cells represent the probabilities of the respective pronunciation appearing at that grid cell.
Now, find the most probable pronunciation in each grid cell. This determines the cell color (hue). Its probability determines the color value (i.e.: the higher, the darker, and vice versa).
Find a way to map that resulting grid / raster – might be tricky because it contains two variables that need to be mapped to visual variables (dominant pronunciation and its probability).

But maybe this is completely off – and also, I wonder whether there are R packages for spatial interpolation of categorical data. So far I have only found knn methods and the like and I suppose I have to build something myself.

Best Answer

In the end, I figured out a method with the kknn package. If somebody is still interested in this answer, I recommend reading my blog post on the topic.

Related Solutions

[GIS] How to perform loop on raster files

I am not sure I completely understand the question, but if it is to sample points from each state raster according to the probability found within the raster, I think this solution should do:

# Bring in libraries and set random seed
library(data.table)
library(raster)
set.seed(123)

# Bring in tifs and subset to your 16 states
tifFiles <- list.files(".", pattern=".TIF$", full.names = TRUE)
tifFiles <- paste0("./Suitability_Prob_v",c(30:39,310:315),".TIF")

# How many points do we want to sample from each state?
# I am not sure what is in your csv file, so I sampled some points from a poisson
sampPoints <- data.table(State=tifFiles, N=rpois(length(tifFiles),6))

# Function to read in raster and sample points
readRasterAndSample <- function(rasFile, count) {

    # Read in raster and deal with 0s and NAs
    ras <- raster(rasFile)
    ras[ras < 0] <- 0

    # Convert to data.table
    rasDT <- as.data.frame(ras, xy=TRUE)
    setDT(rasDT)
    setnames(rasDT, c("x","y","values"))

    # Sample non-zero points according to desired count
    # Use raster probability as sample prob
    outPoints <- rasDT[values > 0]
    outPoints <- outPoints[sample(1:.N, count, prob=values, replace=TRUE)]  

    # Jitter points randomly
    outPoints[, x := x + runif(.N, -0.5*res(ras), 0.5*res(ras))] 
    outPoints[, y := y + runif(.N, -0.5*res(ras), 0.5*res(ras))] 
    outPoints[, state := basename(rasFile)]

    # Combine into list and return
    outList <- list(ras=ras, pts=outPoints)
    return(outList)
}

# Apply function to all states
rasterList <- lapply(1:nrow(sampPoints), function(x) 
    with(sampPoints[x], readRasterAndSample(rasFile=State, count=N)))

# Collect all rasters and merge
rasterAll <- do.call(merge, lapply(rasterList, "[[", "ras"))

# Collect all state points and convert to SpatialPointsDataframe
pointDat <- rbindlist(lapply(rasterList, "[[", "pts"), fill=TRUE)
coordinates(pointDat) <- ~x+y
proj4string(pointDat) <- rasterAll@crs@projargs

# Plot to check
plot(rasterAll)
plot(pointDat, add=TRUE)

[GIS] Accounting for spatial autocorrelation in rasters

The autocorrelation of covariates is not the problem in itself. What may be a problem is if there is correlation of covariates at the position of data observations. In that case, there will be identifiability issues to estimate the parameters of your SDM model.
What people usually do is to test for correlation between covariates at observation points. When two covariates are correlated, you can combine with PCA like you plan to do or choose one of the two covariates (with biological a priori sense).
What I do is to test a model with the first covariate, then the second one and then both together. But the test is a k-fold cross-validation so that useless additional covariates will not be retained, hence if the correlation between covariates is too high, only the best one is retained. But if there is a little information coming from both of them, you may want to keep both of them. Even if there is correlation. Environmental covariates always shows correlation because temperature is linked to altitude, because rain is linked to wind, etc...
Even if you have high resolution covariates, I would recommend to disaggregate them and test combination of all covariates, both with high resolution and lower resolution. Resolution will capture different scale effects and you may be surprised by the outputs.
By the way, the coincidence (but no correlation...) is that I just released a R-package on github that may help you to do that. I present it on my website and the "SDM_Selection" vignette will show you my own way of doing SDM and covariates selection: https://statnmap.com/sdmselect-package-species-distribution-modelling/

Edit

The model will be built on your point dataset. If there are identifiability or multicollinearity issues, it will be because of your dataset, not the external data you do not use for modeling.

You can see it the other way: There may be some environmental data completely not related/correlated, but your sampling plan is biased. e.g. imagine each time you observed in the forest it was a rainy month and each time you went for observations in the swamp it was a sunny month, hence there will be correlation between cover type and rain monthly rate. At the same latitude and altitude, let's say there is few possibility of general correlation between these two covariates, but in your dataset, these two will be correlated. Similarly, you may find correlation between global covariates but not in your dataset because your sampling plan counterbalanced the correlation. Because the correlation is not 100%, this means there are chances for a combination of sampling position for which there is no correlation in the covariates.
Thus, you need to verify this correlation inside your dataset.

There is a recent blog post about multicollinearity: https://datascienceplus.com/multicollinearity-in-r/

Best Answer

Related Solutions

[GIS] How to perform loop on raster files

[GIS] Accounting for spatial autocorrelation in rasters

Edit

Related Question