Working with R and QGIS

Question

I am performing an interpolation in R but i get the following results. What could be the issue? The script

##MySript=group
##rainfall= vector
##output= output raster
    require(sp)
    require(rgdal)
    rainfall
    slotNames(rainfall)
    rainfall@proj4string
    rainfall@data
    att.table = rainfall@data
    class(att.table)
    dim(att.table)
    dimnames(att.table)
    summary(rainfall)
    plot(rainfall)
    require(geoR)
    require(sp)
    require(maptools)
    require(raster)
    require(geoRglm)
    require(rgeos)
    require (rgdal)
    locs=pred_grid(c(35.06,35.72), c(-1.63,-1.0), by = 0.002)
    data1=as.geodata(rainfall)
    plot(variog(data1, uvec=seq(0,1, by = 0.05)))
    LE1=likfit(data1,cov.model="matern",ini=c(0.0002,0.06), 
         nugget =0.000002,kappa=1.5)
    KC = krige.control(type = "sk", obj.mod = LE1)
    sk1 = krige.conv(data1, krige = KC, loc = locs)
    map=image(sk1)

The result

require(geoR)
Loading required package: geoR
--------------------------------------------------------------
Analysis of Geostatistical Data
For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
geoR version 1.7-5.2 (built on 2016-05-02) is now loaded
--------------------------------------------------------------

    require(sp)
    require(maptools)
    Loading required package: maptools
    Checking rgeos availability: TRUE
    require(raster)
    require(geoRglm)
    Loading required package: geoRglm
---------------------------------------------------------

A Package for Generalised Linear Spatial Models


geoRglm version 0.9-11 (2017-10-17) is now loaded

-----------------------------------------------------------



require(rgeos)
Loading required package: rgeos
rgeos version: 0.3-26, (SVN revision 560)
GEOS runtime version: 3.6.1-CAPI-1.10.1 r0
Linking to sp version: 1.2-7
Polygon checking: TRUE

    require (rgdal)
    locs=pred_grid(c(35.06,35.72), c(-1.63,-1.0), by = 0.002)
    data1=as.geodata(rainfall)
    plot(variog(data1, uvec=seq(0,1, by = 0.05)))
    variog: computing omnidirectional variogram
    LE1=likfit(data1,cov.model="matern",ini=c(0.0002,0.06), 
           nugget =0.000002,kappa=1.5)
---------------------------------------------------------------
likfit: likelihood maximisation using the function optim.
likfit: Use control() to pass additional
arguments for the maximisation function.
For further details see documentation for optim.
likfit: It is highly advisable to run this function several
times with different initial values for the parameters.
likfit: WARNING: This step can be time demanding!
---------------------------------------------------------------
likfit: end of numerical maximisation.
KC = krige.control(type = "sk", obj.mod = LE1)
sk1 = krige.conv(data1, krige = KC, loc = locs)
krige.conv: model with constant mean
krige.conv: Kriging performed using global neighbourhood
map=image(sk1)


writeRaster(output,"C:/PHDproposal/R_spatial/R_Spatial Exx/data/map.tif", overwrite=TRUE)
Error in writeRaster(output, "C:/PHDproposal/R_spatial/R_Spatial Exx/data/map.tif", :
object 'output' not found
Execution halted

looks like you haven't defined output. But you really should take the time to properly format your code if you're trying to get good help. — bwp8nt, Aug 28 '18 at 17:02
@bwp8nt thank you for the comment. i have defined but i still get the same response — Eric Wamiti, Aug 29 '18 at 06:56

Spacedman · Accepted Answer · 2018-08-29T07:36:22.780

0

When you have

##output= output raster

in an R processing script in QGIS that means that at the end of the script QGIS will use R to save a raster from an object called output that you have generated in your script. So for example, if you have:

m = matrix(runif(16),4,4)
output = raster(m)

in your script then QGIS will create a 4x4 raster at coordinates (0,0) to (1,1), and that will become a layer in QGIS. This is what the writeRaster line in the error is doing. It fails because you didn't create a raster object called output in your script.

What you did do was map = image(sk1) which draws an image, which isn't what a processing script should do if you want it to output a new raster layer. You possibly want to convert the sk1 object to a raster package "raster" object called output.

You can create a raster from krige.conv output this way - I'm using the example from help(krige.conv):

 > library(raster)
 > output = rasterFromXYZ(cbind(loci, kc$predict))

replace loci with your prediction coordinates (as a 2-column matrix) and kc with your output from krige.conv. If you want the kriging variance then use kc$krige.var.

I think if you want your processing script to return two raster layers to QGIS, then have two ## directives naming two R objects, and then do something like:

  ##pred= output raster
  ##pvar= output raster

   # insert code here...

  pred = rasterToXYZ(cbind(loci, kc$predict))
  pvar = rasterToXYZ(cbind(loci, kc$var))

edited Aug 29 '18 at 07:36

answered Aug 29 '18 at 06:56

Spacedman

63,755
5
81
115

Thank you. sk1 is the prediction data and locs is the prediction grid. Possibly if i can call the matrix that joins the data to the grid, i will achieve the intended result which should be an image. – Eric Wamiti Aug 29 '18 at 07:19
in R image(sk1) gives and image of values in the prediction grid "locs" which is what i would wish to get in QGIS so that i can use this information with other images – Eric Wamiti Aug 29 '18 at 07:22
@ Spacedman When I try the Bayesian method in R for QGIS using “output = rasterFromXYZ(cbind(locss, skb$predict))” in the script below, I get a miss match error. I also realize “locss” has 1056 elements and “sbk ” has 2000 elements. How do I use rasterFromXYZ in this situation? – Eric Wamiti Aug 30 '18 at 07:17
locss=pred_grid(c(35.06,35.7), c(-1.63,-1.0), by = 0.02) MC <- model.control(trend.d = "1st", trend.l = "1st",cov.model="matern",kappa = 1.5) PC <- prior.control(phi.discrete = seq(0, 0.05, l = 40), phi.prior = "reciprocal", tausq.rel.prior = "unif", tausq.rel.discrete = seq(0, 0.005, l = 20)) OC <- output.control(n.post = 2000, moments = T) skb <- krige.bayes(data1, loc = locss,model = MC,prior = PC,output = OC) – Eric Wamiti Aug 30 '18 at 07:18
That sounds like a new question, and hard to show in comments here. Make a new question with a clear R example, take out the references to QGIS and just ask how to create a raster from the output of krige.bayes. – Spacedman Aug 30 '18 at 07:24
i have reposted – Eric Wamiti Aug 30 '18 at 08:35

Working with R and QGIS

1 Answers1