-4

I have conducted a logistic regression in R using two rasters. I want to export it as a raster layer but I keep getting errors. the code is as follows:

S <- stack(sinks, tpi, beer)`
rdata <- data.frame(na.omit(values(s)))
model <- glm(sinks[]~tpi[]+beer[], data=rdata, family=binomial)
p = predict(model, newdata=data.frame(tpi=tpi[],beer=beer[]))
rf <- writeRaster(p, filename="test.tif", format="GTiff")

Error in (function (classes, fdef, mtable)  : 
unable to find an inherited method for function ‘writeRaster’ for signature
‘"numeric", "character"’

Can anybody solve this problem?

closed as unclear what you're asking by Spacedman, aldo_tapia, BERA, whyzar, xunilk Jan 30 '18 at 19:10

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • Please make a reproducible example, perhaps with random data. I've answered some of your questions before and I don't want to go through setting up rasters again, so I'll give a quick answer in comments but not a full one unless you can make it easy for us! – Spacedman Jan 30 '18 at 17:34
  • 2
    p is not a raster - its a vector of values. You need to put those values back into a raster before you can write it but that might be complicated by your omission of NA values. Also your stack is S but your rdata comes from s and you are modelling with sinks[] which might not be the data in rdata so I have no idea what you are doing and this is confusing... – Spacedman Jan 30 '18 at 17:36
4

First, make some sample data. You should be doing this so all we have to do is cut and paste to see what your problem is:

library(raster)
tpi = raster(matrix(runif(50),5,10))
beer = raster(matrix(runif(50),5,10))
sinks = raster(matrix(runif(50)>.5,5,10))

I'll set some values to NA because I suspect you have some:

tpi[2,3]=NA
beer[1,1]=NA

Then you do this:

S <- stack(sinks, tpi, beer)
names(S)
# [1] "layer.1" "layer.2" "layer.3"

Note that the layer names are not the names of the individual rasters.

Then you so (with S replacing s):

 rdata <- data.frame(na.omit(values(S)))

Always check what you get, in this case:

 head(rdata)
 #  layer.1   layer.2    layer.3
 # 1       0 0.2671858 0.80993009
 # 2       1 0.5794600 0.95555964
 # 3       1 0.9979137 0.31638943
 dim(rdata)
 # [1] 48  3

a data frame with those layer names, and 48 rows because we've dropped the two NAs thanks to na.omit.

You then try and model with:

model <- glm(sinks[]~tpi[]+beer[], data=rdata, family=binomial)

which is going to take the values from the rasters you started with, not the rdata data frame you constructed. What you really want at this point is probably:

model <- glm(layer.1~layer.2+layer.3, data=rdata, family=binomial)

referring to the column names in rdata for fitting.

Then you try predictions:

p = predict(model, newdata=data.frame(tpi=tpi[],beer=beer[]))

which should fail because your model is in terms of layer.2 and layer.3 as covariates. So make a data frame with those names:

p = predict(model, newdata=data.frame(layer.2=tpi[],layer.3=beer[]))

Check for surprises:

summary(p)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
# -0.69110  0.06903  0.43210  0.45960  0.81990  1.50900        2 

That's a vector of length 50 with 2 NAs in it. Its not a raster.

However it has the values ready to put into a raster. We first make an empty raster of the same shape as one of your input rasters, and then put the values in:

praster = raster(beer)
praster[] = p

now you can:

writeRaster(praster, "praster.tif")

Appendix

For readability you can rename the stack layers and that follows through into the data frame:

> names(S)=c("sinks","tpi","beer")
> data.frame(na.omit(values(S)))
   sinks        tpi       beer
1      0 0.26718584 0.80993009
2      1 0.57945999 0.95555964
3      1 0.99791372 0.31638943
4      0 0.73093333 0.95536213

but you have to understand when you are getting values from the data frame and when R is getting values from the raster objects with those names.

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