1

I have a DataFrame that resulted from spatially joining a digital elevation map with a square grid map.

enter image description here

This unexpectedly resulted in duplicate rows where two rows will have the same "geometry" but a different "Elevation" value.

How do I get the median of "Elevation" for each unique "geometry"? I'm new to GeoPandas, so I tried the traditional methods of aggregating a DataFrame, but found that "geometry" cannot be operated with the groupby() function.

mrkna_grid.groupby("geometry")['Elevation'].median()

I have also tried using the dissolve() function, but I don't think I'm doing it correctly because the number of rows were reduced to just seventy (70) as opposed to the original two thousand (2000) before the spatial join.

mrkna_grid.dissolve(by="Elevation", aggfunc="median")
2
  • 1
    Do you want to keep the duplicate geometries?
    – BERA
    Nov 5, 2021 at 8:20
  • @BERA no, I just need every unique "geometry" and the median of all its "Elevation" associated with said "geometry"
    – SS-Salt
    Nov 5, 2021 at 8:30

1 Answer 1

1

Let's assume there is a polygon layer 'layer' with its attribute table, see the image below.

input

The expected result of median is:

Polygon 1 | 4
Polygon 2 | 7
Polygon 3 | 15

Using one of the following code:

Solution #1 includes Pandas/GeoPandas native DataFrame.median()

import geopandas as gpd

path_to_layer = "C:/Documents/Python Scripts/median/layer.shp"
layer = gpd.read_file(path_to_layer)

layer_ = layer.groupby(layer["geometry"].to_wkt())['Elevation'].median().reset_index(name='MedianElev')
layer_['geometry'] = gpd.GeoSeries.from_wkt(layer_['index'])
layer_.drop(['index'], inplace=True, axis=1)

output = gpd.GeoDataFrame(layer_, geometry='geometry')
output = output.set_crs(layer.crs)
output.to_file("output.shp")

Solution #2 includes Python's statistics.median()

import geopandas as gpd
from statistics import median

path_to_layer = "C:/Documents/Python Scripts/median/layer.shp"
layer = gpd.read_file(path_to_layer)

layer['geom_wkt'] = layer['geometry'].to_wkt()

layer_ = layer.groupby('geom_wkt')['Elevation'].apply(list).reset_index(name='ElevationList')
layer_['median'] = layer_['ElevationList'].apply(lambda x: median(x))
layer_['geometry'] = gpd.GeoSeries.from_wkt(layer_['geom_wkt'])
layer_.drop(['geom_wkt', 'ElevationList'], inplace=True, axis=1)

output = gpd.GeoDataFrame(layer_, geometry='geometry')
output = output.set_crs(layer.crs)
output.to_file("output.shp")

it is possible to achieve the following output:

result


References:

3
  • Thank you. One question though: the code above returns a dataframe that only has "geometry" and "Elevation". How would you hypothetically retain other attributes?
    – SS-Salt
    Nov 5, 2021 at 13:00
  • I will probably do the join on e.g. the "geom_wkt" attribute. What kind of attributes would you like to retain? Can you share a piece of your data?
    – Taras
    Nov 5, 2021 at 13:02
  • The full picture of the DataFrame in the OP has attributes X, Y and Z (for brevity) alongside the shown "geometry" and "Elevation".
    – SS-Salt
    Nov 5, 2021 at 13:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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