[R-sig-Geo] Example of universal kriging with R/gstat in GRASS needed

Fri Apr 28 21:27:50 CEST 2006

On Fri, 28 Apr 2006, Thomas Adams wrote:

> Roger,
> 
> Your suggestion:
> 
> fullgrid(dem) <- FALSE
> 
> did turn dem into class type SpatialGridDataFrame, but when I tried:
> 
> z <- predict(UK_fit,newdata=dem)
> 
> I got an error:
> 
> Error in model.frame(... :
> invalid variable type.
> 
> I think I should restate the problem:
> 
> I have a file 'temps' which has class SpatialPointsDataFrame read from 
> GRASS 6.1, that looks like:
> 
> 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 I have a a file 'dem' which has class SpatialGridDataFrame which 
> just consists of grid of elevation values read from GRASS 6.1 using 
> dem<-readFLOAT6sp(). (Sorry, I know I'm repeating myself).
> 
> What I want to do is to use the grid of elevation values ('dem') as a 
> proxy in the spatial interpolation of the 'temp' values in my 'temps' 
> file that are located at the coordinates in parentheses(). Notice that 
> the temps file also has 'z' values of elevations. So, is this what you 
> already understood? Converting 'dem' to a SpatialPixelsDataFrame seemed 
> to only leave me with the grid locations and not the elevation values — 
> is this right.

What does:

summary(dem) 

say before and after doing

fullgrid(dem) <- FALSE?

Afterwards it should be a SpatialPixelsDataFrame with 

names(dem)

being "z". Saying summary(dem) will give you an idea of what is inside, 
str() should too.

Roger

PS. This is usually a one-off thing, once it works, you note down how, and 
then it just does from then on.

> 
> Thanks again for your help!
> 
> Regards,
> Tom
> 
> 
> Roger Bivand wrote:
> > On Fri, 28 Apr 2006, Thomas Adams wrote:
> >
> >   
> >> 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 think 
> >
> > fullgrid(dem) <- FALSE
> >
> > should make a SpatialPixelsDataFrame, but you'll have to make sure the 
> > name of the dem variable is the same as in the formula.
> >
> > Roger
> >
> >   
> >> 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
> >>>>
> >>>>
> >>>>     
> >>>>         
> >>>   
> >>>       
> >>
> >>     
> >
> >   
> 
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no