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

Tue May 4 11:49:16 CEST 2010

I just uploaded a new version of automap (1.07) that fixes this problem. 
It deletes fitted variogram models with a negative sill/range/nugget.

cheers,
Paul

Paul Hiemstra wrote:
> 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 at 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