R – Overlay Shapefile on Multiband Raster Using ggRGB and ggplot2 in R

coordinate systemggplot2rraster

I am trying to use ggplot2 to map a polygon feature class over a raster brick (created from multiple Landsat band images). I'm able to map the polygons and the raster separately, but when I try to have them together (following the last code chunk at http://bleutner.github.io/RStoolbox/rstbx-docu/ggRGB.html RStudio crashes (even after restarting the computer).

ggRGB attempt:

fallowFields <- rgdal::readOGR(dsn = fgdb, layer = 'Fallowing2018')
fallowFieldsWGS84 <- spTransform(fallowFields, crs(L8_20180420_brick))

fallowFieldsDataFrame <- fallowFieldsWGS %>% fortify

landsatAprilCropped <- raster::crop(L8_20180420_brick, raster::extent(fallowFields))

p <- ggplot()
p + ggRGB(img = landsatAprilCropped,
          r = 3,
          g = 2,
          b = 1,
          stretch = 'hist',
          ggLayer = TRUE) +
geom_polygon(col = 'Black',
             fill = 'grey70',
             data = fallowFieldsDataFrame,
             aes(x = long, y = lat, group = group))

replacing the data frame (fallowFieldsDataFrame) with the shapefile (fallowFieldsWGS84) works, but it warps the Landsat image.

enter image description here

The cropped image and shapefile are in the same projection, but the shapefile seems to be warping the raster. Is there a way to plot the shapefile over the raster without warping it?

enter image description here

enter image description here

It doesn't get warped with plotRGB, but I was hoping to use ggRGB so that I could make a more elegant map with gglot2.

enter image description here

plotRGB(landsatAprilCropped, 
        r = 3, g = 2, b = 1, 
        stretch = 'hist', axes = TRUE,
        main = 'Fallow Fields Overlaid on Landsat Image: 4/20/2018')

plot(bardFallowFields, col = 'transparent', border = 'red', add = TRUE)

Best Answer

You miss a very important part. First, an example:

library(raster)
library(ggplot2)
library(RStoolbox)

data(lsat)

lsat

poly <- as(extent(622999,625000,-417000,-415000),'SpatialPolygons')

poly <- fortify(poly)

names(poly)[1:2] <- c('x','y')

e <- extent(lsat)


p <- ggplot() + ggRGB(img = lsat,
            r = 3,
            g = 2,
            b = 1,
            stretch = 'hist',
            ggLayer = T) + 
  xlim(e[1:2]) + ylim(e[3:4]) + 
  geom_polygon(aes(x,y), alpha = 0.9, data = poly)

Using ggRBG() output is:

ggRGB(img = lsat,
      r = 3,
      g = 2,
      b = 1,
      stretch = 'hist')

enter image description here

With your procedure I reproduced the same issue:

p

enter image description here

The solution is to add coord_equal()

p + coord_equal()

enter image description here

I recomend you to add ylim() and xlim() with the extent of raster object

Related Question