# Calculate minimum area where x% of raster value total is concentrated within polygons

I have a raster layer (chlorophyll concentration) and a layer of non-overlapping polygons (regions of the ocean). For each polygon, I'd like to calculate the area where the majority of chlorophyll is concentrated – e.g., the smallest area value covered by 80% of the total chlorophyll in that polygon. I cannot think of an elegant way of achieving this, however. I have access to both ArcGIS (10.1) and QGIS (2.8 Wien), so a solution with either is fine.

One crude approach is to start at a high threshold value and calculate the amount of chlorophyll above that threshold for each polygon, then turn that into a percent of the total chlorophyll for that polygon. Incrementally, I decrease the threshold until I reach 80% for each polygon and calculate the corresponding area. This is horrifically cludgey and I know that there must be a more efficient solution. Is there a better way of doing this is ArcGIS or QGIS?

• How would you like to handle areas where chlorophyll conc. is constant throughout the polygon? What does the final product look like to you: raster or vector? Is the "smallest area covered by 80% of the total chlorophyll in that polygon" the final product? – Aaron May 16 '15 at 1:02
• I think that if your "crude" approach will give you the correct result, then you should go for it. I'm assuming that you don't need the result to be contiguous within the polygons? Do you do any python? – user1269942 May 16 '15 at 6:12
• @Aaron Cases where chlorophyll is constant (although not present in my particular polygons) should return the area of the entire polygon. The final product could be a column in the attribute table for the polygon layer. Yes, the smallest area covered by 80% of the total chlorophyll in the final product. – Lyngbakr May 18 '15 at 16:40
• @user1269942 The crude approach is very labour intensive, which is why I'm looking for a more elegant solution. I essentially have to use an approach like the bisection method to hone in on 80% for each polygon and I have >150 polygons. Unfortunately, I don't know Python. – Lyngbakr May 18 '15 at 16:44

Here’s my solution. I gave up with ArcGIS and QGIS, instead opting to output the data to file and perform the calculation in R.

# Load library
library(raster)

chloro <- raster("../data/chloro.tif")

poly <- shapefile("../data/poly.shp")

# Extract data in polygons
poly.chloro <- extract(chloro, poly)

# Determine proportion of polygon where there is specified proportion of data
poly.prop <- function(x, threshold){
foo <- approxfun(1 - cumsum(sort(x, na.last = NA))/sum(x, na.rm=TRUE), sort(x, na.last = NA))
return(sum(x - foo(threshold)>0, na.rm=TRUE)/length(x))
}

# Null array to hold threshold values
threshold.values <- NULL

# Loop through polygons
for(i in poly.chloro) threshold.values <- c(threshold.values, poly.prop(i, 0.8))


The above code determines the proportion of each polygon that is occupied by 80% of the total chlorophyll. Multiplying these proportions by the polygons’ area gives the area occupied by 80% of the chlorophyll, as required.

Initially, I load the raster library, as well as the chlorophyll raster and polygon shape files. Next, I extract all of the raster points that are located in each polygon. The poly.prop function performs the actual calculation. The raster values associated with each polygon are passed to this function, which sorts them into ascending order. A cumulative sum of this sorted vector is calculated, which is then divided by the sum total of the vector. This gives the proportion of the total chlorophyll below each point in the sorted vector. For example, say at the nth element of the vector has the value 0.2. This means that 20% of the chlorophyll is below that point. However, I’m phrasing my question in terms of how much is above that point, so I subtract that value from 1. (In my example, 1 - 0.2 = 0.8 means 80% is above that point.) Next, I use the approxfun function to derive an empirical relationship between the values in the sorted vector (i.e., concentrations of chlorophyll) and the corresponding proportion of total chlorophyll above that point. This function takes a value, like 0.8, and returns the concentration of chlorophyll in the data set where there is 80% of the total chlorophyll above it. Finally, I subtract this threshold concentration from all of the raster values in that polygon and count how many raster values are above zero (i.e., above the threshold). When this is divided by the total number of raster values in that polygon, it gives the proportion of the polygon that is occupied by these points.

I'm not sure if that explanation is clear, but it works fine – I checked it against my initial rough method.

PS Thanks to everyone for your solutions and comments. This solution draws on several on the ideas presented by others above.

I would convert the polygon layer into a raster data set with feature to raster tool and then do a raster overlay analysis using raster calculator. Make sure you are using the same projection and coordinate systems for both or your results will be wrong. Good luck!

• What formula/math would you perform in the raster calculator to solve the specific question of determining which pixels of a zone held 80% of the total zone value in the smallest number of adjacent/contiguous pixels? It seems like you have a start to an answer here, but without further detail I'm not sure how it addresses the issue. – Chris W May 16 '15 at 1:00

