Plotting static base map underneath a sf object

You can also use the package ggspatial which provides the "map tile" annotation layer.

ggplot(nc_centers) +
    annotation_map_tile(zoom = 7) +
    geom_sf(aes(color = SID79, size = BIR74),
            show.legend = "point", inherit.aes = FALSE) +
    coord_sf(datum = NA) +
    theme_minimal()

I've been working on a package recently which someone may find useful. The ggmap answer now require an API key for use with google maps which adds an extra bit of faff.

basemapR also allows you a bit more flexibility in setting the extents for the basemap by using a bounding box.

#devtools::install_github('Chrisjb/basemapR')
library(basemapR)
library(sf)
library(ggplot2)


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

nc_centers <- st_centroid(nc)

# create bbox from our nc layer and expand it to include more area above/below
bbox <- expand_bbox(st_bbox(nc_centers), X = 0, Y = 150000)

ggplot() +
  base_map(bbox, increase_zoom = 2, basemap = 'google-terrain') +
  geom_sf(data = nc_centers, 
          aes(color = SID79, size = BIR74),
          show.legend = "point", inherit.aes = FALSE) +
  coord_sf(datum = NA,
           xlim = c(bbox['xmin'], bbox['xmax']),
           ylim = c(bbox['ymin'], bbox['ymax'])) +
  theme_minimal() +
  labs(caption = 'map data \uA9 2020 Google')

Remember to cite google maps in the caption or elsewhere in the map

You might consider reprojecting your data but the following code seems to work for me.
See here for an explanation about why you need inherit.aes = FALSE and see here for an alternative solution with base plots.

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2
# devtools::install_github("r-lib/rlang")
library(ggplot2)
library(ggmap) 

nc <- st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source `/home/gilles/R/x86_64-pc-linux-gnu-library/3.4/sf/shape/nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs
nc_map <- get_map(location = "North Carolina, NC", zoom = 7)
#> Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=North+Carolina,+NC&zoom=7&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
#> Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=North%20Carolina,%20NC&sensor=false
nc_centers <- st_centroid(nc)
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
#> of_largest_polygon): st_centroid does not give correct centroids for
#> longitude/latitude data

ggmap(nc_map) +
    geom_sf(data = nc_centers, 
            aes(color = SID79, size = BIR74),
            show.legend = "point", inherit.aes = FALSE) +
    coord_sf(datum = NA) +
    theme_minimal()
#> Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Created on 2018-04-03 by the reprex package (v0.2.0).