Как построить интерполяцию данных на проецируемой карте, используя ggplot2 в R

Question

Я хочу нанести некоторые интерполяционные данные на проецируемую карту, используя ggplot2, и я работаю над этой проблемой в течение нескольких недель. Надеюсь, что кто-то может мне помочь, большое спасибо. Шейп-файл и данные можно найти по адресу https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0 и https://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=0

Во-первых, шейп-файл изначально использует проекцию lon-lat, и мне нужно преобразовать его в проекцию Равной площади Альберса (aea).

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

Как мы видим, мы можем правильно построить страну. Затем я хочу интерполировать данные, используя метод Kriging, код взят из 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)

Тогда мы можем увидеть интерполированную карту данных.

Наконец, я хочу объединить эти две цифры, а затем я получаю следующую ошибку: 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 

Мой вопрос: как я могу отобразить интерполяционные данные на карте, а затем добавить линии сетки (долгота и широта). Кроме того, мне интересно, как обрезать интерполированную карту данных, чтобы она соответствовала всей карте Канады. Спасибо за помощь.

Answer 1

После копания немного больше, я думаю, вы можете захотеть это:

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

который дает:

Единственная проблема состоит в том, что у вас все еще есть "пропущенные" интерполированные области на итоговой карте (например, в западной части). Это связано с тем, что, начиная с autokrige Помогите:

new_data: объект sp, содержащий местоположения предсказания. new_data может быть набором точек, сеткой или многоугольником. Не должно содержать NA. Если этот объект не предоставлен, вычисляется значение по умолчанию. Это делается путем взятия выпуклой оболочки input_data и размещения около 5000 ячеек сетки в этой выпуклой оболочке

Поэтому, если вы не предоставите допустимые новые данные в качестве аргумента, интерполированная область будет ограничена выпуклой оболочкой точек входного набора данных (= без экстраполяции). Это может быть решено с помощью spsample вsp пакет:

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

который дает:

Обратите внимание, что маленькие "дыры", которые у вас все еще есть вблизи границ многоугольников, можно удалить, увеличив число точек интерполяции в вызове spsample (Так как это медленная операция, я просил только 4000, здесь)

Более простой и быстрой альтернативой может быть использование пакета mapview

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

(вы можете использовать менее детализированные шейп-файлы границ полигонов, так как это медленно)

HTH!