Change raster values using spatial polygons
Asked Answered
A

1

8

To change raster values under SpatialPoints you can simply use [ .

r <- raster(system.file("external/test.grd", package="raster"))
rTp <- rasterToPoints(r, spatial = T)

set.seed(666)
rTpS <- rTp[sample(1:length(rTp), 500),]
plot(r)
plot(rTpS, add = TRUE, pch = ".")

Now, the raster values beneath the spatial points are changed like this:

r1 <- r
r1[rTpS] <- 99999
plot(r1)

Raster with SpatialPoints (left) and raster with changed values (right): Raster with SpatialPoints (left) and raster with changed values (right)

Thus, I assumed that polygons work the same way. Therefore generate some polygons from this raster data set:

rst <- as.integer(stretch(r, 0, 10))
rTpol <- rasterToPolygons(rst, dissolve = T)
rTpol <- rTpol[rTpol$layer > 3,]

plot(r)
plot(rTpol, add = T)

rst[rTpol[,]] <- 100
plot(rst)

# also tried 
# rst[rTpol] <- 100

However, when i Try to change the raster values below some SpatialPolygons, somehow the ID (or something) of the polygons. Yet, it should be 100.

Raster with SpatialPolygons (left) and raster with strangely changed values (right): Raster with SpatialPolygons (left) and raster with strangely changed values (right)

Therefore, my question is how to change raster values within polygons.

Aramen answered 27/7, 2016 at 9:31 Comment(0)
A
7

Edit (after two years):

This has been fixed as promised by Mr. Hijmans in the comments.
Thanks for this easy and awesome and convenient way of working with geo data!

Original Answer:

Finally, I found the answer myself.

The function rasterize() helps with this problem.

r <- raster(system.file("external/test.grd", package="raster"))
rst <- as.integer(stretch(r, 0, 10))
plot(rst)
rTpol <- rasterToPolygons(rst, dissolve = TRUE)
rTpol <- rTpol[rTpol$layer > 3,]
plot(rTpol)

rpol <- rasterize(rTpol, rst, field = 100, update = TRUE)
plot(rpol)

rasterize() puts the field = value in the rastercells generated from the polygons, while update = TRUE takes the raster and updates all cells wich correspond to the polygons with the field value.

The result looks like this: enter image description here

Aramen answered 27/7, 2016 at 10:45 Comment(2)
That is a good solution, but your original solution was more elegant --- and should have worked. This bug has been fixed in (forthcoming) raster version 2.5-16Helmer
@RobertH: Do you know if the update to the raster package is already on the way?Aramen

© 2022 - 2024 — McMap. All rights reserved.