5

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

enter image description here

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. enter image description here

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.

Community
  • 1
  • 1
Yang Yang
  • 858
  • 3
  • 26
  • 49
  • see answer below. Also note that one of the problems you had was that latitude and longitude were inverted in the `interp_data` object. – lbusett Jan 14 '17 at 15:12
  • 1
    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. – Axeman Jan 17 '17 at 12:43
  • @Axeman Thanks a lot. – Yang Yang Jan 17 '17 at 20:58

1 Answers1

7

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:

![enter image description here

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: enter image description here

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 !

lbusett
  • 5,801
  • 2
  • 24
  • 47
  • 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. – Yang Yang Jan 16 '17 at 18:12
  • 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. – lbusett Jan 16 '17 at 23:10
  • 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. – lbusett Jan 16 '17 at 23:21
  • Thank you so much for your help. – Yang Yang Jan 17 '17 at 01:15
  • Hi. I just revised the answer to solve the "missing data" problem. – lbusett Jan 17 '17 at 14:09
  • Thanks a lot. You really save me! `spsample ` is a critical step in this problem. The points are all yours. – Yang Yang Jan 19 '17 at 09:32
  • Glad it helped. However, be careful with spsample, since "spatial extrapolation" far away from measurement data points can become rapidly unreliable. – lbusett Jan 19 '17 at 17:28
  • 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. – Yang Yang Apr 19 '17 at 17:42