The script shows a R example of geospatial masking in the Maumee River basin, plotting precipitation of a gridded Lambert netCDF file with a watershed in polygon shapefile on Google Map.

Library

xlib  <- c("tidyverse","ncdf4","raster", "ggmap", "rgdal", "rgeos")
lib   <- lapply(xlib, library, character.only = TRUE) # load the required packages

Read gridded NetCDF file

geoFile <- "geo_em.d02.nc"
fname   <- paste0("Precip/", "201401010000.PRECIP_FORCING.nc")
fnc     <- nc_open(fname)

Get Coordinate

lon <- ncvar_get(fnc, varid = "lon")
lat <- ncvar_get(fnc, varid = "lat")
# get varaible precipitation rate
PR  <- ncvar_get(fnc, varid = "precip_rate")

Get Projection

library(rwrfhydro)
proj4 <- GetProj(geoFile)

Assign coordinates

spPR              <- data.frame(lon = as.vector(lon), lat = as.vector(lat), conc = as.vector(PR))
coordinates(spPR) <- ~ lon+lat
proj4string(spPR) = CRS(proj4)

Rasterization

r             <- raster(ncols = ncol(lat),nrows = nrow(lat))
extent(r)     <- extent(spPR)
rPR           <- rasterize(coordinates(spPR), r, spPR$conc)

Read Maumee River Shapfile

rMaumee       <- readOGR(dsn = "MaumeeRiverShapefile", layer = "Maumee_Output")
## OGR data source with driver: ESRI Shapefile 
## Source: "MaumeeRiverShapefile", layer: "Maumee_Output"
## with 1 features
## It has 35 fields
## Integer64 fields read as strings:  OBJECTID GLHDID HydroID COMID WBAREACOMI FCODE GNIS_NBR STRAHLER SHREVE
rMaumee.proj  <- spTransform(rMaumee, CRSobj = CRS(proj4string(rPR)))

Plot raster layers and add Maumee polygon

#par(mar = c(0.1, 0.1, 0.1, 0.1 ))
plot(rPR*1e4, asp = 1)
plot(rMaumee.proj, add = T)

### Raster to Polygon

r   <- projectRaster(rPR, crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
rtp <- rasterToPolygons(r)

Google Basemap

bm <- ggmap(get_map(location = bbox(rtp)))
bm +  geom_polygon(data = rtp,
                   aes(x = long, y = lat, group = group, 
                       fill = rep(rtp$layer*1e4, each = 5)), 
                   size = 0, 
                   alpha = 0.5)  + 
  scale_fill_gradientn(expression("PR\n"~10^{-4}), 
                       colors = topo.colors(255))