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

Question

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?

Answer 1

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)

# Load chlorophyll data
chloro <- raster("../data/chloro.tif")

# Load polygons
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.

Answer 2

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!

Answer 3

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')
chlorophyll_raster = np.array(r.GetRasterBand(1).ReadAsArray())

#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')
polygon_raster = np.array(r2.GetRasterBand(1).ReadAsArray())

#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)
    amount_c_above_threshold_for_polygon = np.sum(chlorophyll_raster[threshold_mask])
    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%
      80_pc_raster[threshold_mask] = 1
      #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[1], 80_pc_raster.shape[0], 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!

Answer 4

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?