DGGS Resample

DGGS resampling utilities (area-weighted and nearest-neighbour transfer).

dggsresample_cli()

Command-line interface for :func:dggsresample.

Source code in vgrid/conversion/dggsresample/dggsresample.py

def dggsresample_cli():
    """Command-line interface for :func:`dggsresample`."""
    parser = argparse.ArgumentParser(
        description=(
            "Resample DGGS cells from one grid to another "
            "(area-weighted overlap or nearest-neighbour transfer; "
            "keep cells by centroid_within or intersects)"
        )
    )
    parser.add_argument(
        "-i",
        "--input",
        type=str,
        required=True,
        help="Source DGGS layer (vector file path, URL, GeoJSON, or GeoDataFrame path)",
    )
    parser.add_argument(
        "--from_dggs",
        "--dggs_from",
        dest="dggs_from",
        type=str,
        required=True,
        help=(
            "Source DGGS type (e.g. h3, s2, rhealpix, dggal_gnosis, dggrid_ISEA7H). "
            "Bare rhealpix is vgrid native; use dggal_rhealpix for DGGAL."
        ),
    )
    parser.add_argument(
        "--to_dggs",
        "--dggs_to",
        dest="dggs_to",
        type=str,
        required=True,
        help="Target DGGS type (same naming as --from_dggs)",
    )
    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        default=None,
        help="Target resolution; omit or pass -1 to match mean source cell area",
    )
    parser.add_argument(
        "-dggs_col",
        "--dggs_col",
        type=str,
        default=None,
        help="Source cell id column (default: same as --from_dggs)",
    )
    parser.add_argument(
        "-resample_col",
        "--resample_col",
        type=str,
        default=None,
        help="Numeric attribute to transfer; omit to return target grid only",
    )
    parser.add_argument(
        "-m",
        "--method",
        type=str,
        default="area_weighted",
        choices=[
            "area_weighted",
            "area",
            "nearest",
            "nn",
            "nearest_neighbour",
            "nearest_neighbor",
        ],
        help="Resampling method (default: area_weighted)",
    )
    parser.add_argument(
        "-p",
        "--predicate",
        type=str,
        default="centroid_within",
        choices=["centroid_within", "intersects"],
        help=(
            "Source-target keep filter (default: centroid_within). "
            "centroid_within: keep target cells that contain a source centroid; "
            "intersects: keep target cells that intersect any source cell"
        ),
    )
    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        default="gpd",
        choices=OUTPUT_FORMATS,
        help="Output format (default: gpd)",
    )
    parser.add_argument(
        "-o",
        "--output_name",
        type=str,
        default=None,
        help="Output base name or file path for file-based formats",
    )
    parser.add_argument(
        "-fix",
        "--fix_antimeridian",
        type=str,
        choices=[
            "shift",
            "shift_balanced",
            "shift_west",
            "shift_east",
            "split",
            "none",
        ],
        default=None,
        help="Antimeridian fixing method for supported grids",
    )
    parser.add_argument(
        "-split",
        "--split_antimeridian",
        action="store_true",
        default=False,
        help="Split geometries at the antimeridian",
    )
    parser.add_argument(
        "-aggregate",
        "--aggregate",
        action="store_true",
        default=False,
        help="For DGGRID targets: dissolve split geometries by global_id",
    )
    parser.add_argument(
        "--dggrid_options",
        type=str,
        default=None,
        help="JSON options for DGGRID grid generation (e.g. '{\"densification\": 2}')",
    )

    parser.add_argument(
        "--a5_options",
        "--a5_options",
        type=str,
        default=None,
        help="JSON options for A5 grid generation (e.g. '{\"segments\": 1000}')",
    )

    add_verbose_argument(parser)
    args = parser.parse_args()

    a5_options = None
    if args.a5_options:
        try:
            a5_options = json.loads(args.a5_options)
        except json.JSONDecodeError as exc:
            print(f"Error: Invalid JSON in a5_options: {exc}", file=sys.stderr)
            sys.exit(1)

    dggrid_options = None
    if args.dggrid_options:
        try:
            dggrid_options = json.loads(args.dggrid_options)
        except json.JSONDecodeError as exc:
            print(f"Error: Invalid JSON in dggrid_options: {exc}", file=sys.stderr)
            sys.exit(1)

    try:
        result = dggsresample(
            source_dggs=args.input,
            dggs_from=args.dggs_from,
            dggs_to=args.dggs_to,
            resolution=args.resolution,
            dggs_col=args.dggs_col,
            resample_col=args.resample_col,
            output_format=args.output_format,
            output_name=args.output_name,
            method=args.method,
            predicate=args.predicate,
            fix_antimeridian=args.fix_antimeridian,
            split_antimeridian=args.split_antimeridian,
            aggregate=args.aggregate,
            dggrid_options=dggrid_options,
            a5_options=a5_options,  # for A5 grid generation
            verbose=args.verbose,
        )
        if args.output_format in STRUCTURED_FORMATS:
            print(result)
    except Exception as exc:
        print(f"Error: {exc}", file=sys.stderr)
        sys.exit(1)

