R raster plotting an image, draw a circle and mask

问题:

Code below plots an image and then plots circle on that image. I want to make all pixels that fall outside that circle black. How could I do that?

library(raster)
library(plotrix)
r1 <- brick(system.file("external/rlogo.grd", package="raster"))
width=50
height=40
x <- crop(r1, extent(0,width,0,height))
plotRGB(x)
circlex=20
circley=15
radius=10
draw.circle(circlex,circley,radius,border="blue")

回答1:

Look at the 'x'-object with str() and you see this:

..@ data    :Formal class '.MultipleRasterData' [package "raster"] with 14 slots
  .. .. ..@ values    : num [1:2500, 1:3] 255 248 221 199 198 210 221 190 104 79 ...
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : chr [1:3] "red" "green" "blue"

....so the 1:50 by 1:50 values are mapped to three columns. The X values are probably 0:2500 %% 50 and the y values are probably 0:2500 %/% 50

So remembering that the "origin" if upper left corner for raster objects but the lower left corner for plot functions, and so the y-value of 20 becomes 50-20 or 30, this gives you close to what you ask for (with apologies for putting the y-sequence first):

x@data@values[( ((1:2500 %/% 50 )- 30)^2 + ((1:2500 %% 50) - 20)^2 ) >=100, 1] <- 0
 x@data@values[( ((1:2500 %/% 50 )- 30)^2 + ((1:2500 %% 50) - 20)^2 ) >=100, 2] <- 0
 x@data@values[( ((1:2500 %/% 50 )- 30)^2 + ((1:2500 %% 50) - 20)^2 ) >=100, 3] <- 0
 plotRGB(x)
 draw.circle(20,20,10,border="blue")

The logic is that the criteria are of the form (x-dx)^2+(y-dy)^2 > r^2 where dx and dy are the center coordinates of the circle and r is the radius == 10.

EDITS after question changed:

For each color layer, the code with named parameters would be similar to this one that makes the darkest "red". This gives a roughly circular mask although getting the centers to line up is not handled correctly:

x@data@values[( ((1:(height*width) %/% (height*5/4) )- (height-circley*5/4) )^2 + 
            ((1:(height*width) %% width) -       circlex )^2 ) >= radius^2, 1] <- 0

Further experimentation delivers this which seems pretty close:

x@data@values[( ((1:(height*width) %/% (height) )- (height-circley) *5/4)^2 + 
                 ((1:(height*width) %% width) -       circlex )^2 ) >= radius^2, 1] <- 0
plotRGB(x, asp=5/4, addfun=function() draw.circle(circlex,circley,radius,border="blue") )

Obviously you could substitute the width/height scaling factor for the new aspect ratio everywhere that 5/4 appears.

回答2:

Here it is a more generic solution that also works for rasters with different coordinate systems.

The function rasterToPoints() converts the raster to points. Following your example, head(rasterToPoints(x)) returns the following:

> head(rasterToPoints(x))
       x    y red green blue
[1,] 0.5 39.5 255   255  251
[2,] 1.5 39.5 204   205  199
[3,] 2.5 39.5 171   172  164
[4,] 3.5 39.5 157   159  148
[5,] 4.5 39.5 162   164  151
[6,] 5.5 39.5 187   189  176

We then need to find which points fall outside the circle and set their values to zero:

is_outside_circ = (rasterToPoints(x)[,1] - circlex)^2 + (rasterToPoints(x)[,2] - circley)^2 >= radius^2
x@data@values[ is_outside_circ,] <- 0
plotRGB(x)
draw.circle(circlex,circley,radius,border="blue")

Result: black points outside circle