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

Tue, 04 May 2010 02:06:27 -0700 

Hi Mark,
Thanks for the reproducible example. The problem is that when I look at the sample variogram, the semivariance values start high and end low. This is best illustrated by: 

plot(variogram(theta_percent~1, sparse))


You see that there are outliers in the data that cause high 
semi-variance at a short distance. I would say that is 'strange' ;). You 
can try and identify which point causes this by:


plot(variogram(theta_percent~1, sparse, cloud = TRUE), identify = TRUE)
# Click on the plot to identify the point pairs


There is not really one value I think that causes this. It might be also 
attributable to the fact that your dataset is somewhat sparse in the 
distance range from 20-40 m.



As a quick fix you can restrict the model selection to "Sph", this 
works. And I would not be enthusiastic about using Ste. This is because 
the main difference between different kappa values if the behavior at 
short distances, but you don't have a lot of data on the short distance 
to fit this value on in a meaningful way.



At this stage I don't see a (relatively quick) fix that could solve this 
problem in an automatic way and on a more fundamental level. Do you have 
any suggestions? From an implementation point of view I can let automap 
discard any model that has negative values in it, this would ensure that 
the user gets at least the Sph model back.


cheers,
Paul

Mark Connolly wrote:

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


_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



--
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo






















					Reply via email to