generate_grid(source_gdf, to_dggs, resolution, *, fix_antimeridian=None, split_antimeridian=False, aggregate=False, dggrid_options=None, a5_options=None, verbose=True)

Build a target DGGS GeoDataFrame over the axis-aligned bounds of source_gdf (source_gdf.total_bounds), using bbox-scoped grid generators only.

Source code in vgrid/conversion/dggsresample/dggsresample.py

def generate_grid(
    source_gdf: gpd.GeoDataFrame,
    to_dggs: str,
    resolution: int,
    *,
    fix_antimeridian: Optional[str] = None,
    split_antimeridian: bool = False,
    aggregate: bool = False,
    dggrid_options: Optional[dict] = None,
    a5_options: Optional[dict] = None,
    verbose: bool = True,
) -> gpd.GeoDataFrame:
    """
    Build a target DGGS GeoDataFrame over the axis-aligned bounds of ``source_gdf``
    (``source_gdf.total_bounds``), using bbox-scoped grid generators only.
    """
    if source_gdf.empty:
        raise ValueError("No features provided for grid generation.")

    bbox = validate_bbox(list(source_gdf.total_bounds))

    if to_dggs == "h3":
        gdf = h3_grid_within_bbox(
            resolution, bbox, fix_antimeridian=fix_antimeridian, verbose=verbose
        )
    elif to_dggs == "s2":
        gdf = s2_grid(
            resolution, bbox, fix_antimeridian=fix_antimeridian, verbose=verbose
        )
    elif to_dggs == "a5":
        gdf = a5_grid(
            resolution,
            bbox,
            options=a5_options,
            split_antimeridian=split_antimeridian,
            verbose=verbose,
        )
    elif to_dggs == "rhealpix":
        gdf = rhealpix_grid_within_bbox(
            resolution, bbox, fix_antimeridian=fix_antimeridian, verbose=verbose
        )
    elif to_dggs == "isea4t":
        if platform.system() != "Windows":
            raise ValueError("isea4t grid generation requires Windows in this build.")
        gdf = isea4t_grid_within_bbox(
            resolution, bbox, fix_antimeridian=fix_antimeridian, verbose=verbose
        )
    elif to_dggs == "qtm":
        gdf = qtm_grid_within_bbox(resolution, bbox, verbose=verbose)
    elif to_dggs == "olc":
        gdf = olc_grid_within_bbox(resolution, bbox, verbose=verbose)
    elif to_dggs == "geohash":
        gdf = geohash_grid_within_bbox(resolution, bbox, verbose=verbose)
    elif to_dggs == "tilecode":
        gdf = tilecode_grid(resolution, bbox, verbose=verbose)
    elif to_dggs == "quadkey":
        gdf = quadkey_grid(resolution, bbox, verbose=verbose)
    elif (dt := _dggal_short_type(to_dggs)) is not None:
        dt = validate_dggal_type(dt)
        bbox_t = (float(bbox[0]), float(bbox[1]), float(bbox[2]), float(bbox[3]))
        use_split = split_antimeridian or bool(fix_antimeridian)
        gdf = dggalgen(
            dggs_type=dt,
            resolution=resolution,
            bbox=bbox_t,
            compact=False,
            output_format="gpd",
            split_antimeridian=use_split,
            verbose=verbose,
        )
        if gdf is None or getattr(gdf, "empty", True):
            raise ValueError(
                f"No DGGAL cells for type {to_dggs!r} ({dt}), resolution {resolution}, bbox {bbox_t}."
            )
    elif (dt := _dggrid_short_type(to_dggs)) is not None:
        dt = validate_dggrid_type(dt)
        dggrid_instance = create_dggrid_instance()
        use_split = split_antimeridian or bool(fix_antimeridian)
        gdf = dggrid_generate_polygons(
            dggrid_instance,
            dt,
            resolution,
            bbox,
            output_address_type="SEQNUM",
            split_antimeridian=use_split,
            aggregate=aggregate,
            options=dggrid_options,
        )
        if gdf is None or getattr(gdf, "empty", True):
            raise ValueError(
                f"No DGGRID cells for type {to_dggs!r} ({dt}), resolution {resolution}, bbox {bbox}."
            )
    else:
        raise ValueError(f"Unsupported DGGS type: {to_dggs}")

    return gdf

