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 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).

share|improve this question
1  
how many classes do you typically have? –  Matthew Snape Mar 22 '12 at 22:19
1  
If each polygon has the same number of classes, then less than frequency will do the job. Otherwise, you can pad the table with artificial (very high) classes and do it. This ought to be efficient unless the number of classes is large (greater than 50 or so). Does this look like a promising approach? –  whuber Mar 23 '12 at 3:54
    
@MatthewSnape, not very many classes perhaps 8, and once settled they'll be stable. –  matt wilkie Mar 23 '12 at 15:35
    
@whuber Thanks, it's the first time I've looked at that tool. I don't see how to apply Less Than Frequency in this case as all inputs are raster. I could convert the polygons to raster, but then each area (poly) would be a constant value. Maybe I need to post some example data or graphics(?) –  matt wilkie Mar 23 '12 at 15:40

1 Answer 1

up vote 2 down vote accepted

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.

share|improve this answer
    
Thanks. This produces the desired result. I confess to not being able to wrap my head around precisely why it works. Even though I've stepped carefully through the process and explanation several times, and I understand each individual step, when I stand back and look at the whole thing my mind goes "huh?". ;-) –  matt wilkie Mar 27 '12 at 21:37
    
I'm glad you were able to make this solution work, Matt. I bet it's fairly efficient, too. Because this is a purely local operation--knowing the stack of values above a single cell suffices to compute the result--if you want to understand it better, just think about the analog for pure numbers. When you reclassify a numeric value x according to set of cutpoints x1 < x2 < ... (think of a legend for a layer), the class x belongs to is determined by the number of cutpoints less than it: exactly like classifying a layer for mapping. That's all that's going on here. –  whuber Mar 27 '12 at 22:22

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.