2

I am trying to create a map out of a file with millions of points. To make this tractable, I am using h3 to aggregate the points. After doing the aggregation using a combination of h3 and pandas in Python (adapted from here), I am converting the h3 indexes to polygon boundaries:

agg1=gdf.groupby(['h3', 'product', 'year']).agg(Total_Count=('product', 'sum')).groupby(level=0).apply(lambda x: 100*x / x.sum()).reset_index()

def add_geometry(row):
    points = h3.h3_to_geo_boundary(
      row['h3'], True)
    return Polygon(points)


agg1['geometry'] = agg1.apply(add_geometry, axis=1)


gdf_op = gpd.GeoDataFrame(agg1, crs='EPSG:4326')
output_filename = 'gridcounts.gpkg'
gdf_op.to_file(driver='GPKG', filename=output_filename)

However, when I open the file in QGIS, some hexagons are streaking across the map: enter image description here

I have a feeling it is because h3 is defined on a globe and as a result the polygons towards the edges translate to coordinate bounds which are at the opposite ends of the projected map. Does anyone have a solution to this apart from deleting the problem hexagons?

alphabetasoup
  • 8,718
  • 4
  • 38
  • 78
DotPi
  • 779
  • 1
  • 8
  • 24
  • Check the code from https://gis.stackexchange.com/questions/345086/gdal-vectortranslate-creates-shards-fragments-along-anti-meridian , you can copy the snippet and apply the idl_resolve function and see if that helps. – spatialthoughts Sep 07 '21 at 15:04

2 Answers2

0

It is a projection issue. The hexagon shape is being split at the antimeridian (180 degrees). Try converting your output to a 3D coordinate system and plot in an orographic view. That should do it!

Tim Loftus
  • 1
  • 1
0

An alternative would be to represent your H3 cells as points, rather than as hexagons. That is, use h3_to_geo rather than h3_to_geo_boundary. Points won't be stretched across the antimeridian like this.

A further option is to drop the affected cells, and draw two maps each with a different central meridian.

alphabetasoup
  • 8,718
  • 4
  • 38
  • 78