Efficiently classify raster with irregular polygons

Question

I have a project on my desk to classify an elevation model using different elevation breaks for each polygon in the area of interest.

In pseudo code:

FOR each polygon IN layer:
    CREATE_MASK FROM this polygon
    CLASSIFY elevmod WITH this polygon's LOOKUP_TABLE
    SAVE_AS elevmod_classified_this_polygonID

MERGE elevmod_classified*
CONVERT from 16bit float to 256 index colour palette

What is the most efficient way to do this?

The workflow I have now is wasteful of processing time and storage. Each classify step creates a new raster file for every polygon, and then all are combined into one at the end. The source raster is 650mb compressed and about 30 zone polygons (so 30 times some portion of 650mb for extra storage). This seems rather silly since I can classify a whole raster at run-time in a second using the standard symbology tools.

Plus the elevation breaks for each polygon haven't been finalized, so I'll need to run through this exercise a number of times before everyone is satisfied.

The classification elevation lookup tables resemble the following. Each number is the upper elevation boundary of that class:

polygonID, class1, class2, class3, ...
169,  850, 1100, 1300, 1450
170, 1000, 1250, 1400, 1600

The final classified raster has cell values of 1,2,3,4,5... corresponding to class1, class2, ... Though this too is open for change.

I have Qgis 1.7 and Arcinfo 10 with Spatial Analyst at my disposal (and Arcview3 for that matter).

Answer 1

One solution is available through the less than frequency operation.

To exploit this, create one polygon-based grid for each field. This will give a [class1] grid, a [class2] grid, ..., etc. The values of each will be constant within each polygon (and presumably all polygons are disjoint, for otherwise this problem is not well-framed). The less-than-frequency operation replaces each value in the DEM by the number of such grids whose values are less than that. Thus, any DEM value which is lower than the smallest value of {[class1], [class2], ...} at its location will get a result of 0; any DEM value which lies between the lowest and second lowest of these class cutpoint will get a result of 1; etc. That's precisely the desired reclassification (well, you might want to add 1 to each of them if you count starting at 1 instead of 0).

Here's a (hand-calculated) example. Let the DEM have values

1 2 3 3 2
3 2 2 3 2
4 3 2 2 1
5 3 1 1 0
6 4 2 2 1

Suppose there are two polygons designated 1 and 2, located like this:

. 1 1 . 1
1 1 1 1 1
1 . 2 1 1
1 2 2 2 .
2 2 2 2 2

(. stands for NoData: these are regions outside all polygons).

Let the cutpoint table have three classes with these values:

Polygon Class1 Class2 Class3
      1      1      2      4
      2      1      3      6

This gives rise to three cutpoint rasters, [Class1], [Class2], and [Class3], shown here side by side:

[Class1]      [Class2]      [Class3]
. 1 1 . 1     . 2 2 . 2     . 4 4 . 4
1 1 1 1 1     2 2 2 2 2     4 4 4 4 4
1 . 1 1 1     2 . 3 2 2     4 . 6 4 4
1 1 1 1 .     2 3 3 3 .     4 6 6 6 .
1 1 1 1 1     3 3 3 3 3     6 6 6 6 6

Applying the less-than-frequency operation (if I have computed it correctly) yields this grid of counts:

. 1 2 . 1
2 1 1 2 1
2 . 1 1 0
3 1 0 0 .
2 2 1 1 0

As a check, consider the elevations for polygon 2 only:

. . 2 . .
. 3 1 1 .
6 4 2 2 1

The cutpoint table asks to reclassify these values at 1, 3, and 6. Using less-than-frequency means the intervals will be from -infinity through 1, from just greater than 1 through 3, from just greater than 3 through 6, and strictly greater than 6. Thus, values of 1 should be reclassified to 0, values of 2 and 3 should be reclassified to 1, values of 4, 5, and 6 should be reclassified to 2, and any value greater than 6 should be reclassified to 3:

. . 1 . .
. 1 0 0 .
2 2 1 1 0

This is precisely what less-than-frequency has done.