title: Removing part of a raster image over a body of water
body: <p>I'm trying to make a map with <code>leaflet</code> that will ultimately be used to display various estimates interpolated over a geographic region.</p>
  
  <p>I apologize for linking to the data.  I wasn't clever enough to figure out how to import the <code>.RData</code> file from Google Docs.  The raster shape I'm working with can be downloaded from <a href="https://drive.google.com/file/d/0Bw3leA1Ef_e5RmdKVDNYX0xmS2s/view?usp=sharing" rel="nofollow noreferrer">https://drive.google.com/file/d/0Bw3leA1Ef_e5RmdKVDNYX0xmS2s/view?usp=sharing</a>
  (This is an updated link with the <code>Coord</code>, <code>rnew</code>, and <code>cleveland</code> objects used below).</p>
  
  <p>I'm able to get the map I want with:</p>
  
  <pre><code>library(tigris) 
  library(leaflet)
  library(raster)
  library(magrittr)
  library(dplyr)
  
  load("Coord.Rdata")
  
  rnew &lt;- 
    rasterFromXYZ(
      Coord,
      crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
    )
  
  pal &lt;- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(rnew),
                      na.color = "transparent")
  
  leaflet() %&gt;% addTiles() %&gt;%
    addRasterImage(rnew, colors = pal, opacity = 0.4)
  </code></pre>
  
  <p>Which produces 
  <a href="https://i.stack.imgur.com/bewbG.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/bewbG.jpg" alt="enter image description here"></a></p>
  
  <p>Now, as beautiful as that is, it doesn't make much sense to me to display estimates out over Lake Erie.  I'd like to remove the parts of the raster image over the Lake, or set the to transparent...something that doesn't leave the impression that I'm calculating risk for marine wildlife.</p>
  
  <p>I figured that if I could find the latitude and longitude of the State of Ohio, I might be able to use just the intersection, thereby removing the points over the Lake.  I tried putting it together using</p>
  
  <pre><code>ohio &lt;- tracts(state = "39")
  
  ohio_raster &lt;- 
    attributes(ohio)$data[, c("INTPTLON", "INTPTLAT")] %&gt;%
    setNames(c("x", "y")) %&gt;%
    mutate(x = sub("^[-]", "", x) %&gt;%
             as.numeric(),
           y = sub("^[+]", "", x) %&gt;%
             as.numeric(),
           layer = 1) %&gt;%
    as.matrix() %&gt;%
    rasterFromXYZ(crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
  </code></pre>
  
  <p>but got the error </p>
  
  <blockquote>
    <p>Error in rasterFromXYZ(., crs = "+init=epsg:4326 +proj=longlat
    +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") :    x cell sizes are not regular</p>
  </blockquote>
  
  <p>This sort of works, but there's no color, so I can't see the intersection</p>
  
  <pre><code>ohio_raster &lt;- 
    raster(ohio, crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
  
  leaflet() %&gt;% addTiles() %&gt;%
    addRasterImage(rnew, colors = pal, opacity = 0.4) %&gt;%
    addRasterImage(ohio_raster, colors = "red", opacity = .4)
  </code></pre>
  
  <blockquote>
    <p>Warning message: 
    In raster::projectRaster(x, raster::projectExtent(x, crs = sp::CRS(epsg3857))) :<br>
    'from' has no cell values</p>
  </blockquote>
  
  <p>(notice: the map is zoomed out to include all of Ohio)
  <a href="https://i.stack.imgur.com/7yDUq.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/7yDUq.jpg" alt="enter image description here"></a></p>
  
  <p>I also thought it might be possible to work out the intersection if I used the spatial package:</p>
  
  <pre><code>spdf &lt;- 
    SpatialPointsDataFrame(Coord[, 1:2], 
                           Coord,
                           proj = CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
  
  ohio_coord &lt;- attributes(ohio)$data[, c("INTPTLON", "INTPTLAT")] %&gt;%
    setNames(c("x", "y")) %&gt;%
    mutate(x = sub("^[-]", "", x) %&gt;%
             as.numeric(),
           y = sub("^[+]", "", x) %&gt;%
             as.numeric(),
           layer = 1) %&gt;%
    as.matrix() %&gt;%
    as.data.frame()
  
  spdf_ohio &lt;- 
    SpatialPointsDataFrame(ohio_coord[, 1:2], 
                           ohio_coord,
                           proj = CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
  
  spdf_merge &lt;- merge(spdf, spdf_ohio)
  
  leaflet() %&gt;% addTiles() %&gt;%
    addRasterImage(raster(spdf_merge), colors = pal, opacity = 0.4) 
  </code></pre>
  
  <p>but I get </p>
  
  <blockquote>
    <p>Warning message: In raster::projectRaster(x, raster::projectExtent(x,
    crs = sp::CRS(epsg3857))) :   'from' has no cell values</p>
  </blockquote>
  
  <p><a href="https://i.stack.imgur.com/qRHfn.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/qRHfn.jpg" alt="enter image description here"></a></p>
  
  <p>So, the actual question(s):</p>
  
  <ol>
  <li>Is this even a reasonable approach to take for this problem?  Is there a better approach that I'm missing?</li>
  <li>How can I go from latitude and longitude to a raster image with a value to give the shading so that I can see what I'm working with?</li>
  </ol>
  
  <p>Even if I don't get a direct answer to the problem at hand, I'll be ecstatic to receive some guidance on where to find examples of similar problems.  I haven't found anything similar yet, which is probably a sign not knowing the right search terms.</p>
  
  <h2>EDIT:</h2>
  
  <p>I feel like I may be one step closer, but still not quite there.</p>
  
  <p>I downloaded the Cleveland geography per @IvanSanchez 's recommendation.  Simply plotting the polygon provides.</p>
  
  <p><a href="https://i.stack.imgur.com/eDeVT.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/eDeVT.png" alt="enter image description here"></a></p>
  
  <p>Which looks like the landspace I want to keep.  That's a good step.</p>
  
  <p>Now I tried intersecting that with my current raster via</p>
  
  <pre><code>cleveland &lt;- 
    readOGR(dsn = "cleveland_ohio.osm2pgsql-shapefiles",
            layer = "cleveland_ohio_osm_polygon",
            p4s = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
  
  leaflet() %&gt;% addTiles() %&gt;%
    addRasterImage(raster::intersect(rnew, cleveland), 
                   colors = pal, 
                   opacity = 0.4)
  </code></pre>
  
  <p>But instead of the intersection of the areas, I seem to have the union.  Or something different.  </p>
  
  <p><a href="https://i.stack.imgur.com/9pKmM.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9pKmM.jpg" alt="enter image description here"></a></p>
  
  <p>I've updated the link above to a file that includes the <code>Coord</code>, <code>rnew</code>, and <code>cleveland</code> objects.</p>