CRS Warning for k_ring_smoothing when using k param

[johnziebro]

What CRS should the GeoDataFrame being worked on have? I couldn't locate this in the docs or provided example notebooks.

When working with a GeoDataFrame with CRS 4326, a warning involving CRS is produced if k_ring_smoothing is applied using k instead of weight coefficients. The warning does not appear if weight coefficients are used.

test_smooth = test.h3.k_ring_smoothing(k=2)

FutureWarning: CRS mismatch between CRS of the passed geometries and 'crs'. Use 'GeoDataFrame.set_crs(crs, allow_override=True)' to overwrite CRS or 'GeoDataFrame.to_crs(crs)' to reproject geometries. CRS mismatch will raise an error in the future versions of GeoPandas.
lambda x: gpd.GeoDataFrame(x, crs="epsg:4326"),

test.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 ensemble
• Ellipsoid: WGS 84
• Prime Meridian: Greenwich

type(test)

geopandas.geodataframe.GeoDataFrame

test = test.h3.h3_is_valid()
test[test['h3_is_valid'] == False]

No rows show as invalid

[DahnJ]

Hi @johnziebro thanks for the question and report!

Coordinates are assumed to be in epsg=4326, in line with the core library. I have now made that more explicit in a2fdc39 which will later make it into the main branch.

I am not yet sure how the warning you posted could happen and did not manage to reproduce it. Could you please share what version of h3pandas and geopandas are you using? You can use e.g. h3pandas.__version__ to obtain the version number.

It would help a lot if you could also provide a small reproducible example for me to test.

[johnziebro]

Great module by the way!

h3pandas.__version__

'0.2.2'

gpd.__version__

'0.10.2'

import sys
print("Python version")
print (sys.version)
print("Version info.")
print (sys.version_info)

Python version
3.8.12 (default, Oct 10 2021, 10:19:55)
[GCC 10.3.0]
Version info.
sys.version_info(major=3, minor=8, micro=12, releaselevel='final', serial=0)

Jupyter Lab Kernel (it also occurs under Xeus kernel so not kernel specific):

{
"argv": [
"python",
"-m",
"ipykernel_launcher",
"-f",
"{connection_file}"
],
"display_name": "Python 3 (ipykernel)",
"language": "python",
"metadata": {
"debugger": true
}
}

Reproducible Example

import h3
import json
import h3pandas
import geopandas as gpd

# example raw data
gjson_index = ['87266c969ffffff', '87266c96dffffff', '87266ca05ffffff']
gjson_columns = ['REDACTED_0', 'REDACTED_1', 'REDACTED_2', 'REDACTED_3', 'REDACTED_4', 'REDACTED_5', 'REDACTED_6', 'REDACTED_7', 'REDACTED_8', 'REDACTED_9', 'REDACTED_10', 'REDACTED_11', 'LON', 'LAT', 'geometry']
gjson = '{"type": "FeatureCollection", "features": [{"id": "87266c969ffffff", "type": "Feature", "properties": {"LAT": 39.27600314991674, "LON": -85.1213943808014, "REDACTED_0": 977.5, "REDACTED_1": 222258708.75, "REDACTED_10": 0.4251362452703561, "REDACTED_11": 0.5354247602765141, "REDACTED_2": 4576063.0, "REDACTED_3": 83.93366542975889, "REDACTED_4": 0.007443973488911269, "REDACTED_5": 85.81456721436454, "REDACTED_6": 0.6711036424897768, "REDACTED_7": 0.3336976347908275, "REDACTED_8": 0.3264427707257549, "REDACTED_9": 0.3082428608333272}, "geometry": {"type": "Polygon", "coordinates": [[[-85.19683353530006, 39.3280507258189], [-85.20935486758405, 39.32032954238419], [-85.20606528270925, 39.30845430223158], [-85.19026103607621, 39.30429854736462], [-85.17774035424408, 39.312015943097265], [-85.18102326718702, 39.323892880693975], [-85.19683353530006, 39.3280507258189]]]}}, {"id": "87266c96dffffff", "type": "Feature", "properties": {"LAT": 39.41170134573524, "LON": -85.0650599012114, "REDACTED_0": 340.0, "REDACTED_1": 301859286.0, "REDACTED_10": 0.4673280928438077, "REDACTED_11": 0.5354862125261872, "REDACTED_2": 13778938.0, "REDACTED_3": 46.13529122132092, "REDACTED_4": 0.0032203361520541218, "REDACTED_5": 116.54852191580137, "REDACTED_6": 0.9122095125445141, "REDACTED_7": 0.35913638873824577, "REDACTED_8": 0.2945114671191823, "REDACTED_9": 0.2832131095896364}, "geometry": {"type": "Polygon", "coordinates": [[[-85.18759397216253, 39.347652019827095], [-85.2001222263731, 39.33992944658035], [-85.19683353530006, 39.3280507258189], [-85.18102326718702, 39.323892880693975], [-85.16849566636272, 39.33161166510846], [-85.17177767885819, 39.34349208277442], [-85.18759397216253, 39.347652019827095]]]}}, {"id": "87266ca05ffffff", "type": "Feature", "properties": {"LAT": 39.01771222409067, "LON": -85.17138141195356, "REDACTED_0": 1195.0, "REDACTED_1": 187713992.0, "REDACTED_10": 0.3725714305458182, "REDACTED_11": 0.4287774007011192, "REDACTED_2": 373883.0, "REDACTED_3": 67.67685000034521, "REDACTED_4": 0.005627416717251525, "REDACTED_5": 72.47677750922183, "REDACTED_6": 0.5664695521911012, "REDACTED_7": 0.30906217864873503, "REDACTED_8": 0.2442274365257504, "REDACTED_9": 0.26295879389349797}, "geometry": {"type": "Polygon", "coordinates": [[[-85.25850094180004, 38.9708121796447], [-85.27091677416291, 38.963127276382664], [-85.2676369585769, 38.95131160796221], [-85.25194787402559, 38.94717914887098], [-85.23953266664613, 38.9548602858906], [-85.24280591745425, 38.96667764754856], [-85.25850094180004, 38.9708121796447]]]}}]}'
gjson_crs = "epsg:4326"

