24

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

csheth
  • 664
  • 1
  • 8
  • 22

1 Answers1

39

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)
Robert Hijmans
  • 10,683
  • 25
  • 35
  • 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 02 '16 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. – Robert Hijmans Apr 02 '16 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 02 '16 at 21:22
  • 1
    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. – Robert Hijmans Apr 03 '16 at 03: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 05 '16 at 10:09
  • The example has 9200 cells and runs in about a second --- which is slow, but I do not see how it would not work for 100 cells. – Robert Hijmans Apr 05 '16 at 15:58