[GRASSLIST:10117] Re: R, kriging and grass6

[Roger Bivand]

On Wed, 1 Feb 2006, ivan marchesini wrote:

> Dear Roger, dear grass/R users...
> we would like to execute a kriging process of points data only inside a
> well defined area.
> we have realised a raster map that we would use as MASK inside R.
> how can we impose to the krige command to perform only into the points
> contained in the MASk???
> we have imported the raster map in R using readFLOAT6sp and the
> resulting map (named "structure") presents some points having value=11
> and a lot of points having value=NA.
> then we have created a
> 
> mask_struttura<-SpatialGridDataFrame(grd,data=structure,,proj4string=CRS(G$proj4))
> 
> later we have used the gstat krige command to perform the prediction,
> defining the parameter newdata=mask_structure
> but we obbtain this error:
> NROWS(locs) != NROWS(X)
> 
> using simply all the region to perform the prediction 
> mask_SG <- SpatialGridDataFrame(grd, data=list(k=rep(1, G$cols*G
> $rows)),proj4string=CRS(G$proj4))
> as you suggested, all is right...  
> where are we wronging???  :-(

Trying too hard? This should be simple enough, by converting the full 
SpatialGrid with NA values to SpatialPixels, and using them to predict to?

Using the meuse data as before:

> library(spgrass6) 
> meuse <- getSites6sp("meuse") 
Exporting 155 points/lines...
 100%
155 features written
> library(gstat)
> cvgm <- variogram(zinc ~ 1, locations=meuse, width=100, cutoff=1000)
> efitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100,
+   nugget=1))
> ffreq <- rast.get6("ffreq1")
# The GRASS raster layer has missing values outside the study area
> ffreq_sp <- as(ffreq, "SpatialPixelsDataFrame")
# so this adds coordinates and drops those that are NA in the layer
> plot(ffreq_sp)
> OK_pred <- krige(zinc ~ 1, locations=meuse, newdata=ffreq_sp, 
+ model=efitted)
> spplot(OK_pred, "var1.pred")

Isn't that nice? That's why Edzer and the gang have spent so long working 
on the sp classes, to make life easier. And writing out to GRASS will put 
the missing values back in the right places too.


> thank you 
> 
> Ivan
> 
> PS: we have started a wiki for the grass/R kriging problem in the grass
> wiki.. as requested from someone.. 
> 

Good idea, please add this.

Roger

PS: A PS from me - rast.get6() and rast.put6() will work faster and not be 
so dependent on the Arc ASCII grid format if we could write either GeoTiff 
or at least the format used by r.in.bin from R - doing what r.out.bin and 
r.in.bin do would be fine. GDAL is OK, but using the GDAL plugin for 
r.out.gdal still needs firm nerves, and it doesn't respect the current 
region. So if anyone could help write the equivalents to r.in/out.bin in 
R, things would get even better!


