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?

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Seen
  • 2,205
  • 7
  • 26
  • 37

2 Answers2

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

whuber
  • 69,783
  • 15
  • 186
  • 281
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.

Robert Hijmans
  • 10,683
  • 25
  • 35
  • is there anyway just perform spatial query rather than cropping the image? Cropping might be very slow for processing. – Seen Sep 23 '12 at 03: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 06: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 05 '12 at 06:10