get_nearest_resolution(source_gdf, from_dggs, to_dggs, from_col=None)

Match mean cell area of the first source cell to the closest to_dggs resolution (same search as the QGIS plugin).

Source code in vgrid/conversion/dggsresample/dggsresample.py

def get_nearest_resolution(
    source_gdf: gpd.GeoDataFrame,
    from_dggs: str,
    to_dggs: str,
    from_col: Optional[str] = None,
) -> int:
    """
    Match mean cell area of the first source cell to the closest ``to_dggs``
    resolution (same search as the QGIS plugin).
    """
    if from_col is None:
        from_col = from_dggs

    if source_gdf.empty or from_col not in source_gdf.columns:
        raise ValueError(f"No features or missing <{from_col}> column.")

    from_dggs_id = source_gdf.iloc[0][from_col]
    if from_dggs_id is None or (
        isinstance(from_dggs_id, float) and pd.isna(from_dggs_id)
    ):
        raise ValueError(f"No valid DGGS IDs found in <{from_col}> column.")

    try:
        if from_dggs == "h3":
            from_resolution = h3.get_resolution(from_dggs_id)
            from_area = h3.average_hexagon_area(from_resolution, unit="m^2")

        elif from_dggs == "s2":
            s2_id = s2.CellId.from_token(from_dggs_id)
            from_resolution = s2_id.level()
            _, _, from_area, _ = s2_metrics(from_resolution)

        elif from_dggs == "a5":
            from_resolution = a5.get_resolution(a5.hex_to_u64(from_dggs_id))
            _, _, from_area, _ = a5_metrics(from_resolution)

        elif from_dggs == "rhealpix":
            rhealpix_uids = (from_dggs_id[0],) + tuple(map(int, from_dggs_id[1:]))
            rhealpix_dggs = RHEALPixDGGS(
                ellipsoid=E, north_square=1, south_square=3, N_side=3
            )
            rhealpix_cell = rhealpix_dggs.cell(rhealpix_uids)
            from_resolution = rhealpix_cell.resolution
            _, _, from_area, _ = rhealpix_metrics(from_resolution)

        elif from_dggs == "isea4t":
            if platform.system() != "Windows":
                raise ValueError(
                    "isea4t area matching is only available on Windows in this build.",
                )
            from_resolution = len(from_dggs_id) - 2
            _, _, from_area, _ = isea4t_metrics(from_resolution)

        elif from_dggs == "qtm":
            from_resolution = len(from_dggs_id)
            _, _, from_area, _ = qtm_metrics(from_resolution)

        elif from_dggs == "olc":
            coord = olc.decode(from_dggs_id)
            from_resolution = coord.codeLength
            _, _, from_area, _ = olc_metrics(from_resolution)

        elif from_dggs == "geohash":
            from_resolution = len(from_dggs_id)
            _, _, from_area, _ = geohash_metrics(from_resolution)

        elif from_dggs == "tilecode":
            match = re.match(r"z(\d+)x(\d+)y(\d+)", from_dggs_id)
            if not match:
                raise ValueError("tilecode id must look like z{x}x{y}y{y}")
            from_resolution = int(match.group(1))
            _, _, from_area, _ = tilecode_metrics(from_resolution)

        elif from_dggs == "quadkey":
            tile = mercantile.quadkey_to_tile(from_dggs_id)
            from_resolution = tile.z
            _, _, from_area, _ = quadkey_metrics(from_resolution)

        elif (dt := _dggal_short_type(from_dggs)) is not None:
            dt = validate_dggal_type(dt)
            cls_name = DGGAL_TYPES[dt]["class_name"]
            dggrs = getattr(dggal, cls_name)()
            zone = dggrs.getZoneFromTextID(str(from_dggs_id))
            from_resolution = int(dggrs.getZoneLevel(zone))
            _, _, from_area, _ = dggal_metrics(dt, from_resolution)

        elif (dt := _dggrid_short_type(from_dggs)) is not None:
            dt = validate_dggrid_type(dt)
            dggrid_instance = create_dggrid_instance()
            if "resolution" in source_gdf.columns:
                rv = source_gdf.iloc[0]["resolution"]
                if rv is None or (isinstance(rv, float) and pd.isna(rv)):
                    raise ValueError(
                        "Missing or invalid <resolution> for DGGRID source."
                    )
                from_resolution = int(rv)
            else:
                pref = int(DGGRID_TYPES[dt]["default_res"])
                _, from_resolution = _resolve_cell_geometry(
                    dggrid_instance,
                    dt,
                    str(from_dggs_id),
                    preferred_resolution=pref,
                )
                if from_resolution is None:
                    raise ValueError(
                        f"Could not resolve DGGRID cell id {from_dggs_id!r} for type {dt}."
                    )
            stats = dggridstats(dggrid_instance, dt, unit="m")
            row = stats[stats["resolution"] == from_resolution]
            if row.empty:
                raise ValueError(
                    f"No DGGRID stats for {dt} at resolution {from_resolution}."
                )
            from_area = float(row["area_m2"].iloc[0])

        else:
            raise ValueError(f"Unsupported from_dggs type for area match: {from_dggs}")

    except Exception as e:
        raise ValueError(f"Failed to calculate area from {from_dggs}: {e}") from e

    nearest_resolution: int
    min_diff = float("inf")

    try:
        if to_dggs == "h3":
            for res in range(16):
                avg_area = h3.average_hexagon_area(res, unit="m^2")
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "s2":
            for res in range(31):
                _, _, avg_area, _ = s2_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "a5":
            for res in range(30):
                _, _, avg_area, _ = a5_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "rhealpix":
            for res in range(16):
                _, _, avg_area, _ = rhealpix_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "isea4t":
            if platform.system() != "Windows":
                raise ValueError(
                    "isea4t nearest resolution is only available on Windows in this build.",
                )
            for res in range(26):
                _, _, avg_area, _ = isea4t_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "qtm":
            for res in range(1, 25):
                _, _, avg_area, _ = qtm_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "olc":
            for res in olc_resolutions:
                _, _, avg_area, _ = olc_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "geohash":
            for res in range(1, 11):
                _, _, avg_area, _ = geohash_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "tilecode":
            for res in range(30):
                _, _, avg_area, _ = tilecode_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif to_dggs == "quadkey":
            for res in range(30):
                _, _, avg_area, _ = quadkey_metrics(res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif (dt := _dggal_short_type(to_dggs)) is not None:
            dt = validate_dggal_type(dt)
            lo = int(DGGAL_TYPES[dt]["min_res"])
            hi = int(DGGAL_TYPES[dt]["max_res"])
            for res in range(lo, hi + 1):
                _, _, avg_area, _ = dggal_metrics(dt, res)
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        elif (dt := _dggrid_short_type(to_dggs)) is not None:
            dt = validate_dggrid_type(dt)
            dggrid_instance = create_dggrid_instance()
            stats = dggridstats(dggrid_instance, dt, unit="m")
            lo = int(DGGRID_TYPES[dt]["min_res"])
            hi = int(DGGRID_TYPES[dt]["max_res"])
            for _, row in stats.iterrows():
                res = int(row["resolution"])
                if res < lo or res > hi:
                    continue
                avg_area = float(row["area_m2"])
                diff = abs(avg_area - from_area)
                if diff < min_diff:
                    min_diff = diff
                    nearest_resolution = res

        else:
            raise ValueError(f"Unsupported to_dggs type for area match: {to_dggs}")

    except Exception as e:
        raise ValueError(
            f"Failed to calculate nearest resolution for {to_dggs}: {e}"
        ) from e

    return nearest_resolution

process_input_data_resample(input_data, dggs_col=None, resample_col=None, crs='EPSG:4326', **kwargs)

Load source DGGS layers for resampling with polygon geometry preserved.

Supports: - GeoDataFrame or DataFrame with a geometry column - GeoJSON dictionary (FeatureCollection) or list of features - Local or remote vector files (GeoJSON, Shapefile, GPKG, etc.) - CSV / Parquet tables that include a geometry column

Parameters

input_data Source DGGS cell layer in any supported form. dggs_col : str, optional If given, column holding DGGS cell identifiers is validated. resample_col : str, optional Numeric attribute column to resample; if given, must be present. crs : str, default "EPSG:4326" Target CRS for the returned GeoDataFrame. **kwargs Passed to :func:geopandas.read_file for file paths.

Returns

geopandas.GeoDataFrame Source cells with geometry and requested attribute columns.

Source code in vgrid/utils/io.py

def process_input_data_resample(
    input_data,
    dggs_col=None,
    resample_col=None,
    crs="EPSG:4326",
    **kwargs,
):
    """
    Load source DGGS layers for resampling with polygon geometry preserved.

    Supports:
    - GeoDataFrame or DataFrame with a geometry column
    - GeoJSON dictionary (FeatureCollection) or list of features
    - Local or remote vector files (GeoJSON, Shapefile, GPKG, etc.)
    - CSV / Parquet tables that include a geometry column

    Parameters
    ----------
    input_data
        Source DGGS cell layer in any supported form.
    dggs_col : str, optional
        If given, column holding DGGS cell identifiers is validated.
    resample_col : str, optional
        Numeric attribute column to resample; if given, must be present.
    crs : str, default ``\"EPSG:4326\"``
        Target CRS for the returned GeoDataFrame.
    **kwargs
        Passed to :func:`geopandas.read_file` for file paths.

    Returns
    -------
    geopandas.GeoDataFrame
        Source cells with geometry and requested attribute columns.
    """
    if isinstance(input_data, str):
        low = input_data.lower()
        if low.endswith(".csv"):
            df = pd.read_csv(input_data, **kwargs)
            if "geometry" not in df.columns:
                raise ValueError(
                    f"CSV file must include a 'geometry' column for resampling: {input_data}"
                )
            gdf = gpd.GeoDataFrame(df, geometry="geometry", crs=crs)
        elif low.endswith(".parquet"):
            df = pd.read_parquet(input_data, **kwargs)
            if "geometry" not in df.columns:
                raise ValueError(
                    f"Parquet file must include a 'geometry' column for resampling: "
                    f"{input_data}"
                )
            gdf = gpd.GeoDataFrame(df, geometry="geometry", crs=crs)
        else:
            gdf = process_input_data_vector(input_data, crs=crs, **kwargs)
    else:
        gdf = process_input_data_vector(input_data, crs=crs, **kwargs)

    if dggs_col is not None and dggs_col not in gdf.columns:
        raise ValueError(f"Missing '{dggs_col}' in input data.")
    if resample_col is not None and resample_col not in gdf.columns:
        raise ValueError(f"Missing '{resample_col}' in input data.")

    if gdf.crs is None:
        gdf = gdf.set_crs(crs)
    elif str(gdf.crs) != crs:
        gdf = gdf.to_crs(crs)

    return _require_resample_geometry(gdf)

resampling(source_gdf, target_gdf, resample_col, method='nearest', verbose=True, predicate='centroid_within')

Transfer resample_col from source cells onto target cells.

Parameters

method "area_weighted" (default) — for each target geometry::

    value = sum(source_value * area(intersection) / area(target_cell))

Only target rows that pass ``predicate`` are returned.

``"nearest"`` — assign the value of the source cell whose centroid is
closest to the target centroid (``geopandas.sjoin_nearest``).
Only target rows that pass ``predicate`` are returned.

predicate Source–target keep filter applied before value transfer:

``"centroid_within"`` (default) — keep a target cell if it contains
at least one source cell centroid.

``"intersects"`` — keep a target cell if it intersects at least one
source cell.

Source code in vgrid/conversion/dggsresample/dggsresample.py

def resampling(
    source_gdf: gpd.GeoDataFrame,
    target_gdf: gpd.GeoDataFrame,
    resample_col: str,
    method: str = "nearest",
    verbose=True,
    predicate: str = "centroid_within",
) -> gpd.GeoDataFrame:
    """
    Transfer ``resample_col`` from source cells onto target cells.

    Parameters
    ----------
    method
        ``\"area_weighted\"`` (default) — for each target geometry::

            value = sum(source_value * area(intersection) / area(target_cell))

        Only target rows that pass ``predicate`` are returned.

        ``\"nearest\"`` — assign the value of the source cell whose centroid is
        closest to the target centroid (``geopandas.sjoin_nearest``).
        Only target rows that pass ``predicate`` are returned.
    predicate
        Source–target keep filter applied before value transfer:

        ``\"centroid_within\"`` (default) — keep a target cell if it contains
        at least one source cell centroid.

        ``\"intersects\"`` — keep a target cell if it intersects at least one
        source cell.
    """
    if resample_col not in source_gdf.columns:
        raise ValueError(
            f"There is no <{resample_col}> column in the source GeoDataFrame."
        )

    _normalize_predicate(predicate)

    norm = method.strip().lower().replace("-", "_")
    if norm in ("area_weighted", "area"):
        return _resampling_area_weighted(
            source_gdf,
            target_gdf,
            resample_col,
            verbose=verbose,
            predicate=predicate,
        )
    if norm in ("nearest", "nn", "nearest_neighbour", "nearest_neighbor"):
        return _resampling_nearest(
            source_gdf,
            target_gdf,
            resample_col,
            verbose=verbose,
            predicate=predicate,
        )

    raise ValueError(
        f"Unsupported resampling method {method!r}; use 'area_weighted' or 'nearest'."
    )