Plotting map resulted from kriging in R

Question

I want to plot a map weighted by prices. I got a data frame of latitudes, longitudes and prices:

> str(data)
'data.frame':   1253 obs. of  3 variables:
 $ lon  : num  11.2 11.1 11 11 11.3 ...
 $ lat  : num  49.6 49.4 49.3 49.4 49.3 ...
 $ price: num  1.4 1.37 1.4 1.36 1.4 ...

> head(data)
       lon      lat    price
1 11.21702 49.58926 1.395000
2 11.07774 49.37093 1.369000
3 11.04533 49.33787 1.397000
4 10.97310 49.43500 1.357333
5 11.34080 49.31650 1.399000
6 11.20708 49.36844 1.399000

My second step is to transform the data frame into a spatial data type:

coordinates(data) = ~lon+lat
proj4string(data) = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

And to plot it:

spplot(data, scales=list(draw=T), sp.layout=list("sp.points", data, pch=""), main="My Data")

(Should look like a small map of Germany)

At this point I got a plot with my points colored from black (lowest price) to yellow (highest price). But there is also a lot of white space in my plot. Now I want to fill the white space with assumptions, so it looks like a heatmap at the end.

I read some articles about "Kriging", but I'm not sure how to code this. I also read, that I need a raster. Therefore I tried the following lines according to this article How to make RASTER from irregular point data without interpolation:

s100 <- as.matrix(as.data.frame(data), ncol=3, byrow=TRUE)
e <- extent(s100[,1:2])
r <- raster(e, ncol=10, nrow=2)
x <- rasterize(s100[,1:2], r, fun='mean')

Answer 1

In order to interpolate prices with kriging you first need to convert your geographic coordinates to projected coordinates. Assuming you have them, below there is a reproducible example, showing a way to accomplish such task.

library(sp)

#Spatial data containing variables which can be interpolated.
#We will use the zinc column; as an equivalent for 'price'. 
data(meuse)  #data.frame with geometries (projected coordinates)
coordinates(meuse) = ~x+y #Convert data.frame object in SpatialPointsDataFrame object
head(meuse[,c(4)],3) #Show first three lines, with projected coordinates and zinc values.

     coordinates zinc
(181072, 333611) 1022
(181025, 333558) 1141
(181165, 333537)  640

#This will be the prediction grid, i.e., where the interpolated values will be printed.
data(meuse.grid)
gridded(meuse.grid) = ~x+y #Convert data.frame object to SpatialPixelsDataFrame object.
plot(meuse.grid) #this is the 'raster'(grid) you said, you'd need.

Now, we will use the function autoKrige to perform an ordinary kriging using zinc values from the meuse dataset, and print the interpolation result to the prediction grid ('meuse.grid').

library(automap)

#Perform ordinary kriging and store results inside object of type "autoKrige" "list" 
kriging_result = autoKrige(zinc~1, meuse, meuse.grid)
plot(kriging_result)

Here you go:

But now, let's suppose you don't have the prediction grid already available as we had with the 'meuse.grid' data. So, here is one example showing how to build the prediction grid:

data(meuse)

#Setting the  prediction grid properties
min_x = min(meuse$x) #minimun x coordinate
min_y = min(meuse$y) #minimun y coordinate
x_length = max(meuse$x - min_x) #easting amplitude
y_length = max(meuse$y - min_y) #northing amplitude
cellsize = 100 #pixel size
ncol = round(x_length/cellsize,0) #number of columns in grid
nrow = round(y_length/cellsize,0) #number of rows in grid

coordinates(meuse) = ~x+y
proj4string(meuse) = CRS("+proj=utm +ellps=WGS84 +datum=WGS84") #assign CRS with projected coordinates.

grid = GridTopology(cellcentre.offset=c(min_x,min_y),cellsize=c(cellsize,cellsize),cells.dim=c(ncol,nrow))

#Convert GridTopolgy object to SpatialPixelsDataFrame object.
grid = SpatialPixelsDataFrame(grid,
                              data=data.frame(id=1:prod(ncol,nrow)),
                              proj4string=CRS(proj4string(meuse)))

plot(grid)

And then, repeat the steps of the ordinary kriging with the autoKrige function. In this case, the interpolation will be performed in the whole graticulate, but the correct thing to do would be to perform the ordinary kriging only within the extents where the data is available. So, the below code shows an example (not reproducible, though) about using a shapefile to clip the extents of the prediction grid.

#Read shapefile and convert it to SpatialPolygons object.
#In the case of OP, the shapefile should be the administrative boundary of Germany.
shp = readOGR(dsn = "...", layer = "shp")
shp = shp@polygons
shp = SpatialPolygons(shp, proj4string=CRS("+proj=utm +ellps=WGS84 +datum=WGS84")) #make sure the shapefile has the same CRS from the data, and from the prediction grid.

#Clip the prediction grid with the shapefile.
grid = grid[!is.na(over(grid, shp)),]

If I had the shapefile with the boundaries of the 'meuse' dataset, this command would produce a grid equal to the command plot(meuse.grid). This post, illustrates how the clipping works.