Help with universal kriging using gstat

On Wed, 15 Jul 2020, Thomas Adams wrote:

Hi all,

It's been some time since I approached universal kriging using gstat (I
struggled with this previously, years ago:
https://stat.ethz.ch/pipermail/r-sig-geo/2006-May/001017.html).

The problem...

Within GRASS GIS, using R, I do this...

(1) read a raster DEM into R from GRASS -- srtm <-
readRAST("mozambique_srtm_patch",cat=FALSE)
(2) read GRASS point data consisting of 4 fields (category, lon, lat,
temperature) -- airtemp <- readVECT("Mozambique_air_temp_2017_ann")
(3) so I have a SpatialGridDataFrame and SpatialPointsDataFrame,
respectively
(4) I can do the following: plot srtm, airtemp, generate an interpolated
grid of air temperatures with krige, using the airtemp
SpatialPointsDataFrame, and overlay the various data for visualization

What I want to do is to use the srtm DEM data as a secondary trend

variable

to spatially interpolate airtemp using universal kriging. I cannot figure
out how to construct the data sets and use krige in gstat to do this. I
have spent several days scouring the internet for an example (including
previous queries of my own, cited above) to no avail.

It seems I should be able to do this, essentially:

dem<-read.asciigrid("gtopo30.dem")
class(dem)

[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"

image(dem)
points(Y ~ X, data=temps)
class(temps)

[1] "data.frame"

coordinates(temps)=~X+Y
dem.ov=overlay(dem,temps)
summary(dem.ov)

vgm <- vgm(psill=8,model="Exp",range=600000,nugget=3.8)
vgm_temps_r<-fit.variogram(variogram(T~gtopo30.dem,temps), model=vgm)
plot(variogram(T~gtopo30.dem,temps),main = "fitted by gstat")
temps_uk<-krige(T~gtopo30.dem,temps,dem, vgm_temps_r)

[using universal kriging]

library(lattice)
trellis.par.set(sp.theme())
spplot(temps_uk,"var1.pred", main="Universal kriging predictions

However, when I take this step:

dem.ov=overlay(srtm,airtemp)

Error in (function (classes, fdef, mtable)  :
 unable to find an inherited method for function ?overlay? for signature
?"SpatialGridDataFrame", "SpatialPointsDataFrame"?

The overlay() method was retired long ago in favour of over(), and the
order of the arguments was standardised. So over(airtemp, strm) should
return the values of strm at the airtemp measuement points.

Roger

airtemp$srtm.dem=dem.ov$srtm  <====== this fails (see below)

vgm <- vgm(psill=8,model="Exp",range=600000,nugget=3.8)
vgm_temps_r<-fit.variogram(variogram(temp~srtm.dem,airtemp), model=vgm)
temps_uk<-krige(temp~srtm.dem,airtemp,srtm, vgm_temps_r)

airtemp$srtm.dem=srtm

Error in `[[<-.data.frame`(`*tmp*`, name, value =
new("SpatialGridDataFrame",  :
 replacement has 72821592 rows, data has 267

Any suggestions?

Best,
Tom

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en