multiple Map with subplots in R

Question

I have a series of rasters plots that I want on a single page. Every raster plot should be together with a subplot of specified region in the raster. I found some related methods that fail to work with rasters plot

This map is a demo made by Arcgis , I want to implement it with R script.

First of all ,I want to try a single plot with sub plot graph. Here is my sample code .

    require(raster)
    require(rgdal)
    op <- par(no.readonly = TRUE)
    oldwd<-getwd()
    setwd(tempdir())
    plotr<-raster(nrow=405,ncol=485,
                  xmn=-2638000,xmx=2212000,ymn=1872000,ymx=5922000,
                  crs='+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs')
    plotr<-setValues(plotr, (coordinates(plotr)[,2]))
    chnshp<-getData('GADM', country="China", level=1)
    chnshp<-spTransform(chnshp,CRS=CRS(projection(plotr)))
    plotr<-mask(plotr,chnshp)
    breaks=seq(min(minValue(plotr)),max(maxValue(plotr)),length.out=10)
    download.file('https://dl.dropboxusercontent.com/s/vzs6166v46de7v7/subplot_frame.zip?dl=1&token_hash=AAF8H3aS64fwU4T_dVBzhQ_OokjGNanvWFTYBZS3IjY4Vg',
                  method='auto',mode="wb",'test.zip')
    unzip('test.zip',exdir=tempdir())
    nnshp<-readShapeSpatial('Nanhaizhudao.shp')

    nanhai<-function()
{
  before_op <- par(no.readonly = TRUE)
  plot(nnshp,add=T)
  par(new=TRUE, oma=c(3.5,3.5,3.4,2.2)
      ,mar=c(2,3,2,1.5)
  )
  layout(matrix(c(0,0,0,0,0,0,
                  0,0,0,0,0,0,
                  0,0,0,0,0,0,
                  0,0,0,0,1,0,
                  0,0,0,0,0,0
  ),nrow=5,byrow=T))

  plot(plotr, xlim =c(1e5, 1.2e+6),ylim =c(1.5e5, 2.475e+06),
       breaks=breaks,box=F,
       legend=F,lwd = 1,axes=FALSE, frame.plot=F,
       col = topo.colors(length(breaks))) 
  par(before_op)
}  
if(as.numeric(dev.cur())>1){graphics.off()}
tiff('test.tif',
     width=12,height=8,units='cm',
     res=600,compression="lzw",pointsize=7)
  par(mar = c(3.5, 3.5, 3, 2.5)) 
  plot(plotr,col = topo.colors(length(breaks)),breaks=breaks,horizontal=F,
       legend.shrink = 0.8,legend.width = 0.6,legend.mar =7,cex=0.1,
       addfun=nanhai,
       xlim=c(xmin(plotr), xmax(plotr)), ylim=c(ymin(plotr), ymax(plotr)))
par(op)
dev.off()
if(as.numeric(dev.cur())>1){graphics.off()}

unlink(tempdir())
setwd(oldwd)

Currently ,I have two questions:

1. How can I modify the display scale of plot area to make sure that nnshp is entirely visible ?(avoid to adjust par parameters manually )

2. How can I adjust the location of sub plot graph precisely to nnshp ? I spent so many time to modify parameters and failed to get it.

3. When I want to get multiple plot (stack raster plot) just like the map demo made by ArcGis, I get the wrong result. How to get the multiple plot of it?

   plot(stack(plotr,plotr),col = topo.colors(length(breaks)),breaks=breaks,horizontal=F,
      legend.shrink = 0.8,legend.width = 0.6,legend.mar =7,cex=0.1,
      addfun=nanhai,
      xlim=c(xmin(plotr), xmax(plotr)), ylim=c(ymin(plotr), ymax(plotr)))