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). can anybody kindly tell me how I can do this?(preferably 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

  • Could you please edit your question to include what you hope to achieve using regression? What does your desired end product/analysis look like? – Aaron Jan 25 '18 at 3:44
  • 1
    @Aaron I edited the question. This region is only a small part of the study area, which I want to build my model based on, and use it for the other parts of the study area. – Arsalan Jan 25 '18 at 4:05

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')
  • Thanks for your reply. I made the model. Then I tried to use the model on another tpi_raster. I get this error: R Error in p[-naind, ] <- predv : number of items to replace is not a multiple of replacement length In addition: Warning message: 'newdata' had 2493378 rows but variables found have 14983552 rows What is the reason? I'm new to R and cant understand this freaking error. – Arsalan Jan 29 '18 at 19:08
  • 1
    I have added some code. Does that help? – Robert Hijmans Jan 29 '18 at 19:53

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 '18 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 '18 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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