# generate gdf
test = gpd.GeoDataFrame.from_features(json.loads(gjson), columns=gjson_columns, crs=gjson_crs)
test.index = gjson_index

# execute smoothing
test = test.h3.k_ring_smoothing(k=2)

# plot
test.plot(column="REDACTED_11", cmap="RdYlGn", legend=True)
plt.title("Hex Grid: Smoothed REDACTED_11")

/home/REDACTED/Projects/h3_study/.venv/lib/python3.8/site-packages/h3pandas/h3pandas.py:166: FutureWarning: CRS mismatch between CRS of the passed geometries and 'crs'. Use 'GeoDataFrame.set_crs(crs, allow_override=True)' to overwrite CRS or 'GeoDataFrame.to_crs(crs)' to reproject geometries. CRS mismatch will raise an error in the future versions of GeoPandas.
lambda x: gpd.GeoDataFrame(x, crs="epsg:4326")

The resulting smoothed GDF does appear to be correct, though the warning appears only when using k and not weights. It should be noted that this is in Jupyter Lab (haven't tested in Jupyter Notebook), though from the error I can't see why that makes a difference.

[johnziebro]

Interesting... running k_ring_smoothing with weights with only 1 coefficient, or multiple coefficients of the same value produces the same error, but using multiple coefficients of differing values does not.

I saw in the module that it defers to the less computationally expensive method of smoothing if possible, so this might help narrow it down to the clause that is less expensive.

[johnziebro]

The following works without a warning, the only difference is that geometry is specified as not returned.

test = (
    test.h3.k_ring_smoothing(3, return_geometry=False)
) 

Since the test gdf contains geometry, this points to a conflict when new geometry overwrites the existing geometry and lines up with the module code cited in the warning: lambda x: gpd.GeoDataFrame(x, crs="epsg:4326")

The issue occurs inside this module call:

# Unweighted case
        if weights is None:
            result = (
                self._df.h3.k_ring(k, explode=True)
                .groupby("h3_k_ring")
                .sum()
                .divide((1 + 3 * k * (k + 1)))
            )

            return result.h3.h3_to_geo_boundary() if return_geometry else result

Specifically from the last line of h3_to_geo_boundary() when the the gdf being smoothed already has geometry.

A workaround for now would be to execute smoothing using k param without returning geometry, and copy the already existing geometry column in the source gdf into the smoothed gdf, OR just extract the smoothed columns and copy them to the original gdf.

[DahnJ]

Great module by the way!

Thank you!

A workaround for now would be to execute smoothing using k param without returning geometry, and copy the already existing geometry column in the source gdf into the smoothed gdf, OR just extract the smoothed columns and copy them to the original gdf.

Thank you for the detailed analysis.

The warning happens because the grouped GeoDataFrame stays a GeoDataFrame even after groupby, but loses crs information after divide:

import geopandas as gpd

df = gpd.GeoDataFrame(
    {'group': [1, 1, 2], 'val': [5, 10, 15]}, 
    geometry=gpd.points_from_xy([15, 16, 17], [50, 51, 52], crs='epsg:4326')
)

# Groupby still retains CRS information
result = df.groupby('group').sum()
type(result)  # geopandas.geodataframe.GeoDataFrame
result.crs  # <Geographic 2D CRS: EPSG:4326>

# Divide afterwards does not
divided = result.divide(1)
type(divided)  # geopandas.geodataframe.GeoDataFrame
divided.crs  # None

I will consider opening an issue in GeoPandas, as I am not completely sure this behaviour makes sense.

A simple and backwards-compatible workaround for h3-pandas would be to cast the result of groupby().sum().divide() to a simple DataFrame, since that's what it effectively is, as there is no geometry.