How to plot the shapefile via ggplot2?

Question

I have a raster data and polygons of parks and I want to overlap it on the raster. When I add the polygon it shows here but on ggplot how I add polygons (polygons of parks is like round shapes)on my raster data through ggplot2,. My code is attached below.

   r <- raster(t((volcano[,ncol(volcano):1] - 94) * 4.95))
   pg <- readOGR("E:/park/1aa.shp") # loadshapfile 
   plot(r)
   plot(pg, add= TRUE,) # it appears here like first picture (left).

But how can I add this polygons o parks in my ggplot 2. My code of ggplot 2 is attached below.

  centile90 <- quantile(r, 0.90)
  df <- as.data.frame(as(r, "SpatialPixelsDataFrame"))
  colnames(df) <- c("value", "x", "y")
  library(ggplot2)

   mybreaks <- seq(0, 500, 50)

   ggplot(df, aes(x, y, z = value)) +
  geom_contour_filled(breaks = mybreaks) +
  geom_contour(breaks = centile90, colour = "pink",
          size = 0.5) +
   scale_fill_manual(values = hcl.colors(length(mybreaks) - 3, "Zissou1", rev = FALSE)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme()

Help is needed how to add ** pg (polygon) ** in my ggplot2 code.

Update 1 Description of polygon data

Answer 1

As explained, it is much handy to work with sf than sp, on top of that sf is meant to superseed sp.

Find here a reproducible example. The first part is just for mocking your file "E:/park/1aa.shp". Since it was not provided it was not possible for me to use your real data, but let's just pretend it is the same dataset...:

library(raster)
library(sf)
library(sp)

r <- raster(t((volcano[,ncol(volcano):1] - 94) * 4.95))

# Let's mock your shapefile
poly <- st_as_sfc(st_bbox(st_as_sf(rasterToPolygons(r))))

# Sample 4 points
set.seed(3456)

sample <- st_sample(poly, 4)
sample <- st_buffer(sample, c(0.01, 0.02, 0.03))
sample <- st_sf(x=1:4, sample)
st_write(sample, "1aa.shp", append = FALSE)
# Mocked data

# Now let's start with your code -------
library(raster)
library(sf)

r <- raster(t((volcano[,ncol(volcano):1] - 94) * 4.95))

# Use sf!!
pg <- st_read("1aa.shp") # loadshapfile 
plot(r)
plot(st_geometry(pg), add= TRUE,) # it appears here like first picture (left).

Now work with geom_sf() on your pg object:

centile90 <- quantile(r, 0.90)
df <- as.data.frame(as(r, "SpatialPixelsDataFrame"))
colnames(df) <- c("value", "x", "y")

library(ggplot2)

mybreaks <- seq(0, 500, 50)

ggplot(df, aes(x, y, z = value)) +
  geom_contour_filled(breaks = mybreaks) +
  geom_contour(breaks = centile90, colour = "pink",
               size = 0.5) +
  # And here we have it
  geom_sf(data=pg, fill="black", inherit.aes = FALSE) +
   scale_fill_manual(values = hcl.colors(length(mybreaks)-1, "Zissou1", rev = FALSE)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme_classic() +
  theme()