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

        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)

1 Answer 1


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




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:


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.