>  
> 
> 
> 
> 
> 
> 
> 
> 
> Il giorno mer, 25/01/2006 alle 10.47 +0100, Roger Bivand ha scritto:
> > On Wed, 25 Jan 2006, ivan marchesini wrote:
> > 
> > > Dear Roger...
> > > we have tried to use your suggestions...
> > > we have tried using ours data...
> > > all seems to work fine, but at the end (writeRast6sp) we have a strange
> > > error regarding the dimension of the cells:
> > > 
> > > error in writeAsciGrid (.........)
> > > Asciigrid does not support grids with non square cells
> > > 
> > > but this is not clear because we have used Gmeta6() to collect the
> > > location data of grass and before we have set "g.region res=20".
> > > 
> > > now I have seen that "g.region -p" gives me a resolution of
> > > ns=20.0003 (or similar)
> > > eo=19.99998 (or similar)
> > > 
> > > do you think is possible that this is the problem???
> > 
> > Yes, Arc ASCII grids have to have exactly square cells. I guess you can 
> > also force it by using:
> > 
> > cellsize=c(20, 20)
> > 
> > in the GridTopology() call.
> > 
> > Roger
> > 
> > > now we are trying to use a perfect integer region (ns=20 and eo=20)
> > > I will let you know the results....
> > > 
> > > best wishes
> > > 
> > > 
> > > 
> > > Il giorno lun, 23/01/2006 alle 21.09 +0100, Roger Bivand ha scritto:
> > > > On Mon, 23 Jan 2006, ivan marchesini wrote:
> > > > 
> > > > > Dear Roger...
> > > > > thank you for your previous answer... 
> > > > > can you fast describe me how to create a SpatialGrid object starting
> > > > > from the data returned from Gmeta6()???
> > > > > only few hints...  :-)
> > > > > I only need this to perform the kriging process, because I have already
> > > > > generated the variogram and I only need the grid over which to perform
> > > > > the prediction.....
> > > > > thank you very much...
> > > > > ivan
> > > > 
> > > > This is the dense version, feedback welcome!
> > > > 
> > > > library(spgrass6) # load package running within GRASS 6
> > > > meuse <- getSites6sp("meuse") # get vector points as 
> > > > SpatialPointsDataFrame
> > > > class(meuse)
> > > > G <- gmeta6() # get region
> > > > grd <- GridTopology(cellcentre.offset=c(G$west+(G$ewres/2), 
> > > >   G$south+(G$nsres/2)), cellsize=c(G$ewres, G$nsres),
> > > >   cells.dim=c(G$cols, G$rows))
> > > > mask_SG <- SpatialGridDataFrame(grd, data=list(k=rep(1, G$cols*G$rows)),
> > > >   proj4string=CRS(G$proj4)) # create a SpatialGridDataFrame
> > > > class(mask_SG)
> > > > cvgm <- variogram(zinc~1, loc=~x+y, data=meuse, width=100, cutoff=1000)
> > > > efitted <- fit.variogram(cvgm, vgm(psill=1, model="Exp", range=100, 
> > > >   nugget=1))
> > > > OK_pred <- krige(id="OK_pred", formula = zinc ~ 1, locations = ~ x + y,
> > > >   data = as(meuse, "data.frame"), newdata=mask_SG, model=efitted)
> > > >   # make the kriging prediction
> > > > names(OK_pred)
> > > > writeRast6sp(OK_pred, "OK_pred", zcol = "OK_pred.pred", NODATA=-9999)
> > > > 
> > > > Roger
> > > > 
> > > > 
> > > > > 
> > > > > 
> > > > > Il giorno lun, 23/01/2006 alle 13.14 +0100, Roger Bivand ha scritto:
> > > > > > On Mon, 23 Jan 2006, ivan marchesini wrote:
> > > > > > 
> > > > > > > Hi to all...
> > > > > > > we are triyng to use kriging using R and grass6,
> > > > > > > We have seen a lot of documentation about R and GRASS5 but not so many
> > > > > > > about grass6...
> > > > > > > the problem is that we are trying to do an interpolation using the
> > > > > > > functions gmeta6 and kriging...
> > > > > > > if we are not wrong, the R library GRASS is developed for grass5 and
> > > > > > > inside there is the function krige.G that use the location data acquired
> > > > > > > using the gmeta() function...
> > > > > > > the problem is that we have found the gmeta6() function but not the
> > > > > > > krige.G6 function and obvously the krige.G doesn't work using the
> > > > > > > location data of grass6 (obtained using gmeta6)  
> > > > > > > 
> > > > > > > how can we simply solve this problem....
> > > > > > 
> > > > > > krige.G was a nasty hack, and should best be forgotten. The direct route 
> > > > > > is to use the sp classes in R and implemented for GRASS in spgrass6 - see 
> > > > > > GRASS News for the GRASS side and R News for the R side, combined with a 
> > > > > > proper R geostats package. Of those available, gstat is tightly bound to 
> > > > > > sp classes, so as far as I know, the only bit that needs doing by hand is 
> > > > > > to create a GridTopology object, or a SpatialGrid object, from the data 
> > > > > > returned by gmeta6() to pass to the kriging prediction function.
> > > > > > 
> > > > > > If you like, we can iterate to a working example from passing the vector 
> > > > > > points to R from GRASS, doing the modelling and kriging predictions, and 
> > > > > > passing the prediction rasters back to GRASS, but I'd like input from 
> > > > > > users to make the description sound.
> > > > > > 
> > > > > > Best wishes,
> > > > > > 
> > > > > > Roger
> > > > > > 
> > > > > > > 
> > > > > > > thank you
> > > > > > > 
> > > > > > > Ivan
> > > > > > > 
> > > > > > > 
> > > > > > > 
> > > > > > > 
> > > > > > 
> > > > > 
> > > > 
> > > 
> > 
> 

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