Given a bounding box, how can I count the number of specific values (say, I am interested in the number of value == 1
) in a raster in 1) arcpy, or 2) raster package in R?
2 Answers
In R
, use crop
to extract the values and (e.g.) table
to count them.
As an example, let's create a 1 degree grid covering the globe:
library(raster)
x.raster <- raster(outer(179:0, 0:359, `+`), xmn=-180, xmx=180, ymn=-90, ymx=90)
The bounding box is converted to an extent
object in order to use crop
:
y.extent <- extent(cbind(c(-125,20), c(-60,50)))
y.raster <- crop(x.raster, y.extent)
Having done that, tabulation is straightforward:
table(getValues(y.raster))
In this output the first row lists the values and the second lists their corresponding counts:
165 166 167 ... 257 258
1 2 3 ... 2 1
As a check we can plot the raster and the extent:
plot(x.raster)
plot(y.extent, add=T)
A minor addition: you could also use (memory-safe) function "freq":
Following the answer by whuber:
library(raster)
x.raster <- raster(outer(179:0, 0:359, '+'), xmn=-180, xmx=180, ymn=-90, ymx=90)
y.extent <- extent(cbind(c(-125,20), c(-60,50)))
y.raster <- crop(x.raster, y.extent)
But now do:
freq(y.raster)
It only matters for very large objects (raster on file). 'freq' returns a two-column matrix (value/count) whereas 'table' returns a table.
-
is there anyway just perform spatial query rather than cropping the image? Cropping might be very slow for processing.– SeenCommented Sep 23, 2012 at 3:49
-
1Cropping ought to be very fast. As a test, I timed the
crop
operation for a one minute grid covering the earth: it has 10800 rows and 21600 columns (233,280,000 cells). The crop was executed in 1.36 seconds total elapsed time.– whuberCommented Sep 23, 2012 at 21:53 -
I also think it should be fast, but you could compare the above with this spatial query:
v <- extract(x.raster, y.extent)
followed bytable(v)
Commented Sep 24, 2012 at 6:02 -
1Robert, what library is
extract
from? It's not part ofraster
and theR
help system (??
) does not find any function with this name, either.– whuberCommented Sep 24, 2012 at 14:54 -
It is a function in the raster package. ?extract shows that (at least for me it does) Commented Oct 5, 2012 at 6:10