[R-sig-Geo] negative range from fit.variogram through autofitVariogram

[Mark Connolly]

I am using autofitVariogram during the process of interpolating a large 
set of daily observations through a volume.  Each volume is decomposed 
into 2D layers prior to selecting a model to use for interpolation.  I 
made it through 2010 interpolations and then ran into a failed 
interpolation when the best model selected by autofitVariogram had a 
negative range.  This was rejected by the krige function.  I see mention 
of negative sills but not of negative ranges.

It appears that autofitVariogram is having some issues with the trial 
arguments sent to fit.variogram.  This is repeatable.  Not sure if this 
is a bug for some package or a data issue.  The data values do not look 
overly strange.



# data
sparse =
structure(list(x = c(740381.862, 740456.052, 740503.958, 740551.752,
740559.502, 740502.995, 740446, 740389.229, 740371.693, 740428.25,
740484.918, 740541.356, 740549.277, 740474.724, 740418.118, 740370.187,
740354.321, 740410.53, 740467.451, 740523.772, 740522.433, 740474.797,
740400.293, 740343.175, 740336.067, 740392.917, 740449.622, 740506.162,
740495.664, 740448.693, 740382.062, 740325.464, 740318.174, 740430.337,
740488.37, 740477.578, 740429.695, 740373.133, 740325.408, 740631.842,
740688.362, 740744.857, 740726.149, 740695.778, 740621.663, 740613.553,
740670.205, 740726.566, 740733.965, 740660.272, 740620.315, 740594.82,
740651.714, 740708.217, 740690.056, 740659.603, 740575.902, 740576.796,
740558.179), y = c(181644.086, 181620.772, 181605.577, 181590.417,
181568.637, 181586.847, 181604.615, 181622.136, 181565.531, 181547.708,
181530.169, 181512.439, 181490.328, 181513.956, 181531.875, 181547.048,
181508.946, 181491.148, 181473.394, 181455.233, 181436.726, 181452.342,
181475.522, 181492.661, 181451.96, 181434.265, 181416.566, 181398.729,
181382.764, 181397.808, 181418.748, 181436.677, 181395.477, 181360.409,
181342.547, 181327.09, 181341.971, 181359.55, 181374.453, 181546.959,
181528.576, 181510.58, 181497.501, 181507.159, 181530.759, 181490.008,
181472.209, 181453.968, 181432.87, 181456.085, 181468.588, 181433.854,
181415.758, 181397.497, 181384.359, 181393.795, 181420.579, 181376.982,
181363.899), depth_cm = c(-8, -8, -8, -8, -8, -8, -8, -8, -8,
-8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8,
-8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8,
-8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8,
-8, -8), theta_percent = c(23.63, 19.68, 23.81, 22.01, 23.98,
12.8, 14.92, 20.49, 22.59, 24.32, 20.24, 23.03, 12.97, 19.09,
39.2, 12.09, 24.52, 25.57, 25.5, 19.76, 19.17, 21.98, 7.5, 22.75,
17.85, 17.75, 17.95, 26.93, 18.84, 22.95, 23.71, 25.03, 40.69,
9.7, 24.66, 17.43, 16.3, 24.13, 19.98, 23.35, 12.16, 17.24, 14.29,
34.42, 21.84, 25.63, 20.51, 25.87, 24.44, 22.35, 8.57, 21.43,
25.63, 21.56, 21.49, 17.66, 25.61, 24.11, 28.31)), .Names = c("x",
"y", "depth_cm", "theta_percent"), row.names = c("1", "5", "10",
"15", "20", "25", "30", "35", "40", "45", "50", "55", "60", "65",
"70", "75", "80", "85", "90", "95", "100", "105", "109", "114",
"119", "124", "129", "134", "139", "144", "149", "154", "159",
"164", "169", "174", "179", "184", "188", "193", "198", "203",
"208", "213", "218", "223", "228", "233", "238", "243", "248",
"253", "258", "263", "268", "273", "278", "283", "288"), class = 
"data.frame")


# the broken fit for best search
require("automap")
coordinates(sparse) = c("x", "y", "depth_cm")
proj4string(sparse) = CRS("+init=epsg:32119")
v.fit <- autofitVariogram(theta_percent~1, sparse)


There were 50 or more warnings (use warnings() to see the first 50)
 > warnings()
Warning messages:
1: In getModel(initial_sill - initial_nugget, m, initial_range,  ... :
   An error has occured during variogram fitting. Used:
         nugget: 34.1432533936652
         model:  Exp
         psill:  13.2004974623731
         range:  53.1549477005646
         kappa:  NA
   as initial guess. This particular variogram fit is not taken into 
account.
Gstat error:
Error in if (direct[direct$id == id, "is.direct"] && any(model$psill <  :
   missing value where TRUE/FALSE needed

2: In fit.variogram(object, model, fit.sills = fit.sills,  ... :
   value out of range in 'bessel_k'
3 ...


# a quick visual of the data in the field
rescale = function(x, to=c(1,10)) (x - min(x)) * ((max(to) - 
min(to))/(max(x) - min(x)))
require("rgl")
sparse_df=as.data.frame(sparse)
spheres3d(sparse_df$x, sparse_df$y, sparse_df$depth_cm, 
radius=rescale(sparse_df$theta_percent))