19

I want to create two polygons.

  1. One is of the rectangular extents of my raster. I know I can obtain the extent of the raster using r <- raster("band5.tif") e <- extent(r) plot(e) gives me this: enter image description here

  2. How can I create ONE polygon of the boundary of the raster as shown below?

This is what my raster looks like:raster("band5.tif")

1 Answer 1

33

Here is an example.

library(raster)
# example data
x <- raster(system.file("external/test.grd", package="raster"))

To get the rectangular extent

e <- extent(x)
# coerce to a SpatialPolygons object
p <- as(e, 'SpatialPolygons')  

To get a polygon that surrounds cells that are not NA

# make all values the same. Either do
r <- x > -Inf
# or alternatively
# r <- reclassify(x, cbind(-Inf, Inf, 1))

# convert to polygons (you need to have package 'rgeos' installed for this to work)
pp <- rasterToPolygons(r, dissolve=TRUE)

# look at the results
plot(x)
plot(p, lwd=5, border='red', add=TRUE)
plot(pp, lwd=3, border='blue', add=TRUE)

Five years later: Nowadays I would use terra, which does this much faster.

library(terra)
z <- rast(system.file("external/test.grd", package="raster"))
pe <- as.polygons(ext(z))
pr <- as.polygons(z > -Inf)

plot(z)
plot(pe, lwd=5, border='red', add=TRUE)
plot(pr, lwd=3, border='blue', add=TRUE)
6
  • I can't seem to find documentation on the r <- r > -Inf part. What is it exactly doing? And how different is it from values(r)[!is.na(values(r))] <- 1 (which sets all locations that are not NA in r to 1).
    – csheth
    Apr 2, 2016 at 18:53
  • 3
    r > -Inf is basic R . Do c(1,NA,3,NA) > -Inf to see how it works. I have added an alternative (reclassify). Your alternative works but it is not a good one for large objects. Apr 2, 2016 at 19:08
  • I'm just waiting for my computer to process the rasterToPolygons() part, and thereupon I shall evaluate your answer. Thanks for the -Inf bit its going to be very useful!
    – csheth
    Apr 2, 2016 at 21:22
  • For a very large raster this will take a while, or may even fail. In that case you could consider first using aggregate as you would not see the difference anyway. Apr 3, 2016 at 3:59
  • Even with only 100 cells after running the aggregate() on r, it has not processed (computer equipped with 16 GB RAM). Problem is I need to keep the boundary as native as I can, since a DEM needs to be extracted from the resultant polygon. If I further down-sample my raster I will loose that boundary. Any work around?
    – csheth
    Apr 5, 2016 at 10:09

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.