title: Overlay raster layer on map in ggplot2 in R?
body: <p>I am trying to overlay a raster layer onto a map in ggplot. The raster layer contains likelihood surfaces for each time point from a satellite tag. I also want to set cumulative probabilities(95%, 75%, 50%) on the raster layer.</p>
  
  <p>I have figured out how to show the raster layer on the ggplot map, but the coordinates are not aligned with one another. I tried making each have the same projection but it does not seem to be working... I want them both to fit the boundaries of my model (xmin = 149, xmax = 154, ymin = -14, ymax = -8.75</p>
  
  <p>Attached is my r code and the figure result:</p>
  
  <pre><code>#load data
  ncname &lt;- "152724-13-GPE3"
  ncfname &lt;- paste(ncname, ".nc", sep = "")
  ncin &lt;- nc_open(ncfname)
  
  StackedObject&lt;-stack("152724-13-GPE3.nc", varname = "monthly_residency_distributions")
  MergedObject&lt;-overlay(StackedObject,fun=mean ) 
  MergedObject[is.na(MergedObject)]&lt;-0 
  Boundaries&lt;-extent(c(149, 154, -14, -8.75)) 
  ExtendedObject&lt;-extend(MergedObject, Boundaries) 
  Raster.big&lt;-raster(ncol=1200,nrow=900,ext=Boundaries) 
  Raster.HR&lt;-resample(x=ExtendedObject, y=Raster.big, method="bilinear") 
  Raster.HR@data@values&lt;- Raster.HR@data@values/sum(Raster.HR@data@values) 
  RasterVals&lt;-sort(Raster.HR@data@values)
  Raster.breaks &lt;- c(RasterVals[max(which(cumsum(RasterVals)&lt;= 0.05 ))], RasterVals[max(which(cumsum(RasterVals)&lt;= 0.25 ))], RasterVals[max(which(cumsum(RasterVals)&lt;= 0.50 ))], 1)
  Raster.cols&lt;-colorRampPalette(c("yellow","orange","red"))
  RasterCols&lt;- c(Raster.cols(3))
  
  #Create Map
  shape2 &lt;- readOGR(dsn = "/Users/shannonmurphy/Desktop/PNG_adm/PNG_adm1.shp", layer = "PNG_adm1")
  map&lt;- crop(shape2, extent(149, 154, -14, -8.75))
  projection(map)&lt;- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  
  p &lt;- ggplot() + geom_polygon(data = map, aes(x = long, y = lat, group = group), color = "black", size = 0.25) + coord_map()
  
  projection(Raster.HR)&lt;- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  
  #plot raster and ggplot 
  
  par(mfrow=c(1,1))
  plot(p)
  par(mfrow=c(1,1), new = TRUE)
  plot(Raster.HR, col=RasterCols, breaks=Raster.breaks, legend = NULL, bbox(map)) 
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/fvQBJ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/fvQBJ.png" alt="enter image description here"></a></p>
  
  <p>Please let me know if there is another package/line of code I should be using to do this! Appreciate any help</p>