7

I found this thread in the old ESRI forums with the same type of problem I am facing today. The solution to the problem is explained in little detail and since I am not an expert, I was wondering if somebody could describe it to me in little more detail.

We are doing a mental map analysis. We have let's say 50 polygons of a certain city square which are all different in shape, size and location. Our goal is to create a single polygon which would represent the answers from our whole input data. So basically, we need a method for generating a single polygon out of these 50 which would be an "average" polygon in size, shape and location.

Can somebody explain the process?

  • 5
    I can't see the thread, only the error page "Thread id is not set" – nadya Aug 11 '13 at 3:02
  • 1
  • The first ESRI forum thread @mkennedy found looks appropriate but it describes two distinct solutions corresponding to radically different interpretations of "average polygon". The second thread offers yet a third interpretation. Which interpretation is intended here? – whuber Aug 12 '13 at 16:29
  • 1
    @mkennedy and nadya - yes, it was the first one from mkennedy post, when admin edited my post it got screwed up @ whuber - The "average polygon" as intended here would have average area of all these 50 (easy to calculate), would be in a location similar to centroid, let's call it "a centroid of all the 50 centroids" and would have a shape which would be as close as possible to represent all input polygons. I hope this clarifies it at least a little bit. If you need any more info please ask. – slowhand Aug 12 '13 at 17:39
  • 1
    I think that unless you had a very smart algorithm, you'll end up with a very rounded-looking polygon, with an "accurate" area but nondescript shape. Then there's the issue of weighting -- are some inputs of higher quality than others? Interesting question, all told. – Martin F Aug 15 '13 at 20:25
11

What if you used a grid approach? Convert each of the polygons to a raster, assign each cell a value of 1, then add all of the rasters together. The resulting polygon would be formed by the highest-value cells whose area equals the average area of the input polygons.

The cell size would have to be small enough to make the difference between the different polygons meaningful.

  • 1
    +1 That's a clever and simple solution. (In the last sentence I think you mean "cell size" instead of "grid size"; the latter could be construed in its usual sense of "grid extent," which would make your statement rather puzzling!) – whuber Aug 14 '13 at 19:23
  • Thanks for pointing that out. The appropriate edits have been made. – Bjorn Aug 19 '13 at 14:32
12

I've figured out an algorithm for the grid approach using several Python tools. Rasterising and polygonising is done with rasterio, which is based on GDAL/OGR. Here are most of the imports:

import rasterio
import numpy as np
from rasterio import Affine, features
from shapely.geometry import mapping, shape
from shapely.ops import cascaded_union
from math import floor, ceil, sqrt

First, get a few polygon shapes to do statistics with. The list of shape must be sort-of overlapping, otherwise this procedure won't work as well. There are many way to get these shape into Python, e.g., read a Shapefile with fiona.

import fiona
shp_fname = 'my_overlapping_polygons.shp'
shapes = []
with fiona.open(shp_fname) as ds:
    minx, miny, maxx, maxy = ds.bounds
    for f in ds:
        shapes.append(f['geometry'])

Or how about these five polygons:

shapes = [
    {'type': 'Polygon', 'coordinates': [[(1095.76, 278.06), (1095.76, 278.06), (1228.25, 301.98), (1377.29, 301.98), (1511.62, 283.58), (1603.62, 254.14), (1669.86, 224.7), (1737.95, 175.02), (1772.91, 129.01), (1791.31, 77.49), (1804.19, -1.63), (1796.83, -53.15), (1776.59, -121.24), (1726.91, -198.52), (1629.38, -303.4), (1491.38, -413.81), (1215.37, -575.73), (764.55, -809.42), (617.34, -883.03), (508.78, -929.03), (431.5, -951.11), (210.69, -965.83), (135.24, -938.23), (111.32, -888.55), (96.6, -783.66), (126.04, -619.9), (194.13, -469.01), (295.33, -296.04), (381.81, -150.68), (501.42, -20.03), (630.22, 83.01), (771.91, 167.66), (924.63, 232.06), (1027.68, 261.5), (1095.76, 278.06)]]},
    {'type': 'Polygon', 'coordinates': [[(1865.28, 145.78), (1865.28, 145.78), (1779.15, 286.31), (1629.55, 381.5), (1425.57, 438.17), (1226.11, 435.9), (1037.99, 404.17), (829.46, 306.71), (657.21, 170.72), (548.41, 32.46), (466.82, -87.67), (328.56, -407.25), (287.76, -559.11), (287.76, -731.37), (321.76, -869.63), (385.22, -944.42), (480.42, -967.09), (729.74, -971.62), (913.33, -917.23), (1144.51, -806.17), (1432.37, -647.51), (1659.02, -482.05), (1819.94, -302.99), (1908.34, -117.14), (1901.54, 14.32), (1865.28, 145.78)]]},
    {'type': 'Polygon', 'coordinates': [[(1175.76, 247.32), (1175.76, 247.32), (1336.5, 258.21), (1450.92, 251.4), (1550.36, 229.61), (1645.71, 195.55), (1724.72, 150.6), (1758.78, 111.1), (1777.85, -19.67), (1765.59, -71.44), (1709.74, -157.25), (1603.49, -258.06), (1463.18, -362.95), (1181.21, -504.61), (524.63, -841.08), (305.32, -965.04), (211.33, -1007.26), (-21.61, -1049.49), (-82.91, -1034.51), (-111.51, -975.93), (-111.51, -857.42), (-86.99, -745.72), (50.59, -505.98), (143.22, -332.98), (290.33, -165.43), (470.14, -30.57), (659.49, 78.41), (881.52, 175.12), (1044.99, 224.16), (1175.76, 247.32)]]},
    {'type': 'Polygon', 'coordinates': [[(886.58, 201.11), (886.58, 201.11), (1106.77, 271.57), (1249.89, 286.98), (1430.44, 286.98), (1531.73, 267.16), (1694.67, 205.51), (1760.72, 152.67), (1789.35, 106.43), (1798.15, 33.77), (1767.33, -107.15), (1613.2, -292.11), (1386.41, -450.64), (1150.81, -569.54), (710.44, -802.94), (441.81, -961.47), (325.11, -1020.92), (223.83, -1045.14), (49.88, -1067.16), (-16.18, -1047.35), (-27.19, -992.3), (-38.2, -913.03), (-16.18, -805.14), (32.26, -655.42), (175.38, -408.81), (340.52, -148.99), (494.65, -10.27), (688.42, 117.44), (813.92, 176.89), (886.58, 201.11)]]},
    {'type': 'Polygon', 'coordinates': [[(802.94, 60.03), (802.94, 60.03), (1012.93, 172.53), (1195.43, 230.02), (1370.42, 257.52), (1510.41, 250.02), (1610.41, 227.52), (1697.91, 195.02), (1755.41, 147.53), (1785.4, 102.53), (1795.4, 32.53), (1800.4, -57.47), (1790.4, -119.96), (1720.41, -227.46), (1585.41, -354.95), (1312.92, -552.45), (1055.43, -707.44), (730.45, -899.93), (540.45, -1009.93), (400.46, -1034.93), (275.46, -1044.93), (225.47, -1024.93), (197.97, -939.93), (200.47, -817.43), (272.96, -632.44), (367.96, -424.95), (472.96, -244.96), (612.95, -84.96), (752.94, 22.53), (802.94, 60.03)]]},
]

