Hi Roger, Thank you... I finally figured out what I was doing wrong and have been successful with gstat universal kriging. I will write-up the details for others later today and plan to write a technical note of such for GRASS GIS users, which I will submit there.
Best, Tom On Thu, Jul 16, 2020 at 4:15 AM Roger Bivand <roger.biv...@nhh.no> wrote: > 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.biv...@nhh.no > https://orcid.org/0000-0003-2392-6140 > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en -- Thomas E Adams, III 1724 Sage Lane Blacksburg, VA 24060 tea...@gmail.com (personal) t...@terrapredictions.org (work) 1 (513) 739-9512 (cell) [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo