11

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 2

10

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)

Map

3

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.

5
  • is there anyway just perform spatial query rather than cropping the image? Cropping might be very slow for processing.
    – Seen
    Commented Sep 23, 2012 at 3:49
  • 1
    Cropping 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.
    – whuber
    Commented 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 by table(v) Commented Sep 24, 2012 at 6:02
  • 1
    Robert, what library is extract from? It's not part of raster and the R help system (??) does not find any function with this name, either.
    – whuber
    Commented 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.