Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have simple question: 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?

share|improve this question
add comment

2 Answers 2

up vote 8 down vote accepted

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

share|improve this answer
add comment

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.

share|improve this answer
    
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) –  RobertH 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) –  RobertH Oct 5 '12 at 6:10
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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