pythongeopandasshapely

Polygons elongation ratio calculation with geopandas - problem with creating a geodataframe with results


I have a polygon shapefile (actually, there are 170 such files) and want to calculate the elongation ratio for each polygon. The elongation ratio is the ratio of the short to the long side of the minimum rotated rectangle around the polygon. I have the following script which works well up to the point I want to create a geodataframe from the rectangles calculated in the for loop (the line causing the error is newgdf = gpd.GeoDataFrame(newData)) where I get this error: GeometryTypeError: Unknown geometry type: 'featurecollection'. It seems that each rectangle is a geoseries with one element. Maybe they should be concatenated somehow into a single series before being passed to GeoDataFrame()? How to create a geodataframe from mrrs?

Further explanations; why use a for loop? Some polygons have corrupted geometry (I got errors when calculating the minimum rotated rectangle in QGIS). I'm afraid that if applying minimum_rotated_rectangle() to the full geodataframe the bad polygons may be skipped. The resulting elongation table will then be shorter than the original shapefile and with shifted row IDs. I will not be able to join correctly the table with the elongations with the original shapefile. That is why I decided to use a for loop and try one polygon at a time to keep track of the id.

Any suggestions on how this can be done more rationally will be also appreciated.

The script:

import geopandas as gpd
from shapely.geometry import Polygon, LineString

gdf = gpd.read_file('original polygons.shp')

# Add field with unique id
gdf["shapeid"] = gdf.index
shapes = list(gdf.shapeid.unique())

mrrs = []
ids = []
elongs = []
for shape in shapes:
    
    gdfSubset = gdf[gdf['shapeid'] == shape]
    try:
        # Create minimum_rotated_rectangle
        mrr = gdfSubset.minimum_rotated_rectangle()
       
        coords = mrr.geometry.get_coordinates()
        x = coords['x'].tolist()
        y = coords['y'].tolist()
        poly = Polygon(list(zip(x, y)))
        
        # get the edges of a polygon as LineStrings based on
        # https://stackoverflow.com/a/68996466/5423259
        b = poly.boundary.coords
        linestrings = [LineString(b[k:k+2]) for k in range(len(b) - 1)]
        
        lengths = [int(ls.length) for ls in linestrings]
        # calculate elongation ratio (E) as:
        # E = 1 - S / L, where S is the short-axis length, and L is the long-axis length
        elong = 1 - (min(lengths) / max(lengths))
        elongs.append(elong)
        mrrs.append(mrr)
        ids.append(shape)
    except:
        pass

newData = {'shapeid': ids, 'elong': elongs, 'geometry': mrrs}
newgdf = gpd.GeoDataFrame(newData)
# export rectangles for checking purposes
newgdf.to_file('mrrs.shp')
# skip geometries and join the elongation field to the original shapes
newgdf = newgdf.drop('geometry', axis=1)
finalgdf = gdf.merge(newgdf, left_on='shapeid', right_on='shapeid', how='left')
# export original shapes with elongation field joined
finalgdf.to_file('final.shp')

Solution

  • A quick fix is to squeeze the geometry :

            elongs.append(elong)
            mrrs.append(mrr)               # (-)
            mrrs.append(mrr.squeeze())     # (+)
            ids.append(shape)
    

    Note that you can streamline the whole process (without the risk of omitting any polygons) :

    from itertools import pairwise
    
    def elong_ratio(poly):
        bnd = poly.boundary.coords
        lengths = list(map(lambda coo: LineString(coo).length, pairwise(bnd)))
        return 1 - min(lengths) / max(lengths)
    
    gdf["elong"] = gdf.make_valid().minimum_rotated_rectangle().apply(elong_ratio)