3

I have a series of many GeoPandas Polygon objects, each with an associated exposure time. My goal is to find the union of all the Polygons in order to create an exposure time map.

For each Polygon in the final exposure map, I want to have the total exposure time, which is defined as the sum of the exposure times of each overlapping Polygon.

I think I have it working for the case of 2 geodataframes:

import geopandas as gpd
import numpy as np
from rtree import index
from shapely.geometry import Polygon
from shapely.ops import unary_union


coordinates1a = [(1,1), (1, 2), (2, 2), (2, 1)]
coordinates1b = [(2.1,2.1), (2.1, 3.1), (3.1, 3.1), (3.1, 2.1)]
coordinates2a = [(0.5,0.5), (0.5, 1.5), (1.5, 1.5), (1.5, 0.5)]
coordinates2b = [(1.6,1.6), (1.6, 2.6), (2.6, 2.6), (2.6, 1.6)]

# Create a Shapely polygon from the coordinate-tuple lists
poly1a = Polygon(coordinates1a)
poly1b = Polygon(coordinates1b)
poly2a = Polygon(coordinates2a)
poly2b = Polygon(coordinates2b)

polys1 = gpd.GeoSeries([poly1a, poly1b])
polys2 = gpd.GeoSeries([poly2a, poly2b])

# Create geodataframes
df1 = gpd.GeoDataFrame({'geometry': polys1, 'df1':[1,2], 'exposure_time': 10.7})
df2 = gpd.GeoDataFrame({'geometry': polys2, 'df2':[1,2], 'exposure_time': 9.7})

res_union = gpd.overlay(df1, df2, how='union')

# Find the list of all exposure time column names
exp_cols = [c for c in res_union.columns if 'exposure_time' in c]

# Sum the contributing exposure times in each polygon of the union
for i in range(len(res_union)):
    exptime = np.sum(res_union.loc[i, exp_cols])
    all_exptimes.append(exptime)
    res_union.loc[i, 'total_exposure_time'] = exptime
    res_union.loc[i, 'area'] = res_union.loc[i, 'geometry'].area

How do I do the same for the case of N geodataframes, where N might be tens or hundreds?

I tried using an rtree index and unary_union(), but couldn't quite get it to work, plus that only seemed to work on the Polygons themselves rather than the geodataframes, meaning that the exposure time information was lost. I am a complete novice with GeoPandas, so maybe I'm missing something obvious, but after a lot of searching and reading documentation yesterday, I'm stumped.

# Try using rtree
df_array = [df1, df2]
polygon_array = [poly1a, poly1b, poly2a, poly2b]

# Need to populate the index with Polygons, rather than geodataframes,
# so how do we keep track of the exposure time information?
intersections = []
idx = index.Index()
pos = 0
for polygons in df_array:
    for index in range(len(polygons)):
        idx.insert(pos, polygons.loc[index, 'geometry'].bounds)
        pos += 1

# This seems like the more natural way to do it, but exposure time
# info is not present
for polygon in polygon_array:
    merged_polygons = unary_union([polygon_array[pos] for pos in idx.intersection(polygon.bounds) if polygon_array[pos] != polygon])
    intersections.append(polygon.intersection(merged_polygons))

intersection = unary_union(intersections)

Here's an image of a case closer to what I eventually want to have. My goal is to have a total exposure time associated with each polygon. In this image, I'm showing the union of three geodataframes, where each geodataframe contains 4 non-overlapping squares. The squares' locations relative to one another are fixed across all geodataframes, but all 4 squares can shift as a unit, from geodataframe to geodataframe.

Union of three geodataframes, where each geodataframe contains 4 non-overlapping squares. The squares' locations relative to one another are fixed, but all 4 squares can shift as a unit, from geodataframe to geodataframe.

Your Answer

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

Browse other questions tagged or ask your own question.