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

[Thomas Adams]

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.Bivand using 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.Bivand using 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
tea3rd using gmail.com (personal)
tea using terrapredictions.org (work)

1 (513) 739-9512 (cell)

	[[alternative HTML version deleted]]