[R-sig-Geo] search,grids,SA significance test

[Roger Bivand]

On Tue, 27 Feb 2007, Trevor Doerksen wrote:

> Hi,
> 
> I'm new to the mailing list and I see there is a large repository of 
> info in the archives. Is there any way to search it, to find answers to 
> my potentially redundant questions?

If going back 18 months is OK, try:

http://news.gmane.org/gmane.comp.lang.r.geo

Alternatively, the RSiteSearch() function inside R takes a search string, 
and searches R documentation including this mailing list, so that:

RSiteSearch("gstat krige")

hits this list a good deal.
> 
> I have some training in IDRISI's GIS and have used R for just over a 
> year. I've started using the gstat package in R and worked through some 
> of the excellent examples in 1) Pebesma and Bivand. 2005. S classes and 
> methods for spatial data: the sp package.  2) Pebesma. 2006. The meuse 
> data set: a tutorial for gstat R package. 3) Pebesma. 2004. 
> Multivariable geostatistics in S: the gstat package. Computers & 
> Geosciences. 30:683-691.
> 
> I have two specific questions:
> 1) I'm having trouble creating a grid of regularly spaced points that 
> gstat will accept for use with the krige function. Also, I don't know 
> how to mask a specific area (like the nice grid in the meuse data set) 
> of interest within that grid. Does anyone have examples/references not 
> listed above? Below is one attempt at some code which didn't work. NOTE: 
> My samples are irregularly spaced and I don't need predictions over the 
> entire rectangular grid.
> 
> #FIRST attempt
> x=seq(246000,606000,1000)
> length(x)    #361
> y=seq(4851000,5084000,1000)
> length(y)    #234
> pts.4.grd=expand.grid(x,y)    # resolution = 1km^2; size = 361*234 = 84474
> names(pts.4.grd)=c('x','y')
> str(pts.4.grd)    # dataframe
> grd.pts = SpatialPixels(SpatialPoints(pts.4.grd))
> str(grd.pts)    #SpatialPixels
> grd = as(grd.pts, "SpatialGrid")
> str(grd)    #X=Var1, Y=Var2
> 
> # variogram fit not shown. Trying 'universal kriging', with a N-S trend 
> removed.
> kr1.f2.v2=krige(bv~UTMN83+UTMN83*UTMN83, 
>     locations=~UTME83+UTMN83, data=all, grd, model=fit2.v2)

You could get to grd quicker by saying:

grdT <- GridTopology(cellcentre.offset=c(246000, 4851000), 
  cellsize=c(1000, 1000), cells.dim=c(361, 234))
grd <- SpatialGrid(grdT) # possibly setting proj4string

I would advise setting up a gstat object first, then predicting from it, 
to keep the steps apart for ease of debugging:

kr1.f2.v2_input <- gstat(bv ~ UTMN83 + I(UTMN83*UTMN83), data=all, 
  model=fit2.v2)

assuming that bv and UTMN83 are columns of the all object, and that the 
all object is a SpatialPointsDataFrame. Please do note that 
I(UTMN83*UTMN83) is of the order of 2.35322e+13, and that this may render 
the UK part of your model unworkable numerically. 

kr1.f2.v2 <- predict(kr1.f2.v2_input, newdata=grd)

will now probably fail because grd does not have a column called UTMN83.

So I suggest you first try setting the variogram model on a gstat object 
including a numerically secure quadratic trend (both directions) - see the 
degree= argument to gstat() - 

quad_obj <- gstat(bv ~ 1, data=all, degree=2)

then get the variogram for quad_obj, fit the variogram, make a new gstat 
object, then predict. When that works, make it more complicated by taking 
out the UTMN83 coordinates as a separate local variable (scaled in the 
same way in all and grd), and it should just work. The source of problems 
here is trying to do too many things in a single step.

> #doesn't work with the following error. Grid and sample data not same size?
> #Error in gstat.formula.predict(d$formula, newdata, na.action = 
> na.action) :
> #        NROW(locs) != NROW(X): this should not occur
> #In addition: Warning messages:
> #1: 'newdata' had 84474 rows but variable(s) found have 440 rows
> #2: 'newdata' had 84474 rows but variable(s) found have 440 rows
> 
> 2) Is there a test to compute the significance of spatial 
> autocorrelation in a variogram? The biological signal I'm modelling in 
> my variograms is very weak, with a huge nugget. How large can the nugget 
> be before it negates the spatially specific weights for kriging, thus 
> rendering the prediction close to an inverse distance weighting scheme?
> 

This is a different question, but the predictions would only be like a 
regional disc if the maximum number of neighbours or the maximum distance 
was set.

Roger

> Thanks,
> 
> Trevor.
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

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