[GRASS-stats] another gstat error

Sat Jun 20 05:01:34 EDT 2009

On Sat, 20 Jun 2009, Jarek Jasiewicz wrote:

>
>
> Roger Bivand pisze:
>> On Fri, 19 Jun 2009, Jarosław Jasiewicz wrote:
>> 
>>> Hi
>>> Another gstat error:
>>> 
>>> var.krigged=krige(z~x+y,gr,grid,model=var.fitted)
>>> [using universal kriging]
>>> 
>>> "chfactor.c", line 130: singular matrix in function LDLfactor()
>>> Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim, 
>>> :
>>> LDLfactor
>> 
>> Jarek:
>> 
>> Shall I repeat what I wrote this morning? 
> ????

My reply is in the list archives and nabble, so I don't know why you 
didn't get it.

>> Did you actually try omitting the x+y in the formula? 
> Yes, before I send post I tried firrst z~(x+y)^2 and if it failed I used z~1 
> to try without drift. The same error ocurs

Googling 'gstat "chfactor.c", line 130: singular matrix in function 
LDLfactor', I got to:

http://osdir.com/ml/lang.r.geo/2006-07/msg00049.html

which says:

HINT: Read the manual at http://www.gstat.org/ ;
look for: Trouble shooting -> Error messages -> From meschach

"chfactor.c", line 130: singular matrix in function LDLfactor()
Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim, :
matrix library error: gstat: matrix library error: singular matrix
---end R output

after consulting the gstat manual I found the following reference:

"These two error messages may occur (a) during simulation, when an
observation falls almost exactly at a simulation location (b) when two
observations occur at identical location occur and noaverage was defined 
in
its data definition. Solution to (a): increase the value of zero, to (b):
remove the noaverage. Read also the next section."

My input file was converted from lat/long to UTM using spTransform. The
sampling density is quite high and it is possible that after conversion
from latlong to UTM because of rounding some measurements fall on 
identical
coordinates. Is this a possible explanation for the error?

The GSTAT online manual suggests to remove "noaverage" or increase "zero"
to overcome the problem. However, I have not figured out how to do this
when running gstat under R.

So it looks as though you may have points that are seen as identical in 
your data set.

Hope this helps,

Roger

>> Including a drift in untransformed coordinates is very likely to lead to 
>> numerical chaos. If z~1 works, then the smoking gun is the x+y, isn't it? 
>> If you want them, transform them (offset with for example elide() in 
>> maptools and transform back afterwards).
>> 
>> Roger
>> 
>>> 
>>> I can't imagine the reason
>>> 
>>> regards
>>> Jarek
>>> 
>>> here is structure of my gr spatialPointsDataFrame used in analysis:
>>> 
>>> str(gr)
>>> Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
>>> ..@ data       :'data.frame':    2201 obs. of  5 variables:
>>> .. ..$ cat: num [1:2201] 2 3 4 5 6 7 8 9 10 11 ...
>>> .. ..$ naz: Factor w/ 874 levels "1000","1001",..: 874 1 2 3 4 5 6 7 8 9 
>>> ...
>>> .. ..$ x  : num [1:2201] 360175 360176 360178 360169 360167 ...
>>> .. ..$ y  : num [1:2201] 512736 512736 512733 512719 512712 ...
>>> .. ..$ z  : num [1:2201] 86.5 87.4 87.7 86.9 86 ...
>>> ..@ coords.nrs : num(0)
>>> ..@ coords     : num [1:2201, 1:2] 360175 360176 360178 360169 360167 ...
>>> .. ..- attr(*, "dimnames")=List of 2
>>> .. .. ..$ : NULL
>>> .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
>>> ..@ bbox       : num [1:2, 1:2] 359907 512481 360492 512778
>>> .. ..- 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 " +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 
>>> +x_0=500000 +y_0=-5300000 +no_defs +a=6378137 +rf=298.257222101 
>>> +towgs84=0.000,0.000,"| __truncated__
>>> 
>>> and newdata grid:
>>> str(grid)
>>> Formal class 'SpatialGridDataFrame' [package "sp"] with 6 slots
>>> ..@ data       :'data.frame':    131352 obs. of  1 variable:
>>> .. ..$ region2: int [1:131352] NA NA NA NA NA NA NA NA NA NA ...
>>> ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
>>> .. .. ..@ cellcentre.offset: Named num [1:2] 359902 512478
>>> .. .. .. ..- attr(*, "names")= chr [1:2] "x" "y"
>>> .. .. ..@ cellsize         : num [1:2] 1.001 0.998
>>> .. .. ..@ cells.dim        : int [1:2] 421 312
>>> ..@ grid.index : int(0)
>>> ..@ coords     : num [1:2, 1:2] 359902 360322 512478 512788
>>> .. ..- attr(*, "dimnames")=List of 2
>>> .. .. ..$ : NULL
>>> .. .. ..$ : chr [1:2] "x" "y"
>>> ..@ bbox       : num [1:2, 1:2] 359901 512477 360323 512789
>>> .. ..- attr(*, "dimnames")=List of 2
>>> .. .. ..$ : chr [1:2] "x" "y"
>>> .. .. ..$ : chr [1:2] "min" "max"
>>> ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>>> .. .. ..@ projargs: chr " +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 
>>> +x_0=500000 +y_0=-5300000 +no_defs +a=6378137 +rf=298.257222101 
>>> +towgs84=0.000,0.000,"| __truncated__
>>> _______________________________________________
>>> grass-stats mailing list
>>> grass-stats at lists.osgeo.org
>>> http://lists.osgeo.org/mailman/listinfo/grass-stats
>>> 
>> 
>

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