[STATSGRASS] kriging, R and GRASS

Mon Jan 17 08:37:44 EST 2005

On Mon, 17 Jan 2005, Roger Bivand wrote:

... and also wanted show that rastput() would scale well to much finer
grids - which it does - and found a bug in rastput() and other functions
that has led to the layer being written to the GRASS database with the
wrong window if system("g.region ... ") is used. This is the same bug that
affected gmeta() until Bodo Ahrens found it 18 months ago. The problem was
that the GRASS library function tries to be too clever (I think disk
access was very slow when it was written), and stores a local version of
the data in memory in a cache of sorts. It does not check, though, whether
the data on disk have changed - it just returns the cached values! Using a 
lower-level function, and conditionally altering G__init_window () seems 
to have fixed it. The new version should reach CRAN soon.

Returning to the example, this is going from 10m x 10m cells to 5m x 5m:

> system("g.region rast=fldsTpsz_5z")
> G <- gmeta()
> summary(G)
Data from GRASS 5.0 LOCATION maas with 518 columns and 464 rows;
UTM, zone: 32
The west-east range is: 269870, 272460,
and the south-north: 5650610, 5652930;
West-east cell sizes are 5 units,
and south-north 5 units.
> system.time(out.p<-predict.surface(fit, grid.list=list(X1=G$xseq, 
+ X2=G$ryseq)))
[1] 18.87  4.50 23.71  0.00  0.00
> system.time(rast.put(G, "fldsTpsz_5x", c(out.p$z)))
[1] 0.30 0.11 0.41 0.00 0.00

So the question was rather useful, too.

Roger

> 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