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: [EMAIL PROTECTED] VOICE: 937-383-0528 FAX: 937-383-0033 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo