Round-trip the data through a flat projection, ideally one which preserves area, such as Equal Area Cylindrical ('+proj=cea'):
df1.to_crs('+proj=cea').centroid.to_crs(df1.crs)
This projects the shapes onto a flat surface, which can then be used to find the centroid, and then converts back into the original coordinate system.
Choosing a projection
Note that the choice of a projection is very important for large shapes, especially near the poles.
As an example - I'll grab Greenland from Natural Earth:
import matplotlib.pyplot as plt
import matplotlib
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import geopandas as gpd
import shapely
f = cfeature.NaturalEarthFeature('cultural', 'admin_0_countries_lakes', '50m')
countries = list(f.geometries())
filter with an arbitrary point in greenland
greenland, = [
c for c in countries
if c.contains(shapely.geometry.Point(-41.4408857, 68.7165341))
]
gdf = gpd.GeoDataFrame({'geometry': [greenland]}, crs='epsg:4326')
Depending on the projection chosen, you can get a wide range of centroids:
>>> default_centroid = gdf.centroid
<ipython-input-123-350c3321d4e0>:2: UserWarning: Geometry is in a
geographic CRS. Results from 'centroid' are likely incorrect. Use
'GeoSeries.to_crs()' to re-project geometries to a projected CRS
before this operation.
default_centroid = gdf.centroid
>>> default_centroid
0 POINT (-41.34191 74.71051)
dtype: geometry
>>> centroid_mercator = gdf.to_crs('epsg:3785').centroid.to_crs(gdf.crs)
>>> centroid_mercator
0 POINT (-41.08900 77.31320)
dtype: geometry
>>> centroid_equal_area = gdf.to_crs('+proj=cea').centroid.to_crs(gdf.crs)
>>> centroid_equal_area
0 POINT (-41.55424 71.98766)
dtype: geometry
These three points are hundreds of km away from one another! So what's going on here?
The issue is that the different projections distort distances near the poles in different ways. The default centroid relies on math based on simple cartesian coordinates. This is equivalent to projecting the shape onto a simple Cylindrical (PlateCarree) projection - you get the same answer as the default if you use a cylindrical projection (e.g. "epsg:4087"). This inflates distances between longitudes, causing the centroid to be too-heavily influenced by the expanded northern parts of the shape. The mercator projection is even worse, and dramatically over-emphasizes higher-latitude areas. To correct this, use an area-preserving projection, such as Equal Area Cylindrical ('+proj=cea').
Visualizing projection effects
PlateCarree (WGS84, or simple lat/lon cartesian projection)
This is the default projection, and is equivalent to using a simple Cylindrical projection. The red centroid dot appears in the center of the shape, because northern latitudes have been inflated.
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.PlateCarree()})
gdf.plot(ax=ax, alpha=0.6, transform=ccrs.PlateCarree())
ax.set_extent((-120, 40, 55, 85), crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='k', alpha=0.5, linestyle='dotted')
l1 = ax.scatter(default_centroid.x, default_centroid.y, color='red', transform=ccrs.PlateCarree())
l2 = ax.scatter(centroid_mercator.x, centroid_mercator.y, color='blue', transform=ccrs.PlateCarree())
l3 = ax.scatter(centroid_equal_area.x, centroid_equal_area.y, color='k', transform=ccrs.PlateCarree())
ax.legend([l1, l2, l3], ['default', 'mercator', 'equal area'], bbox_to_anchor=(0.5, -0.2), loc='lower center', ncol=3)
Equal area cylindrical projection
We can view the placement of the equal-area-based centroid by projecting the shape on a globe using the Orthographic projection. This allows us to see that using an equal area projection offers the truest centroid.
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.Orthographic(centroid_equal_area.x[0], centroid_equal_area.y[0])})
gdf.plot(ax=ax, alpha=0.6, transform=ccrs.PlateCarree())
ax.set_extent((-120, 40, 55, 85), crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='k', alpha=0.5, linestyle='dotted')
l1 = ax.scatter(default_centroid.x, default_centroid.y, color='red', transform=ccrs.PlateCarree())
l2 = ax.scatter(centroid_mercator.x, centroid_mercator.y, color='blue', transform=ccrs.PlateCarree())
l3 = ax.scatter(centroid_equal_area.x, centroid_equal_area.y, color='k', transform=ccrs.PlateCarree())
ax.legend([l1, l2, l3], ['default', 'mercator', 'equal area'], bbox_to_anchor=(0.5, -0.2), loc='lower center', ncol=3)
Mercator
The mercator projection amplifies the distortion near the poles, giving a dramatic bias toward northern latitudes:
fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.Mercator()})
gdf.plot(ax=ax, alpha=0.6, transform=ccrs.PlateCarree())
ax.set_extent((-120, 40, 55, 85), crs=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=1, color='k', alpha=0.5, linestyle='dotted')
l1 = ax.scatter(default_centroid.x, default_centroid.y, color='red', transform=ccrs.PlateCarree())
l2 = ax.scatter(centroid_mercator.x, centroid_mercator.y, color='blue', transform=ccrs.PlateCarree())
l3 = ax.scatter(centroid_equal_area.x, centroid_equal_area.y, color='k', transform=ccrs.PlateCarree())
ax.legend([l1, l2, l3], ['default', 'mercator', 'equal area'], bbox_to_anchor=(0.5, -0.2), loc='lower center', ncol=3)
Note that most web maps use a version of the Mercator projection. Therefore, if the goal is to have a visual centroid, Mercator may be the correct choice.
- you can do all your calculations in a projected CRS and only return to wgs84 at the end
- or project your polygons to projected CRS, get your centroids, reproject you centroids and polygons to wgs84, but if you need to do something else with either you might run into the same problem
– Kalak Aug 27 '20 at 12:58GeoPandas. Thanks – user88484 Aug 27 '20 at 13:02ps: might probably be easier to look far and wide if somethink like that doesn't exist somewhere x)
– Kalak Aug 27 '20 at 13:20