0

I am trying to do kriging with my data. But there is a problem here. I want to use krige command in gstat package. But I do not have "Spatial object with prediction/simulation locations".

How I can make SpatialPolygons?

I have a shapefile but the problem is that I want to do kriging outside of my shapefile.

  • It sounds like you are asking how to load a shapefile into R as a SpatialPolygons or SpatialPolygonsDataFrame. In which case the answer is to use the maptools or rgdal package. But this doesn't really have anything to do with kriging (or any particular use you will make of the SpatialPolygons), so I want to make sure I understand the question – Lee Hachadoorian Apr 27 '16 at 12:45
  • No. I want to create a grid to visualise the output of kriging results. – mohammad gerami Apr 27 '16 at 12:47
  • The question in this comment ("create a grid to visualise the output of kriging results") is not the same as the question you have posted ("Creating SpatialPolygons in R"). In addition, the gstat package has a lot of options in it. I suggest you edit your question and include a code sample so we can see how far you have got and what you want to accomplish next. – Lee Hachadoorian Apr 27 '16 at 13:11
5

You are confusing terms and thus, confusing us. The expected input for kriging prediction in the gstat krige function is a systematic array of points and not polygons. It would also be nice if you provided a reproducible code example of what you have tried.

You can use the extent of an sp object to create an array of points for the kriging prediction using the extent, raster and rasterToPoints functions in the raster package.

First, we add our libraries and example data (meuse SpatialPointsDataFrame data from sp library).

library(gstat)                                         
library(sp)
library(raster)                                            

data(meuse)
coordinates(meuse) <- ~x + y
proj4string(meuse) <- CRS("+init=epsg:28992")

Here we create an extent polygon using a bounding extent of any relevant sp or raster class object. You can then use this extent polygon to create a dummy raster object that is coerced into a SpatialPoints object using rasterToPoints. The resulting sp SpatialPoints object can act as the kriging prediction grid.

ext <- as(extent(meuse), "SpatialPolygons")
r <- rasterToPoints(raster(ext, resolution = 30), spatial = TRUE)
proj4string(r) <- proj4string(meuse)  

Now we can specify a kriging model using our point (meuse) and prediction grid (r) objects.

# krige model log(copper):                                              
v1 <- variogram(log(copper) ~ 1, meuse)                  
x1 <- fit.variogram(v1, vgm(1, "Sph", 800, 1))           
G1 <- krige(copper ~ 1, meuse, r, x1, nmax = 30)
gridded(G1) <- TRUE                                      
G1@data = as.data.frame(G1@data[,-2])
plot(G1)

For some reason the raster function was not always honoring the crs argument. I fixed it by assigning the projection, using proj4string, after the points were created.

  • Thank you very much. That was the issue that I wanted to solve it. – mohammad gerami Apr 28 '16 at 1:50
  • Could you please describe the "proj4string" command? what does it for? – mohammad gerami Apr 28 '16 at 1:56
  • Another problem: When I wanted to run <G1 <- krige(zinc ~ 1, meuse, r, x1, nmax = 30)> this error appeared: Error in predict.gstat(g, newdata = newdata, block = block, nsim = nsim, : var1 : data item in gstat object and newdata have different coordinate reference systems > – mohammad gerami Apr 28 '16 at 2:02
  • The proj4string defines the projection. If this response answered your question pleas mark it as answered so everybody knows that it has been resolved. – Jeffrey Evans Apr 28 '16 at 2:06
  • I resolved the error you were receiving. For some reason the projection was not consistently assigning within the raster function. Assigning the projection so they will match, is a critical step. This is done using the proj4string as illustrated in the code example. The "meuse" SpatialPointsDataFrame object is acting in place of whatever data you are using. – Jeffrey Evans Apr 28 '16 at 2:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.