[STATSGRASS] kriging, R and GRASS

[Roger Bivand]

On Sun, 16 Jan 2005, A Horacio Samaniego wrote:

> I haven't use GRASS for a little while and seem to have forgotten how the 
> tricks work...
> 
> I wonder is there is a way to do the kriging interpolation in some 
> efficient way. I remember that you have to actually create a grid in R to 
> export to GRASS later on... I also, remember that this is terribly slow 
> and unstable. Is this still the case?

Using the UTM maas location, (see:

http://grass.itc.it/statsgrass/grass_r_interface.html#Using_GRASS_R_Maas)

I get:

> system.time(G <- gmeta())
[1] 0 0 0 0 0
> summary(G)
Data from GRASS 5.0 LOCATION maas with 259 columns and 232 rows;
UTM, zone: 32
The west-east range is: 269870, 272460,
and the south-north: 5650610, 5652930;
West-east cell sizes are 10 units,
and south-north 10 units.
> system.time(maas <- sites.get(G, "ex.utm.maas"))

## your sites file name may vary, in the tar archive it is called "Zn" and 
## only has the easting, northing, and single variable, use str() to see 
## what they are called

[1] 0.04 0.01 0.05 0.00 0.00
> str(maas)
`data.frame':   98 obs. of  15 variables:
 $ east   : num  271674 271930 272146 272391 271789 ...
 $ north  : num  5652772 5652738 5652728 5652704 5652714 ...
 $ elev   : num  8.06 7.51 6.98 7.8 7.94 7.98 8.8 7.9 7.74 6.68 ...
 $ id     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ x      : num  1637 1894 2110 2356 1755 ...
 $ y      : num  2651 2630 2630 2618 2599 ...
 $ d.river: num  10 170 340 550 100 10 60 10 10 70 ...
 $ Cd     : num  8.2 2.4 3 1.6 4.2 17 2.4 12 9.4 10.9 ...
 $ Cu     : num  47 32 32 27 51 128 47 117 104 90 ...
 $ Pb     : num  191 102 97 82 281 405 297 654 482 541 ...
 $ Zn     : num  812 298 321 213 746 ...
 $ LOI    : num  11.1 1.4 1.6 3.1 5.1 12.3 10 16.5 13.9 10.2 ...
 $ Fldf   : num  1 1 1 1 1 1 2 1 1 1 ...
 $ Soil   : num  1 2 2 2 2 1 1 1 1 1 ...
 $ Zn.o   : Factor w/ 5 levels "crisis","high",..: 2 4 4 4 2 1 2 1 1 1 ...
 - attr(*, "nncols")= Named int  3 0 10 1
  ..- attr(*, "names")= chr  "dims" "cattype" "dbl.att" "str.att"

> system.time(fit <- Tps(cbind(maas$east, maas$north), maas$Zn))
[1] 0.42 0.05 0.47 0.00 0.00
Warning message:
GCV search gives a minumum at the endpoints of the grid search in: 
Krig.find.gcvmin(info, lambda.grid, gcv.grid$"-Log Profile",  

> system.time(out.p<-predict.surface(fit, grid.list=list(X1=G$xseq, 
+  X2=G$ryseq)))
[1] 4.60 0.92 5.51 0.00 0.00

## using X2=G$ryseq to reverse the row order

> system.time(rast.put(G, "fldsTpsz0", c(out.p$z)))
[1] 0.11 0.03 0.15 0.00 0.00

(Timings on RHEL 3, Xeon 1.5GHz, R-devel, R GRASS package as released on 
CRAN (0.2-22), and GRASS 5.4, 5.0.* will do the same if that is what you 
are running)

I think your recollection is from several years ago, before the compiled 
interface was written. Please let the list know how you get on!

Best wishes,

Roger 

> 
> I have gone trhough the geostat litterature in R and GRASS and are getting 
> comfortable with the fields package to do all my variograms. However I've 
> come to the stage where I need to predict ....
> 
> What is the best way to do this so that I can end up with a GRASS raster 
> map? I'll  appreciate any hints on how to go about this.
> 
> 
> thank you very much!
> 
> H
> 
> Horacio Samaniego
> Department of Biology
> University of New Mexico
> Albuquerque, NM. 87131
> 
> _______________________________________________
> statsgrass mailing list
> statsgrass at grass.itc.it
> http://grass.itc.it/mailman/listinfo/statsgrass
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no