[R-sig-Geo] Example of universal kriging with R/gstat in GRASS
Roger Bivand
Roger.Bivand at nhh.no
Sat May 13 10:16:47 CEST 2006
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
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>
> >>>>>
> >>>>>
> >>>>
> >>>>
> >>>
> >>>
> >>
> >>
> >
> >
>
>
>
--
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
More information about the R-sig-Geo
mailing list