[STATSGRASS] gioo78@libero.it

Edzer J. Pebesma e.pebesma at geo.uu.nl
Thu Aug 3 12:27:10 EDT 2006


Well, you ignored a warning message after the fit.variogram() call; 
verifying its result might help. Otherwise, I have no idea, except that 
I would try leaving out the periodic component from the variogram model.
--
Edzer

george wrote:
> Hello all!
> I have this problem with local kriging in Grass6.0 interface R:
>
> Display my region in location XY
> projection: 0 (x,y)
> zone:       0
> north:      1032
> south:      687
> west:       248.22524977
> east:       395.35876476
> nsres:      1
> ewres:      1.00090827
> rows:       345
> cols:       147
>
>
> GRASS 6.0.0 (XY):~ > R
>   
>> library(spgrass6)
>>     
> Loading required package: sp
> Loading required package: maptools
> Loading required package: foreign
> Loading required package: rgdal
> Loading required package: abind
> Loading required package: pixmap
> Geospatial Data Abstraction Library extensions to R successfully loaded
>   
>> library(gstat)
>>     
> # import my set data x|y|z 
>   
>> impvett<-read.table(file="/home/giorgio/Desktop/s50_p2_ba_bb.txt",
>>     
> sep="|")
>   
>> str(impvett)
>>     
> `data.frame':   1605 obs. of  3 variables:
>  $ V1: num  249 251 253 255 257 ...
>  $ V2: num  1032 1032 1032 1032 1032 ...
>  $ V3: num  -0.117 -0.230  0.897  2.252  2.912 ...
>   
>> coordinate<-data.frame(impvett$V1,impvett$V2)
>>  data<-data.frame(impvett$V3)
>>  rilievo<-SpatialPointsDataFrame(coordinate,data,proj4string =
>>     
> CRS(as.character(NA)),match.ID = TRUE)
>   
>> str(rilievo)
>>     
> Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
>   ..@ data       :Formal class 'AttributeList' [package "sp"] with 1
> slots
>   .. .. ..@ att:List of 1
>   .. .. .. ..$ impvett.V3: num [1:1605] -0.117 -0.230  0.897  2.252
> 2.912 ...
>   ..@ coords.nrs : num(0)
>   ..@ coords     : num [1:1605, 1:2] 249 251 253 255 257 ...
>   .. ..- attr(*, "dimnames")=List of 2
>   .. .. ..$ : NULL
>   .. .. ..$ : chr [1:2] "impvett.V1" "impvett.V2"
>   ..@ bbox       : num [1:2, 1:2]  249  687  395 1032
>   .. ..- attr(*, "dimnames")=List of 2
>   .. .. ..$ : chr [1:2] "impvett.V1" "impvett.V2"
>   .. .. ..$ : chr [1:2] "min" "max"
>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>   .. .. ..@ projargs: chr NA
>
> #import region setting and create grid
>   
>> G<-gmeta6()
>> 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)))
>   
>> str(mask_SG)
>>     
> Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
>   ..@ data       :Formal class 'AttributeList' [package "sp"] with 1
> slots
>   .. .. ..@ att:List of 1
>   .. .. .. ..$ k: num [1:50715] 1 1 1 1 1 1 1 1 1 1 ...
>   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3
> slots
>   .. .. ..@ cellcentre.offset: num [1:2] 249 688
>   .. .. ..@ cellsize         : num [1:2] 1 1
>   .. .. ..@ cells.dim        : int [1:2] 147 345
>   ..@ grid.index : int(0)
>   ..@ coords     : num [1:2, 1:2]  249  395  688 1032
>   .. ..- attr(*, "dimnames")=List of 2
>   .. .. ..$ : NULL
>   .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
>   ..@ bbox       : num [1:2, 1:2]  248  687  395 1032
>   .. ..- attr(*, "dimnames")=List of 2
>   .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
>   .. .. ..$ : chr [1:2] "min" "max"
>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>   .. .. ..@ projargs: chr NA
>
> # variogram
>   
>> cvgm<- variogram(impvett.V3~1,location=rilievo,cutoff=300,width=1)
>> efitted <-
>>     
> fit.variogram(cvgm,vgm(29.5,model="Exp",range=12,nugget=0,add.to=vgm(0.7,"Per",118,0)))
> Warning: singular model in variogram fit
>
> #ordinary kriging all data 
>   
>>  kr<-krige(impvett$V3~1,location=rilievo, newdata=mask_SG
>>     
> model=efitted)
> [using ordinary kriging]
>   
>> # ok!!!!
>>     
>
> # !!!!!!!!!!!!!!!!ORDINARY LOCAL KRIGING!!!!!!!!!!!!!!!!!!!!
>
> #case A) 
>   
>> kr<-krige(impvett$V3~1,location=rilievo, newdata=mask_SG,
>>     
> model=efitted,maxdist=50)
> [using ordinary kriging]
>
> "chfactor.c", line 130: singular matrix in function LDLfactor()
>
> gstat caught an error that occurred in the matrix library,
> the reason for it was: singular matrix
>
> HINT: Read the manual at http://www.gstat.org/ ;
> look for: Trouble shooting -> Error messages -> From meschach
>
> Error in predict.gstat(g, newdata = newdata, block = block, nsim =
> nsim,  :
>         matrix library error: gstat: matrix library error: singular
> matrix
>
> #case B)
>   
>> kr<-krige(impvett$V3~1,location=rilievo, newdata=mask_SG,
>>     
> model=efitted,maxdist=300)
> [using ordinary kriging]
>
> "chfactor.c", line 130: singular matrix in function LDLfactor()
>
> gstat caught an error that occurred in the matrix library,
> the reason for it was: singular matrix
>
> HINT: Read the manual at http://www.gstat.org/ ;
> look for: Trouble shooting -> Error messages -> From meschach
>
> Error in predict.gstat(g, newdata = newdata, block = block, nsim =
> nsim,  :
>         matrix library error: gstat: matrix library error: singular
> matrix
>
> #case C)
>   
>> kr<-krige(impvett$V3~1,location=rilievo, newdata=mask_SG,
>>     
> model=efitted,nmax=50) 
> [using ordinary kriging]
>
> "chfactor.c", line 130: singular matrix in function LDLfactor()
>
> gstat caught an error that occurred in the matrix library,
> the reason for it was: singular matrix
>
> HINT: Read the manual at http://www.gstat.org/ ;
> look for: Trouble shooting -> Error messages -> From meschach
>
> Error in predict.gstat(g, newdata = newdata, block = block, nsim =
> nsim,  :
>         matrix library error: gstat: matrix library error: singular
> matrix
>
>
> Where is the problem? 
> I have followed examples on package gstat: with data meuse function
> krige() work correctly!
>
> Thanks 
> George
>
>
>
>
>
>
>
> _______________________________________________
> statsgrass mailing list
> statsgrass at grass.itc.it
> http://grass.itc.it/mailman/listinfo/statsgrass
>   




More information about the grass-stats mailing list