I have a global coastline grid made up of polygons (hexagons) and I am trying to aggregate global data onto these hexagons. For example, I am trying to compute the distance of each hexagon to the nearest coral reef (also defined as a polygon). Since I use global datasets, which are quite large, I am using GeoPandas sjoin_nearest() function (both datasets are geodataframes), which links each hexagon to the nearest coral polygon (typically the distance will be in not be more than a few hundred km).
Do the distances that GeoPandas uses for this actually make sense? The CRS I am using for both datasets is EPSG:6933, which is an ellipsoidal, equal area projection with unit meters. However, the documentation of sjoin_nearest states that "Since this join relies on distances, results will be inaccurate if your geometries are in a geographic CRS. Every operation in GeoPandas is planar, i.e. the potential third dimension is not taken into account."
Does this mean I cannot use sjoin_nearest at all for geographic calculations, even with an area-preserving projection? And if not, does someone know a better way to do this?
I know that I will never get perfectly accurate distances, but I only need an accuracy in the order of 10km. Also, the analysis is limited to between 47.5 and -37.5 degrees latitude, so relatively far from the poles.