Làm thế nào để đặt một bản đồ geom_sf sản xuất trên đỉnh của một raster ggmap sản xuất

Question

Tôi đã thử các mã sau:

library(ggplot2) 
library(ggmap) 
library(sf) 
nc <- st_read(system.file("shape/nc.shp", package = "sf")) 
str(nc) 

Classes ‘sf’ and 'data.frame': 100 obs. of 15 variables: 
$ AREA  : num 0.114 0.061 0.143 0.07 0.153 0.097 0.062 0.091 0.118 0.124 ... 
$ PERIMETER: num 1.44 1.23 1.63 2.97 2.21 ... 
$ CNTY_ : num 1825 1827 1828 1831 1832 ... 
$ CNTY_ID : num 1825 1827 1828 1831 1832 ... 
$ NAME  : Factor w/ 100 levels "Alamance","Alexander",..: 5 3 86 27 66 46 15 37 93 85 ... 
$ FIPS  : Factor w/ 100 levels "37001","37003",..: 5 3 86 27 66 46 15 37 93 85 ... 
$ FIPSNO : num 37009 37005 37171 37053 37131 ... 
$ CRESS_ID : int 5 3 86 27 66 46 15 37 93 85 ... 
$ BIR74 : num 1091 487 3188 508 1421 ... 
$ SID74 : num 1 0 5 1 9 7 0 0 4 1 ... 
$ NWBIR74 : num 10 10 208 123 1066 ... 
$ BIR79 : num 1364 542 3616 830 1606 ... 
$ SID79 : num 0 3 6 2 3 5 2 2 2 5 ... 
$ NWBIR79 : num 19 12 260 145 1197 ... 
$ geometry :sfc_MULTIPOLYGON of length 100; first list element: List of 1 
..$ :List of 1 
.. ..$ : num [1:27, 1:2] -81.5 -81.5 -81.6 -81.6 -81.7 ... 
    ..- attr(*, "class")= chr "XY" "MULTIPOLYGON" "sfg" 
    - attr(*, "sf_column")= chr "geometry" 
    - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ... 
    ..- attr(*, "names")= chr "AREA" "PERIMETER" "CNTY_" "CNTY_ID" ... 

map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google") 
a <- unlist(attr(map,"bb")[1, ]) 
bb <- st_bbox(nc) 
ggplot() + 
    annotation_raster(map, xmin = a[2], xmax = a[4], ymin = a[1], ymax = a[3]) + 
    xlim(c(bb[1], bb[3])) + ylim(c(bb[2], bb[4])) + 
    geom_sf(data = nc, aes(fill = AREA)) 

Hai lớp không chồng chéo đúng cách; Tôi đã thử thay đổi chiếu với coord_sf() nhưng tôi đã không thành công.

bất kỳ đề xuất nào? cảm ơn

Answer 1

Bạn có thể sử dụng phương pháp vẽ đồ thị từ gói sf để thực hiện việc này. Bạn sẽ cần phải xác định hệ quy chiếu tọa độ, mà chúng tôi sẽ cần phải giả định (và có vẻ như đúng như vậy) là 3857.

library(ggplot2) 
library(ggmap) 
library(sf) 

nc_shp <- st_read(system.file("shape/nc.shp", package = "sf")) 
nc_map <- get_map("north carolina", maptype = "satellite", zoom = 6, source = "google") 

st_crs(nc_map) 
# Coordinate Reference System: NA 

# assume the coordinate refence system is 3857 
plot(st_transform(nc_shp, crs = 3857)[1], bgMap = nc_map)