[R-sig-Geo] Help with universal kriging using gstat

[Roger Bivand]

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 using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en