0

I started using the automap package for automated variogram modelling two days ago and played around with it (doing it manually was tedious and doesn't come with too many advantages...although I find it strange that automap seems really to love the Matérn function for modelling the variogram). What came to my mind was the presence of negative values in the dataset while those were not shown on the map that is produced by this package.

Here is a figure I just glued together: Nr. 1 is the output after using autoFitVariogram and autoKrige. It appears to have no negative values. As I didn't like the graphical output, I performed ordinary kriging with gstat that uses the variogram from autoFitVariogram (from automap) (Nr. 2) and it clearly shows areas of negative precipitation. I exported the map produced from autoKrige to QGIS and double checked and there were also lots of negative values (Nr. 3; blackened areas).

enter image description here

So my questions:

1.) The kriging prediction map produced by autoKrige apparently produces a map that follows a scheme such as "if var.pred < 0 then 0". I can see this is usefull for producing such a map, but how does this make sense with the actual data? I want to use my dataset for further analyses and negative rainfall values appear a bit...non-realistic. Would it be mathematically sound to just apply this "if var.pred < 0 then 0" scheme also to the data if this would make it physically plausible?

2.) You may determine variogram parameters additionally with Generalized Least Squares (GLS) as an own argument in the automap package. However, I want to use a REML fit from gstat (fit.variogram.reml). I tried substituting variogram$var_model with my REML estimates (which also produced my modified variogram), but when performing kriging with autoKrige the original variogram model is used again. Is there a way to fit the model by other means than those provided by the package? (there is an argument 'data_variogram' in the autoKrige manual, but this seems to do something else).

Thanks!

GeoEki
  • 437
  • 1
  • 7
  • 20
  • I am quite sure automap does not restrict predictions to above zero. This is probably a by product of the variogram you are using. In addition, you cannot use fit.variogram.reml right now. – Paul Hiemstra Jan 22 '16 at 13:11

1 Answers1

0
  1. I don't think this is the case; you seem to conclude this from looking at a graph, but please verify the numbers produced. Graph 1 does not appear to me that it does not contain negative values.
  2. autoKrige uses fit.variogram, and cannot use fit.variogram.reml; that is how it's coded, but feel free to modify it.
Edzer Pebesma
  • 3,814
  • 16
  • 26
  • 1. Graph 1 and 3 are from the identical output (Nr. 1 from `autoKrige` output, Nr. 3 the same map exported as an ascii file and imported in a GIS where I colored negative values in black, whereas graph 2 (not based on `autoKrige` but on `autofitVariogram`, so the same data) exposes areas with negative values. If graph 1 also contains negative values, the color scheme or the signature is really misleading then. I tried substituting all negative values for graph 2 with 0 (`ifelse(prediction_spdf$var1.pred<0,0,prediction_spdf$var1.pred)`) and it looked identical to graph 1. – GeoEki Jan 22 '16 at 11:40
  • 2. Such a pitty. Unfortunately I’m not capable of writing my own packages (yet), but Paul may take this as a possible future implementation :-) – GeoEki Jan 22 '16 at 11:40