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?

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.

  • is there anyway just perform spatial query rather than cropping the image? Cropping might be very slow for processing. – Seen Sep 23 '12 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 Sep 23 '12 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) – Robert Hijmans Sep 24 '12 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 Sep 24 '12 at 14:54
  • It is a function in the raster package. ?extract shows that (at least for me it does) – Robert Hijmans Oct 5 '12 at 6:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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