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

Thomas Adams Thomas.Adams at noaa.gov
Mon May 15 21:47:14 CEST 2006


Roger,

That did it! The sessionInfo() command identified that I had older 
versions than what were available on CRAN, namely:

lattice -- 0.13-8
gstat -- 0.9-29 (was older: 0.9-17)
sp -- 0.8-14 (was older: 0.8-8)

On my Mac, I have:

lattice -- 0.11-8
gstat -- 0.9-21
sp -- 0.7-10

which works…

BTW, I also updated to R v. 2.3 on Linux. Thanks again for your help!

Regards,
Tom



Roger Bivand wrote:
> Thomas:
>
> Please do:
>
> sessionInfo()
>
> on both machines, and see if any of the package versions are different.
>
> It often also helps to name function arguments explicitly, to be sure that 
> they are being understood correctly inside the functions. Depending on the 
> class of the data= arugment, as far as I recall different versions of the 
> gstat package treat the locations= argument differently.
>
> Roger
>
>
> On Fri, 12 May 2006, Thomas Adams wrote:
>
>   
>> Roger,
>>
>> I have made considerable progress in what I was trying to do with gstat 
>> & Universal Kriging using R. I was fortunate enough to find a detailed 
>> example of what I was trying to do: http://spatial-analyst.net/RKguide.php.
>>
>> I exported from GRASS my elevation grid, "gtopo30.dem", and temperature 
>> point file, "temps.txt", to ascii files to assure I followed a parallel 
>> path with my data.
>>
>> So, with my data, following Tomislav Hengl's example, I have:
>>
>>  > temps<-read.delim("temps.txt",sep=" ")
>>  > summary(temps)
>>
>> cat lon lat z T name
>> Min. : 1.00 Min. :-90.05 Min. :35.82 Min. : 99.0 Min. :61.00 AGC : 1
>> 1st Qu.: 31.75 1st Qu.:-86.63 1st Qu.:38.27 1st Qu.: 197.0 1st Qu.:64.00 
>> AID : 1
>> Median : 60.50 Median :-84.06 Median :40.09 Median : 255.5 Median :66.00 
>> ALN : 1
>> Mean : 60.37 Mean :-84.14 Mean :39.86 Mean : 314.2 Mean :66.83 AOH : 1
>> 3rd Qu.: 89.25 3rd Qu.:-81.47 3rd Qu.:41.58 3rd Qu.: 360.0 3rd Qu.:69.00 
>> AOO : 1
>> Max. :118.00 Max. :-78.32 Max. :42.49 Max. :1155.0 Max. :73.00 ARB : 1
>> (Other):110
>> X Y
>> Min. :-341460 Min. :-342057
>> 1st Qu.: -53847 1st Qu.: -71919
>> Median : 163999 Median : 119784
>> Mean : 154474 Mean : 99916
>> 3rd Qu.: 371950 3rd Qu.: 282632
>> Max. : 632335 Max. : 398186
>>
>>  > library(sp)
>>  > 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)
>>
>> Object of class SpatialPointsDataFrame
>> Coordinates:
>> min max
>> X -341459.8 632334.6
>> Y -342056.9 398185.9
>> Is projected: NA
>> proj4string : [NA]
>> Number of points: 116
>> Data attributes:
>> gtopo30.dem
>> Min. : 115.3
>> 1st Qu.: 196.9
>> Median : 245.4
>> Mean : 306.9
>> 3rd Qu.: 331.0
>> Max. :1064.5
>>
>>  > temps$gtopo30.dem=dem.ov$gtopo30.dem
>>  > library(lattice)
>>  > plot(T~gtopo30.dem, as.data.frame(temps))
>>  > abline(lm(T~gtopo30.dem, as.data.frame(temps)))
>>  > library(gstat)
>>
>>  > 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 
>> TEMPERATURE")
>>
>>
>> Which works perfectly on my Macintosh running Mac OS X 10.4 and using R 
>> 2.2.1. (see attachment, temperatures in deg. F) However, following the 
>> *identical* steps with the identical data on Linux, at the step:
>>
>> temps_uk<-krige(T~gtopo30.dem,temps,dem, vgm_temps_r)
>>
>> I get the error:
>>
>> Error in eval(expr, envir, enclos) : object "gtopo30.dem" not found
>>
>> This has me baffled; any thoughts? I could send you my files if you 
>> would like to see what happens for you…
>>
>> Regards,
>> Tom
>>
>> BTW, the grid spacing on my DEM is coarse (9 km) and I will probably do 
>> my final analyses at 1-km.
>>
>>
>> Roger Bivand wrote:
>>     
>>> 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
>>>>>>>>
>>>>>>>>
>>>>>>>>     
>>>>>>>>         
>>>>>>>>             
>>>>>>>>                 
>>>>>>>   
>>>>>>>       
>>>>>>>           
>>>>>>>               
>>>>>>     
>>>>>>         
>>>>>>             
>>>>>   
>>>>>       
>>>>>           
>>>>     
>>>>         
>>>   
>>>       
>>
>>     
>
>   


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




More information about the R-sig-Geo mailing list