How to create geometry created by geom_sf on top of the generated ggmap raster

Question

I tried the following code:

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)) 

Two layers do not overlap properly; I tried changing the projection using coord_sf() , but I had no success.

any suggestion? thanks

Answer 1

To do this, you can use the graph plotting method from the sf package. You will need to specify the coordinate system that we will need to read (and this looks right) - 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)