1

I have a problem with autoKrig function and I try to make a reproducible example here:

library(automap)
library(raster)
library(dismo)

bio <- getData("worldclim", var="bio", res=10)

bio1 <- raster(bio, layer=1)
bio12 <- raster(bio, layer=12)
predictors <- stack(bio1, bio12)

bg <- randomPoints(bio1, 50)

data <- extract(predictors, bg)
data <- cbind(bg,data)
data <- data.frame(data)

coordinates(data)=~x+y
proj4string(data) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
kg <- autoKrige(bio1~bio12, data, new_data=predictors)

This will result in:

Error in autoKrige(bio1 ~ bio12, data, new_data = predictors) : 
  Either input_data or new_data is in LongLat, please reproject.
   input_data:  +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
   new_data:    +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

I get the same error with my original data. I appreciate any help.

Geo-sp
  • 1,704
  • 3
  • 19
  • 42

1 Answers1

3

If you read the help file, it will tell you why it is throwing that error.

autoKrige performs some checks on the coordinate systems of input_data and new_data. If one or both is NA, it is assigned the projection of the other. If they have different projections, an error is raised. If one or both has a non-projected system (i.e. latitude-longitude), an error is raised. This error is raised because 'gstat does use spherical distances when data are in geographical coordinates, however the usual variogram models are typically not non-negative definite on the sphere, and no appropriate models are available' (Edzer Pebesma on r-sig-geo).

It looks like you need to project your data before calling autoKrige.

ialm
  • 8,510
  • 4
  • 36
  • 48
  • Thanks, but as the error shows both my input and output have the same projection. – Geo-sp Jul 19 '13 at 19:02
  • @N16 your data is in the `longlat` projection, which is what the help file details as a "non-projected system (i.e. latitude-longitude)". – ialm Jul 19 '13 at 19:06
  • So in this case what do you recommend? My data sets are global. – Geo-sp Jul 19 '13 at 19:12
  • Sorry, I only have experience with using kriging to solve very elementary problems. You can try projecting your data before calling `autoKrige`. For example, you can try the Robinson projection by calling `projectRaster(bio1, crs=CRS("+proj=robin"))`, though I suspect it may take some time. – ialm Jul 19 '13 at 19:27
  • @Ialm is right, you need to project your data from latlon to a projected system. – Paul Hiemstra Jul 20 '13 at 09:46