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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
