Example of universal kriging with R/gstat in GRASS needed

Roger,

This got me further along, but I am encountering a problem with:

z <- predict(UK_fit, newdata=BMcD_SPx)

The gstat step works for me, where I have:

UK_fit<-gstat(formula=temps$temp~dem,data=temps,model=efitted)

temps has class SpatialPointsDataFrame:

              coordinates cat      x     y    z temp name
1     (-341460, -2154.42)   1 -90.05 38.90  166   63  ALN
2       (-198769, 301388)   2 -88.47 41.77  215   67  ARR
3       (-334899, -40321)   3 -89.95 38.55  140   66  BLV
4       (-240028, 163910)   4 -88.92 40.48  268   69  BMI
5       (-187957, 114806)   5 -88.27 40.04  229   64  CMI
6     (-351730, -37305.9)   6 -90.15 38.57  126   65  CPS
7      (-242424, 98244.7)   7 -88.92 39.87  204   66  DEC
8       (-179844, 315889)   8 -88.24 41.91  232   69  DPA
9     (-136093, -24538.2)   9 -87.61 38.76  131   68  LWV
10     (-278964, -126152)  10 -89.25 37.78  125   66  MDH
11      (-140792, 302011)  11 -87.75 41.79  187   73  MDW
12      (-364737, 274189)  12 -90.51 41.45  180   73  MLI
13     (-190503, 54493.9)  13 -88.28 39.48  219   64  MTO

and dem has class SpatialGridDataFrame and just consists of grid values.

I tried to create a SpatialPixelsDataFrame for predict(), but with (for 
example):

m = SpatialPixelsDataFrame(points=meuse.grid[c("x","y")],data=meuse.grid)

I have nothing like meuse.grid, so this does not work. I can use 
image(dem), which produces a plot of elevation values. My point is that 
meuse.grid and my dem files have very different structures.

I'm not sure where to go to from here.

Regards,
Tom


Roger Bivand wrote:

On Thu, 27 Apr 2006, Thomas Adams wrote:

  

List:

I can not seem to work out the syntax for using R/gstat within a GRASS 
6.1 session to do universal kriging. I have a DEM (elevation data on a 
grid) and point data for temperature; theoretically, the temperatures 
should relate to elevation. So, I am trying to spatially interpolate the 
temperature data based on the elevations at the grid points. How do I 
setup the gstat command in R/gstat (and using spgrass6, of course)? I 
have no trouble reading in my elevation data (DEM) from GRASS and I have 
no problem doing ordinary kriging of my temperature data using 
GRASS/R/gstat.
    

What do the data look like? Do you have temperature and elevation at the
observation points and elevation over the grid? If temperature is the 
variable for which you want to interpolate, then the formula argument in 
the gstat() function would be temp ~ elev, data=pointsdata (if a 
SpatialPointsDataFrame no need for location= ~ x + y). Then the predict() 
step would need a SpatialGridDataFrame object as newdata, with elev as 
(one of) the columns in the data slot.

An example for the Meuse bank data in Burrough and McDonnell:

cvgm <- variogram(Zn ~ Fldf, data=BMcD, width=100, cutoff=1000)
uefitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100, 
  nugget=1))
UK_fit <- gstat(id="UK_fit", formula = Zn ~ Fldf, data = BMcD, 
  model=uefitted)
z <- predict(UK_fit, newdata=BMcD_SPx)

where BMcD_SPx is a SpatialPixelsDataFrame (the grid has ragged edges) 
with flood frequencies in Fldf (actually a factor, but works neatly).

Hope this helps,

Roger

  

Regards,
Tom