0

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
Cœur
  • 37,241
  • 25
  • 195
  • 267
R.B
  • 1
  • 1

1 Answers1

1

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.

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
  • Thank you for your answer. I 've tried with utm coordinates. Now i do not get the error , but the final plot is not interpolated the whole map appears in one single color . This is my new code . Do yoy have any idea of the cause of this problem? Thank you very much – R.B May 01 '16 at 17:38
  • My first guess would be that the observations are outside the area where you want to predict. If the observations are far away the kriging prediction converges to a mean value, which is what you observe. Check the coordinates, maybe you did something wrong with the coordinate transform, for example by setting the wrong starting coordinate system before transforming to utm. Also note that utm works with zones, please also use the correct zone for your study area. – Paul Hiemstra May 01 '16 at 18:02
  • I added the zones again , I swear I have already written them maybe with all of this tryal and error I erase them , my mistake . By the way I used this line: result1<-automapPlot(kriging_result$krige_output,"var1.pred",sp.layout = list("sp.points", data1.utm));result1 to see where my points are situated and all of them appear inside the map . So I do not know where the error can be. And once again thank you – R.B May 01 '16 at 18:36