GeoPandas – How to Compute Centroids Based on Field While Keeping Aggregated Values

aggregatecentroidsgeopandaspythonshapely

I have a shapefile containing 85.000 points representing houses. I want to create a new shapefile containing the centroid of each subset of houses with the same postal code. This centroids file should also contain the mean values for attributes of the square meters and building year of the houses in that postal code. This last part is not working.

My workflow is as follows:

Dissolve the points shapefile into a multipoints file, where each multipoint is a group of points with the same postal code. Aggregate mean values for other attributes in the multipoints.
Compute the centroid of each multipoint
Write to new shapefile

This is my code for dissolving and computing centroids:

gpd = geopandas.read_file(my_file)
dissolved_gpd = gpd.dissolve(by='postcode',aggfunc='mean')
centroids = dissolved_gpd.centroid

When I print dissolved_gpd and centroids they look like this (just 1 line per example).
dissolved_gpd:

                   geometry                                  ...  bouwjaar
postcode                                                     ...             
6511AA    MULTIPOINT (187217.617 428676.815, 187254.576 ...  ...  1983.692308
[5048 rows x 4 columns]

centroids:

postcode
6511AA    POINT (187242.378 428870.156)
Length: 5048, dtype: geometry

As you can see, the 'bouwjaar' aggregated attribute is lost when computing the centroids. Is there any way to pass this on when calculating centroids in GeoPandas? I can't find anything about this in the docs. Or is there another method of computing centroids where it's possible to pass on aggregated values?

Best Answer

The way geopandas works is indeed a little confusing. When you call centroid it is returning a GeoSeries (just the geometries) and not the full GeoDataFrame (geometries + attributes).
Do this to retain the full GeoDataFrame and replace the MultiPoint geometries with their centroid.

gpd = geopandas.read_file(my_file)
dissolved_gpd = gpd.dissolve(by='postcode',aggfunc='mean')
dissolved_gpd['geometry'] = dissolved_gpd.geometry.centroid

Related Solutions

Geometry – Calculating Spherical Polygon Centroid

It won't work consistently even when you perform all triangulations relative to a single, fixed point. The problem is that spherical and Euclidean calculations are being mixed without any consideration of what they might mean.

One way to make this obvious is to consider a rather extreme triangle, such as almost one-half of a hemisphere. For instance, starting at (lon,lat) = (-179, 0), run along the equator to (0, 0), then up to the north pole at (0, 90), then back to the beginning at (-179, 0). This is a 90-179-90 triangle comprising most of the northern half of the western hemisphere. The problem is that its endpoints (shown as white dots in the figure) lie practically in a plane: one is at the pole and the other two are almost on opposite sides of it. Thus their average, projected back to the sphere (the red dot), is almost at the pole--but that's about as far from any reasonable center as one can get:

Large spherical triangle

As another example, let's triangulate a polygon representing the upper hemisphere relative to its center, the North Pole. We will always divide the Western hemisphere into two equal halves, each of them a 90-90-90 triangle (thereby avoiding any problems with huge, hemisphere-spanning triangles). The Eastern hemisphere, however, will be divided into n equal semi-lunes. The vertices of lune k (k = 1, 2, ..., n) have (lon, lat) coordinates

((k-1) * 180/n, 0),  (k * 180/n, 0),  (k * 180/n, 90).

Lunes for k=8

This figure shows the setup for k=8. The red dots are the individual triangle "centers" computed according to the "Tools for Graphics and Shapes" document, pp 65-67.

Doing the calculations, I find that with k = 2, the area-weighted center indeed is at the North Pole (as would be indicated by symmetry considerations), but as n increases, the result quickly shifts into the Western hemisphere and, in the limit, approaches a latitude of 89.556 degrees along longitude -90 degrees. This is approximately 50 kilometers south of the North Pole itself.

Admittedly, a +/-50 kilometer error for a polygon spanning 20,000 kilometers is small; the total amount of arbitrary variation due to different triangulations in this case is only 0.5%. Evidently the relative errors can be made arbitrarily large by including negative triangles (just add and subtract some really large triangles relative to a small triangle). Regardless, anyone going to the effort of doing spherical computations evidently is trying to avoid projection errors, so they are looking for high accuracy. This triangulation method cannot be recommended.

GeoPandas – Fixing AttributeError When Field Names Are Tuples

Convert field names to string

dissolved.columns = [field_name[0] + "_" + field_name[1] if type(field_name) is tuple else field_name for field_name in list(dissolved.columns.values)]

Reproducible example

import geopandas
from shapely.geometry import Point
d = {("col1", "mean"): ['name1', 'name2'], 'geometry': [Point(1, 2), Point(2, 1)]}
gdf = geopandas.GeoDataFrame(d, crs="EPSG:4326")
gdf.columns = [field_name[0] + "_" + field_name[1] if type(field_name) is tuple else field_name for field_name in list(gdf.columns.values)]

Best Answer

Related Solutions

Geometry – Calculating Spherical Polygon Centroid

GeoPandas – Fixing AttributeError When Field Names Are Tuples

Related Question