Hot answers tagged

14

This error is commonly returned because you have duplicate locations. You can check this using the sp::zerodist function. To remove duplicate locations you call sp::zerodist within a bracket index. WeatherData <- WeatherData[-zerodist(WeatherData)[,1],]


7

Some has already been said by Spacedman in the comments. This warning may not pose a problem, if the variogram looks good. A good option might be to initialize some of the variogram parameters in vgm("Sph"). I usually take these values as default: vario <- variogram(Temperature ~1, WU_data_spatial) vario.fit <- fit.variogram(vario, vgm(psill=max(...


5

When the prediction error of a model based on log(Y) is Gaussian like here, you normally should use the correction of Laurent (1963) to estimate back to the original scale. In your case: MeanY = exp(a$var1.pred + 0.5 * (a$var1.var)) and: SdY = exp(2*a$var1.pred + a$var1.var)*(exp(a$var1.var)-1) You can verify on yourself with simulation of a Gaussian ...


5

After the study of https://rpubs.com/nabilabd/118172 and https://rpubs.com/nabilabd/134781 I assembled the following solution. library(sf) library(sp) library(gstat) library(tidyverse) Here are points I would like to play with. dt<-tibble::tribble( ~id, ~lon, ~lat, ~z, "a", 1.1, 1.3, 12.5, "b", 2.4, 1.7, 13.0, "c", 3.2, 1.4, ...


5

In sp, SpatialPoints*, SpatialPixels* and SpatialGrid* (with * omitted or replaced by DataFrame) do support more than 2 spatial dimensions, as OP has done, but SpatialPolygons* and SpatialLines* do not. With gstat you can do 3-D block kriging with 3-D blocks (using block = c(10,10,10)), but you cannot do this for non-rectangular blocks, as OP wants. It is ...


3

Once you've got a variogram fitted then kriging is trivially parallelizable. Kriging predictions are independent of each other. So, divide your prediction points (grid) into N sets, where N is your number of cores, and do the predictions for each of the point sets on a separate core. Merge the predictions afterwards. You can use any of the paralleling ...


3

You really need to do a bit more research on your methodology and read the documentation to understand the structure of sp S4 class objects and interaction of sp objects with relevant gstat functions. In the sp Vignette there is a detailed explanation of the difference between SpatialPolygons (only polygon topology) and SpatialPolygonDataFrame (polygons with ...


2

You should remove your NA in your data frame first using the function na.omit (see this: remove-rows-with-nas-in-data-frame) and then apply the variogram function. # Load libraries library('gstat') # Example data newdf <- data.frame("id" = 1:10, "yield_cleaned" = runif(n = 10, min = 0, max = 1), "x" = sample(x = ...


2

It depends on the fit.method chosen in the function fit.variogram. If you use: fit.method = 0, no fit fit.method = 1, weights are equal to Nj fit.method = 2, weights are equal to Nj / {γ(hj)}2 fit.method = 5, REML (restricted maximum likelyhood) fit.method = 6, no weights (OLS = Ordinary Least Squares) fit.method = 7, Nj / hj2 Please see this: Table 4.2: ...


2

A common problem is when you have duplicate locations. When I tried to add the following line of code before creating a spatial object from t1: ... t1 = t1[which(!duplicated(t1[1:2])), ] coordinates(t1) <- ~ x+y ... (which subsets rows from data frame t1 which have unique x and y values) the kriging process runs just fine and produces results. However, ...


2

This is a simple and reproducible example with meuse dataset of gstat package based on @Spacedman answer and using the parallel package: More info and help here: Parallelizing and clustering in R # Load libraries library('sp') library('gstat') library('parallel') # Load example data: meuse dataset data("meuse") # data data("meuse.grid") # 40m x 40m ...


2

A common misconception is that kriging estimates may be simply exponentiated to recover the field values. Sebastien Rochette's suggests a back-transformation for field values y following Laurent (1963): Because the prediction of log(y) is based on a Gaussian distribution, in many cases an additional correction factor is needed because the expected value of ...


2

I've done it like this: first make a grid to which you'll interpolate your met station values, library(gstat) library(ggplot2) library(rgdal) grd <- expand.grid(x=seq(from=bbox(met_stations)[1,1], to=bbox(met_stations)[1,2], by=0.05), y=seq(from=bbox(met_stations)[2,1], to=bbox(met_stations)[2,2], by=0.05)) coordinates(grd) <- c("...


2

It is a consequence of using a Gaussian variogram model with zero nugget; this easily leads to near-singular covariance matrices even if points don't (exactly) overlap. Solutions: choose a different variogram model (Matern, with strong smoothing?) or add a small nugget effect (very small may already help).


2

stAni: A spatio-temporal anisotropy; the number of space units equivalent to one time unit. and it is only used in spatio-temporal variograms that have a metric to them. By this it is meant that you can compute a distance between two points S at (x,y,t) and S' at (x',y',t'). For the "separable" variograms, the semi-variance is a product of a ...


2

I think there's a difference in the variogram between the model fitted here with gstat and the one in the linked question. In the linked question, the model is expressed in terms of the covariance as C: C(h) = σ² = 4, if |h| = 0 C(h) = (σ² - a)exp(-3|h|/r) = (4 - 1)exp(-3|h|/6), if |h| > 0 In R, C = function(h){3*exp(-3*h/6)} And the variogram from ...


2

Have you read the help for gstat::krige: newdata: data frame or Spatial object with prediction/simulation locations; should contain attribute columns with the independent variables (if present) and (if locations is a formula) the coordinates with names as defined in ‘locations’ You create a newdata data frame with the ...


1

The proj string for your newdata appears to be coordinate reference system 4326, which is in units degrees. If your data is, as you say, in units meters, then they are indeed not in identical coordinate reference systems. Try reprojecting your data to that string returned in the error first and see if that solves the problem: library(rgdal) ...


1

At this point dat isn't a spatial object, so when predict refers back to it it can't get any spatial data: idw_ff = gstat(formula = ff ~ 1, data = dat, nmax = neighbors, set = list(idp = beta)) make it spatial: coordinates(dat)=~dist+rel_alt and it works. idw_ff = gstat(formula = ff ~ 1, ...


1

From the help file: ?gstat::idw nmin = 0 nmax = Inf maxdist = Inf


1

Almost. Since the output from kriging is a distribution at every prediction point, you want to divide the sd by the prediction at each point: krig1$cov = 100 * sqrt( krig1$var1.var) / krig1$var1.pred Your code which looks like this: krig1$cov = 100 * sqrt( krig1$var1.var) / mean(krig1$var1.pred) is the variance multiplied by a constant everywhere.


1

You will find all you need in the excellent (and didactic) technical note from Rossiter (2012)*: Technical Note: Co-kriging with the gstat package of the R environment for statistical computing. Co-kriging will use different functions from those with univariate kriging (for example, ordinary kriging). The datasets (target and co-variables) should ...


1

This was reported as https://github.com/edzer/gstat/issues/33 and fixed at source. Before a new release becomes available, you need to install the package from source: library(devtools) install_github("edzer/gstat")


1

Well cross validation is a method to help test your models when you do not have enough data to separate into train and validation datasets. So by doing k-fold cross validation you can see how the performance of your model varies as the training/validation data partitions are changed. If you are happy with the CV performance, then just implement using the ...


1

Using an amalgamation of code and examples, I came up with something that looks correct to me: library(ggplot2) # start needed libraries library(gstat) library(sp) library(maptools) my_location<-c(lon= -84.94, lat = 45.72) basemap <- get_map(location = my_location, zoom = 14, maptype = "hybrid") plot (basemap) #load data field_data<-read.csv ("...


1

I found a very simple solution for obtaing the fitted gamma values within gstat package. It's the variogramLine() function. A simple code is attached here. # Empirical variogram ev = variogram("pH", data=data,....) fv = fit.variogram(q, "Sph", ....) fitted=variogramLine(v.fit, maxdist=max(q$dist), dist_vector=q$dist) fitted # see what are the values of ...


1

I couldn't answer your question using gstat package. However, you can also use geoR package to fit a variogram model to an empirical variogram and analyse fitted and residuals values. I give you a reproducible example below: # Load libraries ---------------------------------------------------------- library("sp") library("geoR") # Load data --------------...


1

Nine points is a quite small number, so I would use some arbitrary boundary instead of trying to build a complex algorithm that might "go wild". I suggest that you use r.grow.distance in grass to create a distance layer around your points, and to set a threshold that would constraint the size of your study area (for instance, the largest distance value ...


Only top voted, non community-wiki answers of a minimum length are eligible