Help with universal kriging using gstat

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Help with universal kriging using gstat

Thomas Adams-2
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
[hidden email] (personal)
[hidden email] (work)

1 (513) 739-9512 (cell)

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Help with universal kriging using gstat

Roger Bivand
Administrator
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: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Help with universal kriging using gstat

Thomas Adams-2
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 <[hidden email]> 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: [hidden email]
> 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
[hidden email] (personal)
[hidden email] (work)

1 (513) 739-9512 (cell)

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo