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 generated zonal statistics using feature zone data on a value raster. This shapefile contains about 20000 polygons and the analysis went without a hitch.

While testing quality, I ran the same operation on 10 polygons taken from the same data. The values in zonal max changed this time. I compared the values by overlaying the polygons on the value raster and found that the pixel values belong in the polygon, but just not the maximum pixel.

I am aware that an internal raster conversion takes place for computing stats. I want to know what exactly happens during the process of zonal statistics and whether sample size affects the output.

I have ArcGIS 10 with an Arcview license and spatial analyst. Any guidance is highly appreciated.

share|improve this question
add comment

2 Answers

up vote 11 down vote accepted

Zonal stats ultimately is a raster-on-raster operation. When the polygons are provided in vector format, they will automatically be converted to raster format. I would hope that this would be the same cellsize and registered with the value raster, but I'm not sure of that. (It would take some difficult reverse engineering to check. Possibly the polygons are converted using whatever spatial analysis properties are in effect and then later are resampled for the zonal stats operation. That would be silly, but it's not hard to see how the software might be engineered that way.)

(The conversion to raster format explains why, by the way, you cannot directly perform zonal stats with overlapping polygons. Each raster cell can represent at most one polygon at a time.)

Historically, Spatial Analyst has had a hard time with zonal stats. Some relatively recent versions (as recently as 9.3, I recall) raised a lot of questions on the ESRI forums concerning wrong answers, missing values, and outright failures when more than about 2000 polygons were involved. Some of these (if not all) I suspect may have been due to misunderstandings about how to use zonal stats and what it does. In particular, the conversion from vector polygon to raster assigns a cell to a polygon if and only if the cell's center is contained in the polygon. (How the decision is made when the center lies on the polygon's boundary is another undocumented mystery requiring some tedious reverse engineering.) When a tortuous polygon manages to thread its way around all the nearby cell centers, that polygon simply disappears from the results altogether!

There's even more going on under the hood and it, too, has been changing from version to version of the software. When on-the-fly projection is involved, the polygons might or might not be reprojected before conversion to raster. So might the raster itself, for that matter. Precisely what happens may depend on the environment: whether this is in the Raster Calculator, the command line, a Python script, a tool, or whatever.

Being aware of all this is more than half the battle, because we can protect ourselves from the software. The best way is to create a raster of the polygons that has the same cellsize, with registered cells, as the value raster, and has the same projection. Disable any automatic resampling or on-the-fly projection. In this way you minimize the invisible "help" the software is giving you and have a reasonable chance of knowing precisely what the inputs to the calculation are. If you do that and still find discrepancies, then (IMHO) you have ample basis for a formal bug report.

(By following these procedures religiously, I have had no problems with zonal stats in any version of Spatial Analyst (1.0 through 9.3), even with very large grids and large numbers of polygons. For example, I obtained correct values for a set of Census blocks throughout a large US state at a resolution of 10 meters. (That would have been over 10^5 polygons, I recall.) However, I sometimes take extraordinary steps to preclude problems: I accomplished that Census block calculation by automatically breaking the state into its 50 or so counties, performing the calculation for each county, and then reassembling the results. Among other things, this allowed parallelization of the computation. It also allowed me to handle the tiny polygons, which weren't captured during the rasterization process, in a separate step whereby the value grid was queried at the polygon centroids. These values served as surrogates for the zonal means being computed.)

Incidentally, cell size ("sample size"?) affects some output in important ways. When you resample the value raster,

  1. The number of value cells within a polygon will change.

  2. Individual values are interpolated, causing a slight change in their overall statistical distribution. This can include changes in maxima and minima, which tend to shift towards the mean.

  3. Statistics related to commonality and frequency (such as a mode) can change entirely.

  4. Some geometric summaries, especially perimeter, can change arbitrarily.

  5. Many fancier statistical calculations, such as standard errors of estimates or regressions, can be entirely misleading due to the arbitrary changes in numbers of cells (see #1).

share|improve this answer
    
Thank you so much for your patient reply. I am still baffled about how the result might change for the same raster when processed with a smaller number of polygons.I used a cell size of 0.001 for raster generation.What confuses me is, if the same polygon is being rasterized internally and same number of pixels are contained in it, why does zonal stats show a value different than previous ones. –  Janus Jul 28 '12 at 20:02
    
I processed 20000 polygons of the whole country initially and then picked up 5 random polygons to check if the results are same. The unique pixel values falling in polygon A are, say, 314,624,456,234. Zonal stats returned 314 as max for polygon A, which is the majority but not max. This was the case when I processed it along with 20000 other polygons. Interestingly, when the same polygon was processed with just 5 other polygons, it returned 410. Both of them are incorrect. –  Janus Jul 28 '12 at 20:10
1  
Further, after reading your reply, I changed projections just to check and it returns different values for different projections. For UTM, it returned 624, for LCC it returned 456. Now, I am even more confused :). When the projection of both zones and value raster is same, should it not return the same maximum each time?. And in above case, why would it not return 624. I hope I am not doing anything wrong. –  Janus Jul 28 '12 at 20:26
    
Reprojection necessarily involves resampling and interpolation. When the maximum is isolated and surrounded by substantially smaller values, the interpolation can cause these kinds of discrepancies. The best procedure is to leave the value grid alone so its values don't get resampled and register the zone grid with the value grid. Then, assuming zonal stats is working correctly, you will get consistent, predictable results. –  whuber Jul 28 '12 at 20:47
    
Thank you so much :) –  Janus Jul 28 '12 at 21:01
add comment

Hey I'm trying to get an Ersi "Idea" promoted so they'll fix some of these issues, specifically the overlapping polygon issue.

Help out if you can. http://ideas.arcgis.com/ideaView?id=0873000000086QNAAY

share|improve this answer
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.