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

[Thomas Adams]

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"’

> 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

-- 
Thomas E Adams, III
1724 Sage Lane
Blacksburg, VA 24060
tea3rd using gmail.com (personal)
tea using terrapredictions.org (work)

1 (513) 739-9512 (cell)

	[[alternative HTML version deleted]]