Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Change raster values using spatial polygons

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.

like image 525
loki Avatar asked Jul 27 '16 09:07

loki


People also ask

How do you convert the pixels of a raster to a polygon?

Use the Raster to Point tool to convert each pixel to a point. Navigate to ArcToolbox > Conversion Tools > From Raster > Raster to Point. These points are later used to label polygons.

How do you extract raster from a polygon?

In the Clip section, under Clipping Geometry/Raster, browse to the desired polygon features, and check the Use Input Features for Clipping Geometry check box. Click OK > OK. Export the raster to save it permanently.


1 Answers

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

like image 86
loki Avatar answered Oct 09 '22 10:10

loki