Checkout antemeridian line for h3 tiles

[GriffinBabe]

It was notified by Daniele that the h3 library has issues on the hexagonal tiles close to the antemeridian line. There should be an investigation of the split_h3_tiles function as well a fix to implement if necessary.

[GriffinBabe]

Daniele managed to solve this issue previously by reconstructing the grid with this solution: https://github.com/VITO-RS-Vegetation/veg-workflows/blob/main/workflows/grids/h3_grid/build_h3_grid.py#L17

[dzanaga]

# build h3 grid
import geopandas as gpd
import h3
import hvplot.pandas
import shapely
from shapely.geometry import LineString, MultiPolygon, Polygon
from shapely.ops import split, unary_union
from tqdm import tqdm

# Get the boundary coordinates for each cell
def get_boundary_coordinates(h3_index):
    lat_lon_boundary = h3.h3_to_geo_boundary(h3_index, geo_json=True)
    return Polygon(lat_lon_boundary)

def check_and_fix_antimeridian_crossing(polygon):
    # A line that we will use to split the polygon at the antimeridian
    eps = 0
    splitter = LineString([(-180, 90), (-180, -90)])

    # Check if the polygon crosses the antimeridian by looking for a large jump in longitude
    coords = list(polygon.exterior.coords)
    crosses_antimeridian = any(
        (x1 > 100 and x2 < -100) or (x1 < -100 and x2 > 100)
        for (x1, y1), (x2, y2) in zip(coords, coords[1:] + [coords[0]])
    )

    if crosses_antimeridian:
        # Normalize coordinates to split across the antimeridian
        normalized_coords = [(x - 360 if x > 0 else x, y) for x, y in coords]
        normalized_polygon = Polygon(normalized_coords)

        # Split the polygon using the splitter line
        split_result = split(normalized_polygon, splitter)

        # # Filter out non-polygon geometries (like LinearRings or Points that might be created)
        # split_polygons = [geom for geom in split_result if isinstance(geom, Polygon)]

        # Normalize the polygons to have their coordinates within the range (-180, 180)
        fixed_polygons = []
        for poly in split_result.geoms:
            new_coords = []
            if poly.centroid.x < -180:
                new_coords = [(x + 360, y) for x, y in poly.exterior.coords]
                new_coords = [(x + eps, y) if x == -180 else (x, y)
                              for x, y in new_coords]
                new_coords = [(x - eps, y) if x == 180 else (x, y)
                              for x, y in new_coords]
                fixed_polygons.append(Polygon(new_coords))
            else:
                new_coords = [(x, y) for x, y in poly.exterior.coords]
                new_coords = [(x + eps, y) if x == -180 else (x, y)
                              for x, y in new_coords]
                new_coords = [(x - eps, y) if x == 180 else (x, y)
                              for x, y in new_coords]
                fixed_polygons.append(Polygon(new_coords))
                # fixed_polygons.append(poly)

        # Return the resulting polygons as a MultiPolygon
        new_polygon = MultiPolygon(fixed_polygons)

        # Return the resulting polygons as a MultiPolygon
        return new_polygon
    else:
        return polygon  # Return the original polygon if it doesn't cross the antimeridian

def get_cell_id(centroid):
    x, y = centroid.x, centroid.y
    if x < 0:
        east_west = 'W'
    else:
        east_west = 'E'

    if y < 0:
        north_south = 'S'
    else:
        north_south = 'N'

    return f"{east_west}{abs(round(x)):03d}{north_south}{abs(round(y)):02d}"

def _cast_to_multi(geom):
    if isinstance(geom, Polygon):
        return MultiPolygon([geom])
    else:
        return geom

if __name__ == '__main__':
    import argparse

    from loguru import logger

    parser = argparse.ArgumentParser()
    parser.add_argument('resolution', type=int)

    args = parser.parse_args()

    # Get all H3 cells at resolution
    resolution = args.resolution

    logger.info(f"Building H3 grid at resolution {resolution}")
    h3_grid_fn = f'/vitodata/vegteam/auxdata/grid/h3/h3_grid_res{resolution}.fgb'

    logger.info(f"Saving H3 grid to {h3_grid_fn}")

    all_cells = list(h3.get_res0_indexes())
    all_cells_res2 = [h3.h3_to_children(cell, resolution)
                      for cell in tqdm(all_cells)]
    # Flatten the list of sets into a unique set
    all_cells_res2 = set().union(*all_cells_res2)

    # Create a GeoDataFrame
    gdf = gpd.GeoDataFrame({
        'geometry': [get_boundary_coordinates(cell)
                     for cell in tqdm(all_cells_res2)],
        'h3_index': list(all_cells_res2)
    }, crs=4326)

    # Example usage:
    polygons = gdf.geometry.values  # Define your polygon here
    new = [check_and_fix_antimeridian_crossing(
        p) for p in tqdm(polygons.tolist())]

    gdf_fix = gdf.copy()
    gdf_fix['geometry'] = gdf.apply(
        lambda row: check_and_fix_antimeridian_crossing(row.geometry), axis=1).values
    # cell_id from centroid

    gdf_fix['cell_id'] = gdf_fix.centroid.apply(get_cell_id)
    gdf_fix.shape[0], gdf_fix.value_counts('cell_id').max()
    gdf_fix = gdf_fix[['cell_id', 'h3_index', 'geometry']]
    gdf_fix['geometry'] = gdf_fix.geometry.apply(_cast_to_multi)

    gdf_fix.to_file(h3_grid_fn)