5

I have this geopandas dataframe:

df1 = gpd.GeoDataFrame(df,crs=rasterio.crs.CRS.from_string("epsg:4326"), geometry=df['geojson'].apply(json_load))
df1 = df1.set_crs(epsg=4326)
print (df1.shape)
df1.head()

 account_id area    geojson                                             geometry
0   1365    1.7     {"type": "MultiPolygon","coordinates":[[[[35.4...   MULTIPOLYGON (((35.48632 32.78037, 35.48600 32...
1   3019    0.6     {"type": "MultiPolygon","coordinates":[[[[4.92...   MULTIPOLYGON (((4.92786 43.68305, 4.92951 43.6...



df1.crs
>> <Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

I would like to get the coordinates of the centroid of each polygon. I used to be able to get it by df1['geometry'].centroid without any warning, but now when I try that line I get:

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.

  """Entry point for launching an IPython kernel.

I don't understand why because my CRS looks fine, also the coordinates of the centroid look fine, should I just ignore this warning?

My packages versions:

print(pyproj.__version__)
print(rasterio.__version__)
print(gpd.__version__)

>> 2.6.1.post1
>> 1.1.5
>> 0.8.1
2

I am no export on geopandas, but to me the warning only tells you that you should use a crs with coordinates in meters,

Probably because the function calculating the centroid doesn't support ellipsoidal calculations or something.

So if your polygons are relatively small, wich they seem to be, the 'curvature' of the ellipsoid doesn't have much of an impact and treats your coordinates as if they are in meters and does cartesian calculations, hence why your centroid is 'correct'.

EDIT:

This answer to a similar question seems to point to my answer being right

4
  • So I need to project it, then ask for the centroid coordinates, then project both the coordinates of the centroid and the polygon back to the original CRS?
    – user88484
    Aug 27 '20 at 12:53
  • 1
    i believe you have two choices: - 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 Aug 27 '20 at 12:58
  • If I do that, I will use the UTM projection, and since I have many polygons that can be anywhere in the world I'll have also to find the correct UTM zone.. so I think I'll look for a solution outside GeoPandas. Thanks
    – user88484
    Aug 27 '20 at 13:02
  • pages 65 to 75 might help you if can implement the mathematical formulas: jennessent.com/downloads/Graphics_Shapes_Online.pdf ps: might probably be easier to look far and wide if somethink like that doesn't exist somewhere x) Aug 27 '20 at 13:20
2

Use df1 = df1.set_crs(epsg=3395) to project to World Mercator before calculating the centroid.

You can then re-project back as necessary.

2
  • Mercator seems like a really bad projection to use when calculating centroids, since meridonal distances at the poles are amplified, no? At the very least something like WGS84 Equidistance Cylindrical ('epsg:32663') would be preferable. Jun 17 at 19:30
  • actually looks like "epsg:4087" is the updated standard Jun 17 at 19:36
0

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)

Cylindrical, or PlateCarree projection, showing the default centroid appearing correct, below the Mercator centroid and above the equal area centroid

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)

The equal area cylindrical projection preserves area. When viewed on a globe, it becomes apparent that the equal area projection correctly identifies the centroid

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)

The Mercator projection expands the distance between latitudes close to the poles, inflating areas asymptotically close to either pole

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.

Your Answer

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

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