ggproj <- function(r, lat, lon, tolerance=0.001, plot.borders=TRUE, border.color='black') {
    #r, lat, lon must be rasterLayers
    #Adapted from
    #http://stackoverflow.com/questions/43612903/how-to-properly-plot-projected-gridded-data-in-ggplot2
    if (!compareRaster(r, lat, lon)) {stop('Inputs should have the same size and projection')}
    require(ggplot2) #needs ggplot >= 2.2.1.9000
    require(sf)
    require(sp)
    crs <- projection(r)

    coords = data.frame(lat = values(lat), lon = values(lon))
    spPoints <- SpatialPointsDataFrame(coords,
                                      data = data.frame(data = values(r)),
                                      proj4string = CRS("+init=epsg:4326"))
    orig_grid = spTransform(spPoints, crs)
    polys = as(SpatialPixelsDataFrame(orig_grid, orig_grid@data, tolerance = tolerance),"SpatialPolygonsDataFrame")
    polys_sf = as(polys, "sf")

    gg <- ggplot(polys_sf) +
        geom_sf(aes(fill = data), color = "transparent") +
        theme_bw() +
        theme(panel.ontop=TRUE, panel.background=element_blank()) +
        scale_x_continuous(expand = c(0, 0)) +
        scale_y_continuous(expand = c(0, 0))

    if (plot.borders) {
        require(maps)
        borders <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))
        gg <- gg +
               geom_sf(data = borders, fill = "transparent", color = border.color) +
               coord_sf(crs = st_crs(crs),
                        xlim = st_bbox(polys_sf)[c(1,3)],
                        ylim = st_bbox(polys_sf)[c(2,4)]
                )
    }
    return(gg)
}