I have two rasters. The first raster is a region showing sinkholes(depressions) of the region(cell values are 1 for sinkholes and 0 for non-sinkhole cells). The second raster is the TPI(Topographic Positioning System) of the same region. I want to conduct a logistic regression analysis between the two rasters(Sinkhole~TPI).

How can I do this using R?

(The purpose is to build a sinkhole detection model using TPI index so that it can use the model to detect sinkholes in other parts of the study area. )

enter image description here enter image description here


3 Answers 3


To eleborate on Spacedman's answer (and to make it more clear where the problem may reside:

Can you combine the data?

s <- stack(sinks, tpi)

If so, you should be able to do:

v <- data.frame(na.omit(values(s)))
m <- glm(sinks~tpi, data=v, family=binomial)


p <- predict(tpi, m, type="response")

If your rasters are very large, you may use a sample instead. I think it is better to use a regular (not a random) sample to minimize redundancy in your data.

v <- data.frame(sampleRegular(s, 10000))  
m <- glm(sinks~tpi, data=v, family=binomial)

I would not be much concerned about spatial autocorrelation in residuals, as your purpose is prediction, not inference.

Now with some example data:

elevation <- getData('alt', country='CHE')
tpi <- terrain(elevation, opt='tpi', unit='degrees')
sinks <- elevation < 500
names(sinks) <- "sinks"
s <- stack(sinks, tpi)

v <- data.frame(na.omit(values(s)))
m <- glm(sinks~tpi, data=v, family=binomial)

p <- predict(tpi, m, type='response')

It is a terrible model fit and prediction, but it shows that the mechanism works.

Here is how you can predict to another raster:

elevation2 <- getData('alt', country='LUX')
tpi2 <- terrain(elevation2, opt='tpi', unit='degrees')
p <- predict(tpi2, m, type='response')

If your rasters have the same grid size and position, and the sinks raster is 0s and 1s (or TRUEs and FALSEs) then its just like doing a logistic regression any other way in R, but getting the values out of the raster using []:

m = glm(sinks[]~tpi[],family=binomial)

BUT this model assumes independent errors, which is patently not true for most spatial data. What sort of spatial model would you like to fit? Once you've figured that out, ask another question.

  • Thanks for your reply. I tried it but I get the error: Error in model.frame.default(formula = boun[] ~ tpi[], drop.unused.levels = TRUE) : variable lengths differ (found for 'tpi[]') both rasters have the same spatial reference and the same size, and also both have the same cell size.
    – Arsalan
    Jan 25, 2018 at 22:08
  • Without further evidence I can't believe that. When I fit a model with two objects of different size and I get that exact error. What is ncell(tpi) and ncell(boun) - that will tell you.
    – Spacedman
    Jan 25, 2018 at 22:30

To supplement the answers you already got, I believe that one viable strategy would entail:

  • drawing a number of random sampling points within your study area;
  • drawing them with a subjective minimum distance in order to possibly alleviate spatial autocorrelation of your model residuals (more on this later on);
  • extract the values of your Dependent Variable and of your Predictor to your sampling points;
  • perform Logistic Regression in R (you already got some advices on how to go about this);
  • checking if the model is significant;
  • using the estimated constant and predictor's coefficient (Beta) to plot a fitted model raster;
  • getting the model residuals;
  • plotting a correlogram to see if spatial autocorrelation exists among the model's residuals.

As for random sampling points, you could use the sp R package which has a facility to draw random points within a given study plot. Using a minimum distance may potentially reduce the autocorrelation among the residuals, which you will nontheless formally check that afterward by using the correlogram.

I happen to describe a similar procedure, even though from a ArcGIS perspective, in an earlier answer of mine in this same forum (also providing a couple of references to literature): How to explain the spatial relation of forest loss hotspots with several others variables

Not the answer you're looking for? Browse other questions tagged or ask your own question.