Removing part of a raster image over a body of water

Question

I'm trying to make a map with leaflet that will ultimately be used to display various estimates interpolated over a geographic region.

I apologize for linking to the data. I wasn't clever enough to figure out how to import the .RData file from Google Docs. The raster shape I'm working with can be downloaded from https://drive.google.com/file/d/0Bw3leA1Ef_e5RmdKVDNYX0xmS2s/view?usp=sharing (This is an updated link with the Coord, rnew, and cleveland objects used below).

I'm able to get the map I want with:

library(tigris) 
library(leaflet)
library(raster)
library(magrittr)
library(dplyr)

load("Coord.Rdata")

rnew <- 
  rasterFromXYZ(
    Coord,
    crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  )

pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(rnew),
                    na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(rnew, colors = pal, opacity = 0.4)

Which produces

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.

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

ohio <- tracts(state = "39")

ohio_raster <- 
  attributes(ohio)$data[, c("INTPTLON", "INTPTLAT")] %>%
  setNames(c("x", "y")) %>%
  mutate(x = sub("^[-]", "", x) %>%
           as.numeric(),
         y = sub("^[+]", "", x) %>%
           as.numeric(),
         layer = 1) %>%
  as.matrix() %>%
  rasterFromXYZ(crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

but got the error

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

This sort of works, but there's no color, so I can't see the intersection

ohio_raster <- 
  raster(ohio, crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

leaflet() %>% addTiles() %>%
  addRasterImage(rnew, colors = pal, opacity = 0.4) %>%
  addRasterImage(ohio_raster, colors = "red", opacity = .4)

Warning message: In raster::projectRaster(x, raster::projectExtent(x, crs = sp::CRS(epsg3857))) :
'from' has no cell values

(notice: the map is zoomed out to include all of Ohio)

I also thought it might be possible to work out the intersection if I used the spatial package:

spdf <- 
  SpatialPointsDataFrame(Coord[, 1:2], 
                         Coord,
                         proj = CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

ohio_coord <- attributes(ohio)$data[, c("INTPTLON", "INTPTLAT")] %>%
  setNames(c("x", "y")) %>%
  mutate(x = sub("^[-]", "", x) %>%
           as.numeric(),
         y = sub("^[+]", "", x) %>%
           as.numeric(),
         layer = 1) %>%
  as.matrix() %>%
  as.data.frame()

spdf_ohio <- 
  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 <- merge(spdf, spdf_ohio)

leaflet() %>% addTiles() %>%
  addRasterImage(raster(spdf_merge), colors = pal, opacity = 0.4)

but I get

Warning message: In raster::projectRaster(x, raster::projectExtent(x, crs = sp::CRS(epsg3857))) : 'from' has no cell values

So, the actual question(s):

Is this even a reasonable approach to take for this problem? Is there a better approach that I'm missing?
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?

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.

EDIT:

I feel like I may be one step closer, but still not quite there.

I downloaded the Cleveland geography per @IvanSanchez 's recommendation. Simply plotting the polygon provides.

Which looks like the landspace I want to keep. That's a good step.

Now I tried intersecting that with my current raster via

cleveland <- 
  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() %>% addTiles() %>%
  addRasterImage(raster::intersect(rnew, cleveland), 
                 colors = pal, 
                 opacity = 0.4)

But instead of the intersection of the areas, I seem to have the union. Or something different.

I've updated the link above to a file that includes the Coord, rnew, and cleveland objects.

jbaums jbaums · Accepted Answer · 2016-04-30T00:09:38

You can mask the raster with a shapefile representing the land areas. You could use your Ohio tract shapefile, but Natural Earth's lakes data might simplify things. For example:

library(leaflet)
library(raster)
library(magrittr)
library(rgdal)
library(rgeos)

load('Coord.Rdata')  
rnew <- rasterFromXYZ(Coord, crs='+init=epsg:4326')
pal <- colorNumeric(c('#0C2C84', '#41B6C4', '#FFFFCC'), values(rnew),
                    na.color = 'transparent')

# Download and extract the Natural Earth data
download.file(
  file.path('http://www.naturalearthdata.com/http/',
            'www.naturalearthdata.com/download/10m/physical/ne_10m_lakes.zip',
            f <- tempfile())
unzip(f, exdir=tempdir())

# Read in the lakes shapefile    
lakes <- readOGR(tempdir(), 'ne_10m_lakes')

# Create a "land" shapefile by taking the difference between the 
# raster's extent and the lakes polys
land <- gDifference(as(extent(rnew), 'SpatialPolygons'), lakes)

# Mask rnew
rnew2 <- mask(rnew, land)

leaflet() %>% addTiles() %>%
  addRasterImage(rnew2, colors = pal, opacity = 0.4)

It doesn't perfectly follow the coastline, but it's not too bad.

Removing part of a raster image over a body of water

EDIT:

2 Answers