There are a couple of ways of tackling this problem. The first is to try and find a tool to do the job. The second is to use a bit of Python scripting to do it manually. Since finding tools seems to be hard, here's a python solution.

You'll need to figure out how to install python/gdal/numpy yourself, but once you do you won't be disappointed. It's easiest on linux but there are plenty of 'how tos' out there.

I disagree that the 'brute force' approach is too complicated. With only 150 polygons, you likely will not even be able to take a sip of your coffee before this finishes running...unless there are millions of pixels per polygon, in which case you'll get a few sips in.

Here's a snippet I put together...not tested but I looked it over and I'm pretty sure it'll work for you.

[edit: print out area as it goes] [another edit: gdal ImportError in python on Windows may be helpful for getting python+gdal running]

import gdal
import numpy as np

area_per_pixel = 100 #??? you fill in

#chlorophyll_raster  is a numpy array made from using gdal to open the raster
r = gdal.Open('somefilename')

#polygon_raster is a raster with the exact same shape/size/resolution as chlorophyll_raster, with values equal to the polygon-id
#this gdal_rasterize commandline utility may be useful here.
r2 = gdal.Open('somefilename2')

#max/min used in thresholding
max_c = np.max(chlorophyll_raster)
min_c = np.min(chlorophyll_raster)

#make blank raster with same shape as others...nodata default
#in the end, this raster will have values equal to "1" in pixel
#locations that contribute to 80% per polygon
80_pc_raster = np.ones(chlorophyll_raster.shape) * -9999

#tune this(0.01) to match the scale of your chorophyll concentration
increments = (max_c - min_c) / .01

for pid in np.unique(polygon_raster):#for each polygon
for threshold in np.linspace(max_c, min_c, increments):#go through each level (brute force). start with max
total_c_for_polygon = np.sum(chlorophyll_raster[polygon_raster == pid])
threshold_mask = (polygon_raster == pid) & (chlorophyll_raster >= threshold)
pc_above_threshold = amount_c_above_threshold_for_polygon / total_c_for_polygon
if pc_above_threshold > 0.8: #if your increment is small, as soon as you cross the 80%
#fill in the 80_pc_raster raster with the pixels that comprise the 80%
#print polygon-id, area
print pid, ', ', np.sum(threshold_mask) * area_per_pixel
break #break out of inner loop and on to the next polygon

new_raster = gdal.GetDriverByName('GTiff').Create('80pc_raster.tif', 80_pc_raster.shape, 80_pc_raster.shape, 1, gdal.GDT_Float32)
#new_raster.SetGeoTransform(geo_transform) #should do this
#new_raster.SetProjection(projection)  #this too
new_raster.GetRasterBand(1).SetNoDataValue(-9999)
new_raster.GetRasterBand(1).WriteArray(80_pc_raster)

#all done.


good luck!

• @lyngbakr I've been trying to think of a tool or combination of tools that would do this and can't get any farther than Zonal stats and a field calc to come up with the 80% values for each poly. From there I know that, as user1269942's answer suggests, you start with the max cell and work your way down the value list, adding until you hit the 80% mark - preferably within some fudge factor in case a group of cells with the same value pushes you over the edge. And I don't know an existing tool/method to do that. If this script works, it's got my vote. – Chris W May 18 '15 at 23:32

Could you calculate quintiles for the raster within each polygon and use these values to create contour lines where the lower quintile represents the 80/20 boundary?

• How would I go about doing that? – Lyngbakr May 18 '15 at 16:46
• I'm not familiar with ArcGIS so can't comment, and I've been trying to find a tidy way to do it in QGIS and so far haven't found one. – Adrian May 18 '15 at 17:05
• I think your general approach is right, though. Perhaps the following workflow: 1) find the 20th percentile for each polygon, 2) mask the data below that value, 3) sum the area covered by the remaining data. – Lyngbakr May 18 '15 at 17:10
• I'm not familiar with ArcGIS so can't comment but thinking about how it might be possible in QGIS (I've not had time to try it), I wonder if it may be possible to use the Processing toolbox in QGIS to create a model using an R script to calculate the 80/20 quartile boundary with the raster library eg. quantile(myraster, probs = c(0.20), type=7,names = FALSE). This would give you the 80/20 boundary as an input to the QGIS Raster/Extract/Contours tool to generate a contour line around the area. Equally you could use the raster calculator to set any value less than the generated value to 0. – Adrian May 18 '15 at 17:17
• I just found a plug-in for QGIS called "RasterPixelCountStat" that will calculate the breaks for you - appears that's all it does so you'd still need create a Processing model or similar to use the value it generates – Adrian May 18 '15 at 17:30