How to plot interpolating data on a projected map using ggplot2 in R

Question

I want to plot some interpolating data on a projected map using ggplot2 and I have been working on this problem for a few weeks. Hope someone can help me, thanks a lot. The shapefile and data can be found at https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0 and https://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=0 .

First, the shapefile is originally using "lon-lat" projection and I need to convert it to Albers Equal Area (aea) projection.

library(automap)
library(ggplot2)
library(rgdal)
load("Mydata.rdata",.GlobalEnv)
canada2<-readOGR("gpr_000b11a_e.shp", layer="gpr_000b11a_e")
g <- spTransform(canada2, CRS("+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
Borders=ggplot() +geom_polygon(data=g,aes(x=long,y=lat,group=group),fill='white',color = "black")
Borders

As we can see, we can plot the country correctly. Then I want to interpolate the data using Kriging method, the code is taken from Smoothing out ggplot2 map.

coordinates(Mydata)<-~longitude+latitude
proj4string(Mydata)<-CRS("+proj=longlat +datum=NAD83")
sp_mydata<-spTransform(Mydata,CRS(proj4string(g)))
Krig=autoKrige(APPT~1,sp_mydata)
interp_data = as.data.frame(Krig$krige_output)
colnames(interp_data) = c("latitude","longitude","APPT_pred","APPT_var","APPT_stdev")
interp_data=interp_data[,1:3]
ggplot(data=interp_data, aes(x=longitude, y=latitude)) + geom_tile(aes(fill=APPT_pred),color=NA)

Then we can see the interpolating data map.

Finally I want to combine these two figures and then I get the following error: Error: Don't know how to add o to a plot

ggplot(data=interp_data, aes(x=longitude, y=latitude)) + geom_tile(aes(fill=APPT_pred),color=NA)+Borders 

My question is: how can I plot the interpolating data on the map and then add grid lines (longitude and latitude). Also, I wonder how to clip the interpolating data map to fit the whole Canada map. Thanks for the help.

Tonedeaf

see answer below. Also note that one of the problems you had was that latitude and longitude were inverted in the interp_data object. – Raimondo

The be clear, the error is raised because you can't + two objects created with ggplot(), so no ggplot() + geom_x() + ggplot() + geom_y(). You can only add new layers. – Alvarez

@Alvarez Thanks a lot. – Tonedeaf

Answer 1

After digging a bit more, I guess you may want this:

Krig = autoKrige(APPT~1,sp_mydata)$krige_output
Krig = Krig[!is.na(over(Krig,as(g,"SpatialPolygons"))),]  # take only the points falling in poolygons
Krig_df = as.data.frame(Krig)
names(Krig_df) = c("APPT_pred","APPT_var","APPT_stdev","longitude","latitude")
g_fort = fortify(g)
Borders = ggplot() +
  geom_raster(data=Krig_df, aes(x=longitude, y=latitude,fill=APPT_pred))+
  geom_polygon(data=g_fort,aes(x=long,y=lat,group=group),
               fill='transparent',color = "black")+
  theme_bw()
Borders

which gives:

Only problem is that you still have "missing" interpolated areas in the resulting map (e.g., on the western part). This is due to the fact that, as from autokrige help:

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

Therefore, if you do not provide a feasible newdata as argument, the interpolated area is limited by the convex hull of the points of the input dataset (= no extrapolation). This can be solved using spsample insp package:

library(sp)
ptsreg <- spsample(g, 4000, type = "regular")   # Define the ouput grid - 4000 points in polygons extent
Krig = autoKrige(APPT~1,sp_mydata, new_data = ptsreg)$krige_output
Krig = Krig[!is.na(over(Krig,as(g,"SpatialPolygons"))),]  # take only the points falling in poolygons
Krig_df = as.data.frame(Krig)
names(Krig_df) = c("longitude","latitude", "APPT_pred","APPT_var","APPT_stdev")
g_fort = fortify(g)
Borders = ggplot() +
  geom_raster(data=Krig_df, aes(x=longitude, y=latitude,fill=APPT_pred))+
  geom_polygon(data=g_fort,aes(x=long,y=lat,group=group),
               fill='transparent',color = "black")+
  theme_bw()
Borders

which gives:

Notice that the small "holes" that you still have near polygon boundaries can be removed by increasing the number of interpolation points in the call to spsample (Since it is a slow operation I only asked for 4000, here)

A simpler quick alternative could be to use package mapview

library(mapview)
m1 <- mapview(Krig)
m2 <- mapview(g)
m2+m1

(you may want to use a less detailed polygon boundaries shapefiles, since this is slow)

HTH !

Raimondo

Thanks, the map looks nice but there is still a question. Some parts of the country are not covered by the "APPT" value, such as the western part. Is there any solution? Thanks a lot. – Tonedeaf

from autokrige help: 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. – Raimondo

Therefore, if you do not provide a feasible newdata as argument, the interpolated area is limited by the convex hull of the points of the input dataset (= no extrapolation). You could create the necessaary dataset using somthing like SpatialPixels or SpatialGrid. I'll try to amend the answer ASAP. – Raimondo

Thank you so much for your help. – Tonedeaf

Hi. I just revised the answer to solve the "missing data" problem. – Raimondo

Thanks a lot. You really save me! spsample is a critical step in this problem. The points are all yours. – Tonedeaf

Glad it helped. However, be careful with spsample, since "spatial extrapolation" far away from measurement data points can become rapidly unreliable. – Raimondo

Hi, sorry to trouble you again. I wonder if I want to use a color gradient such as scale_fill_gradientn(colours = heat.colors(100), space = "Lab") to show the gradient of the data,how I can do this? I have tried many times, but failed. Thank you for your time. – Tonedeaf