rings

In order to grid the vectors, some appropriate raster resolutions need to be determined, and applied to the data extents. E.g., I'm doing everything on a 1 metre grid, but it can be modified to any number. Then build an affine geotransform from the upper-left corner.

max_shape = cascaded_union([shape(s) for s in shapes])
minx, miny, maxx, maxy = max_shape.bounds
dx = dy = 1.0  # grid resolution; this can be adjusted
lenx = dx * (ceil(maxx / dx) - floor(minx / dx))
leny = dy * (ceil(maxy / dy) - floor(miny / dy))
assert lenx % dx == 0.0
assert leny % dy == 0.0
nx = int(lenx / dx)
ny = int(leny / dy)
gt = Affine(
    dx, 0.0, dx * floor(minx / dx),
    0.0, -dy, dy * ceil(maxy / dy))

Loop through each polygon, and rasterise them to a grid, where each will be 0 outside the polygon and 1 inside. With each result, accumulate a values to the array pa, which we'll normalise to values 0.0 (zero coverage) to 1.0 (all covered). This raster is basically the probability that a polygon covers the grid.

pa = np.zeros((ny, nx), 'd')
for s in shapes:
    r = features.rasterize([s], (ny, nx), transform=gt)
    pa[r > 0] += 1
pa /= len(shapes)  # normalise values

pa

The next step is really only needed if you want a smoothish result from only a few overlapping polygons. It will blur the sharp edges of probability array. You will need a modified gaussian_blur function with full convolution to avoid edge distortions, which modifies the geotransform size.

from scipy.signal import fftconvolve

def gaussian_blur(in_array, gt, size):
    """Gaussian blur, returns tuple `(ar, gt2)` that have been expanded by `size`"""
    # expand in_array to fit edge of kernel; constant value is zero
    padded_array = np.pad(in_array, size, 'constant')
    # build kernel
    x, y = np.mgrid[-size:size + 1, -size:size + 1]
    g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
    g = (g / g.sum()).astype(in_array.dtype)
    # do the Gaussian blur
    ar = fftconvolve(padded_array, g, mode='full')
    # convolved increased size of array ('full' option); update geotransform
    gt2 = Affine(
        gt.a, gt.b, gt.xoff - (2 * size * gt.a),
        gt.d, gt.e, gt.yoff - (2 * size * gt.e))
    return ar, gt2

Now use the function to do a Gaussian blur on a radius of 100 grid cells (or 100 m in this example), and return the larger-sized array and geotransform results:

spa, sgt = gaussian_blur(pa, gt, 100)

spa

Now the "average" result can be extracted by selecting a quantile threshold. So, using 0.1 will select a larger area that 10% or more of the polygons cover, while 0.95 will choose a smaller area with 95% coverage. Using 0.5 is the median quantile threshold close to what could be called "average". Because this example has an odd-number of shapes we get a nice result, but with even-numbers of samples the result is really sensitive.

thresh = 0.5  # median
pm = np.zeros(spa.shape, 'B')
pm[spa > thresh] = 1

pm

Convert the grid result back to vector, and simplify-away the jagged pixel shapes using a minimum distance of a grid diagonal (or more).

poly_shapes = []
for sh, val in features.shapes(pm, transform=sgt):
    if val == 1:
        poly_shapes.append(shape(sh))
if not any(poly_shapes):
    raise ValueError("could not find any shapes")
avg_poly = cascaded_union(poly_shapes)
# Simplify the polygon
simp_poly = avg_poly.simplify(sqrt(dx**2 + dy**2))
simp_shape = mapping(simp_poly)

result

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.