title: Plotting points and multipolygon objects in R with sf
body: <p>I've been trying to plot both <code>points</code> and <code>multipolygon</code> <code>sf</code> objects with no success. Following is a reproducible example of my problem.</p>
  
  <pre><code>library(sf)
  library(magrittr)
  library(RColorBrewer)
  
  nc &lt;- st_read(system.file("shape/nc.shp", package = "sf")) %&gt;% 
      st_transform(crs = 4326)
  points &lt;- data.frame(p = seq(15, 75, 15), 
                       long = c(-85, -80, -78, -75, -82), 
                       lat = c(34, 36, 37, 38, 35)) %&gt;% 
          st_as_sf(coords = c('long', 'lat'), crs = 4326) 
  points$p_cut &lt;- cut(points$p, seq(0, 100, 20))
  
  #plot1
  plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
       pal = brewer.pal(5, 'Paired'))
  plot(st_geometry(nc), axes = TRUE, graticule = TRUE, add = TRUE)
  </code></pre>
  
  <p><img src="https://i.imgur.com/Z7NzPgl.png" alt=""></p>
  
  <p>The <code>multipolygon</code> doesn't appear but if I inverse the order:</p>
  
  <pre><code>#plot2
  plot(st_geometry(nc), axes = TRUE, graticule = TRUE)
  plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
       pal = brewer.pal(5, 'Paired'), add = TRUE)
  </code></pre>
  
  <p><img src="https://i.imgur.com/Nc5sjbs.png" alt="">
  Two points and the key are missing from the final plot.</p>
  
  <p>How can I overcome these issues? I thought of changing the <code>bbox</code> of the <code>multipolygon</code> object in order to make the missing points appear in the second plot but I couldn't find a way for this. </p>
  
  <p>Any help is appreciated.</p>
  
  <p><strong>EDIT.</strong></p>
  
  <p>Modified the epsg of the <code>multipolygon</code> in order to match the <code>points</code> epsg. The issues remains though.</p>
  
  <p><strong>EDIT2.</strong>
  It's weird but running the code randomly seems to produce the "correct plot". I used this particular order (copied from the console):</p>
  
  <pre><code>&gt; plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
  +      pal = brewer.pal(5, 'Paired'))
  &gt; plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
  +      pal = brewer.pal(5, 'Paired'), add = TRUE)
  &gt; plot(st_geometry(nc), axes = TRUE, graticule = TRUE, add = TRUE)
  &gt; #plot2
  &gt; plot(st_geometry(nc), axes = TRUE, graticule = TRUE)
  &gt; plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
  +      pal = brewer.pal(5, 'Paired'), add = TRUE)
  &gt; # plot1
  &gt; plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
  +      pal = brewer.pal(5, 'Paired'))
  &gt; plot(st_geometry(nc), axes = TRUE, graticule = TRUE, add = TRUE)
  &gt; plot(points['p_cut'], axes = TRUE, graticule = TRUE, pch = 16, 
  +      pal = brewer.pal(5, 'Paired'), add = TRUE)
  </code></pre>
  
  <p>and got this:
  <img src="https://i.imgur.com/pfmfEUz.png" alt=""></p>
  
  <p>One of the missing points appears (-85,34) but there are repeated points. Also the key is not shown properly. It seems too that the location of the repeated points and the size of the multipolygon depend on the size of the plot viewer.</p>