Hfrhyu

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.

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?

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.

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)

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')

With your procedure I reproduced the same issue:

p

The solution is to add coord_equal()

p + coord_equal()

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