Can HRA Reclassification Percentages be Derived from pygeoprocessing.zonal_statistics outputs?

Background

Among its outputs, HRA produces a raster where the calculated floating-point risk score is reclassified into the following classifications:

1. 0 representing no risk
2. 1 representing low risk
3. 2 representing medium risk
4. 3 representing high risk

The current iteration of HRA does a lot of extra work to tally up the counts of pixels under each AOI polygon to determine the percentage of pixels under the polygon that are in each risk classification category because this output is not directly supported by HRA.

For reference, the structure of the return dictionary from zonal_statistics is:

{FID: {
    'min': 0,
    'max': 1,
    'sum': 1.7,
    'count': 3,
    'nodata_count': 0,
}}

Thus the question becomes: can the breakdown of high/medium/low classifications be derived from the outputs of zonal_statistics?

Mathy Approach: this can't be done with just these values

At first glance, we appear to be able to construct a system of linear equations:

0*x_0 + 1*x_1 + 2*x_2 + 3*x_3 = stats_{sum}

x_0 + x_1 + x_2 + x_3 = stats_{count}

It's worth noting that althoug this appears plausible, it is impossible to find a single solution to all of the unknowns: to solve 4 unknowns, we need 4 equations. Thus, solving a system of linear equations is only possible when there are 2 or fewer unknowns.

But in the general case, we end up with 2 unknowns that cannot be accounted for, effectively a whole plane of possible solutions

Demonstration

By way of demonstration, I have added a program here, test_lulc_pixelcounts.py, that:

1. Generates a set of 3 random integers representing pixel counts under a hypothetical area of interest
2. Uses a brute-force method to determine if there are any alternate solutions.

Example output:

In this example, this hypothetical polygon has:

• 0 pixels with risk classification 0
• 11 pixels with risk classification 1
• 62 pixels with risk classification 2
• 4 pixels with risk classification 3

Given the sum of these risk classifications ( 11+(2*62)+(3*4)=147 ) and the count of valid pixels (11+62+4=77), we can see that there are at least 36 other possible solutions:

Testing 0, 11, 62, 4, 36 solutions
    Solutions:
        0 7 70 0
        0 8 68 1
        0 9 66 2
        0 10 64 3
        0 11 62 4
        0 12 60 5
        0 13 58 6
        0 14 56 7
        0 15 54 8
        0 16 52 9
        0 17 50 10
        0 18 48 11
        0 19 46 12
        0 20 44 13
        0 21 42 14
        0 22 40 15
        0 23 38 16
        0 24 36 17
        0 25 34 18
        0 26 32 19
        0 27 30 20
        0 28 28 21
        0 29 26 22
        0 30 24 23
        0 31 22 24
        0 32 20 25
        0 33 18 26
        0 34 16 27
        0 35 14 28
        0 36 12 29
        0 37 10 30
        0 38 8 31
        0 39 6 32
        0 40 4 33
        0 41 2 34
        0 42 0 35

Conclusion

The only way to be completely sure about the per-class breakdown per AOI polygon is to inspect them. zonal_statistics cannot be used at this time.