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

Fri Apr 28 21:13:12 CEST 2006

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.

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
>>>>
>>>>
>>>>     
>>>>         
>>>   
>>>       
>>
>>     
>
>   

-- 
Thomas E Adams
National Weather Service
Ohio River Forecast Center
1901 South State Route 134
Wilmington, OH 45177

EMAIL:	thomas.adams at noaa.gov

VOICE:	937-383-0528
FAX:	937-383-0033