Unable to properly use autokrige from automap package (R cannot read the prediction locations well)

1k Views Asked by At

I'm trying to use R to perform a map of interpolated frequencies of data collected from Iberian Peninsula. (something like this https://gis.stackexchange.com/questions/147660/strange-spatial-interpolation-results-from-ordinary-kriging )

My problem is that the plot is not showing the interpolated data, due to some kind of error in the atributte new_data of the autokrige function.

https://cran.r-project.org/web/packages/automap/automap.pdf new_data: A sp object containing the prediction locations. new_data can be a points set, a grid or a polygon. Must not contain NA’s. If this object is not provided a default is calculated. This is done by taking the convex hull of input_data and placing around 5000 gridcells in that convex hull.

I think the problem is that R is not reading well the map transformed to poligons because if I avoid this new_data attribute I get a propper plot of the krigging values. but I do not obtain a good shape of the iberian peninsula. Can someone help me please? I would be truly grateful

here you can see my data: http://pastebin.com/QHjn4qjP

Actual code: now since I transform my data coordinates to UTM projection i do not get error messages but the last plot is not interpolated, the whole map appear of one single color :(

 setwd("C:/Users/Ruth/Dropbox/R/pruebas")

#Libraries

library(maps)
library(mapdata)
library(automap)
library(sp)
library(maptools)
library(raster)
library(rgdal)

####################MAPA#############

#obtain the grid of the desired country
map<-getData('GADM', country='ESP', level=0)

#convert the grid to projected data
mapa.utm<-spTransform(mapa3,CRSobj =CRS(" +proj=utm +zone=29 +zone=30 +zone=31  +ellps=WGS84"))

###############################Datos#######################

#submit the data
data1<-read.table("FRECUENCIASH.txt",header=T)
head(data1)
attach(data1)
#convert longlat coordinates to UTM
coordinates(data1)<-c("X","Y")
proj4string(data1) = CRS("+proj=longlat +datum=WGS84") 
data1.utm=spTransform(data1, CRS("+proj=utm +zone=29 +zone=30 +zone=31 +ellps=WGS84 "))

######################Kriging interpolation #####################


#Performs an automatic interpolation

kriging_result<-autoKrige(Z~1,data1.utm,mapa.utm,data_variogram = data1.utm) 

#Plot the kriging

result1<-automapPlot(kriging_result$krige_output,"var1.pred",sp.layout = list("sp.points", data1.utm));result1

result2<-plot(kriging_result);result2
1

There are 1 best solutions below

3
On

The error you get is related to the fact you are using unprojected data as input for automap. Automap can only deal with projected data. Googling for map projections should give you some background information. To project your data to a suitable projection, you can use spTransform from the sp package.

The fact that it works without newdata is because in that object the projection is not set, so automap cannot warn you. However, the results of automap with latlon input data is not reliable.