ESA CCI land cover → WGS84 ellipsoidal HEALPix#
Convert ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif using nearest-neighbor, NESTED indexing and level 15. Class IDs are categorical: never use bilinear, PSF, averaging or any other interpolation. Both the output and chunk grids explicitly use the WGS84 ellipsoid.
Install this repository and the visualization dependencies in the notebook kernel environment: python -m pip install -e . rioxarray xdggs healpix-plot. Keep the TIFF locally in the repository root; it is ignored by git. Run from the repository root or notebooks/.
The default is a small regional trial using the actual level-15 grid. Set FULL_GLOBE = True for the entire TIFF. A global level-15 grid has 12,884,901,888 cells: about 12 GiB of uint8 classes plus 96 GiB of explicit uint64 cell IDs before compression, in addition to the staged source. Full conversion takes substantial disk space and time. The code stages and converts in chunks; it never loads the whole TIFF into RAM.
from pathlib import Path
import numpy as np
import rioxarray
import shapely
import xarray as xr
import zarr
from healpix_convert.cache import create_staging_cache
from healpix_convert.convert import convert_group_to_healpix, prepare_healpix_dataset
from healpix_convert.settings.common import ConvertSettings
Configuration#
ROOT = Path.cwd().parent if Path.cwd().name == "notebooks" else Path.cwd()
SOURCE = ROOT / "ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif"
FULL_GLOBE = False
# The region controls output chunk coverage, not an exact polygon mask.
REGION = shapely.box(2.0, 48.0, 2.02, 48.02)
SUFFIX = "global" if FULL_GLOBE else "trial"
STAGED = ROOT / f"esacci-lc-2015-{SUFFIX}-source.zarr"
OUTPUT = ROOT / f"esacci-lc-2015-{SUFFIX}-healpix-l15.zarr"
GROUP = "measurements/land_cover"
settings = ConvertSettings.from_dict(
{
"group_settings": {
GROUP: {
"healpix": {
"refinement_level": 15,
"indexing_scheme": "nested",
"ellipsoid": {"name": "wgs84"},
},
"chunk": {
"method": "healpix_cell_dense",
"healpix": {
"refinement_level": 8,
"indexing_scheme": "nested",
"ellipsoid": {"name": "wgs84"},
},
"chunk_buffer_width": 1000.0, # metres; halo beyond each chunk
},
"resampler": {"name": "nearest"},
}
}
}
)
Step 1 — Stage the TIFF as a geographic Zarr group#
Keep no-data code 0 as an integer category during resampling, so holes are not filled from adjacent classes. Do not enable masking or scale/offset decoding. Pixel-centre coordinates come from the TIFF affine transform. The trial retains a one-degree input margin, covering the selected level-8 chunks and their 1 km buffer.
if STAGED.exists() or OUTPUT.exists():
raise FileExistsError("Choose fresh STAGED and OUTPUT paths before running")
raster = rioxarray.open_rasterio(
SOURCE, chunks={"x": 2048, "y": 2048}, masked=False, mask_and_scale=False
)
assert raster.sizes["band"] == 1
assert raster.rio.crs.to_epsg() == 4326
assert raster.dtype == np.uint8
assert raster.rio.nodata == 0
if not FULL_GLOBE:
west, south, east, north = REGION.bounds
raster = raster.rio.clip_box(west - 1, south - 1, east + 1, north + 1)
bounds = raster.rio.bounds()
land_cover = raster.squeeze("band", drop=True).drop_vars("spatial_ref")
land_cover = land_cover.rename({"x": "longitude", "y": "latitude"})
land_cover.attrs = {"long_name": "ESA CCI LCCS land-cover class", "units": "1"}
land_cover.encoding = {"_FillValue": None}
land_cover.longitude.attrs = {"standard_name": "longitude", "units": "degrees_east"}
land_cover.latitude.attrs = {"standard_name": "latitude", "units": "degrees_north"}
ds = land_cover.to_dataset(name="lccs_class")
footprint = shapely.box(*bounds)
root = xr.Dataset(
attrs={
"stac_discovery": {
"type": "Feature",
"stac_version": "1.0.0",
"stac_extensions": [],
"id": SOURCE.stem,
"bbox": list(bounds),
"geometry": shapely.geometry.mapping(footprint),
"properties": {
"datetime": None,
"start_datetime": "2015-01-01T00:00:00Z",
"end_datetime": "2015-12-31T23:59:59Z",
"proj:epsg": 4326,
},
"links": [],
"assets": {},
}
}
)
xr.DataTree.from_dict({"/": root, GROUP: ds}).to_zarr(STAGED, mode="w-")
raster.close()
Step 2 — Prepare the output and convert all spatial chunks#
# Specify the exact world extent to avoid floating-point TIFF bounds near +/-180.
extent = shapely.box(-180, -90, 180, 90) if FULL_GLOBE else REGION
cache = create_staging_cache([str(STAGED)], settings, output_extent=extent)
prepared = prepare_healpix_dataset(cache, settings, str(OUTPUT))
prepared.close()
convert_group_to_healpix(
GROUP,
cache=cache,
settings=settings,
output_path=str(OUTPUT),
load_input_data=False,
)
Step 3 — Add categorical metadata and consolidate#
The generic converter currently writes grid metadata only. Restore the class description and declare 0 as no-data in the output. Uncovered output cells also use this code. Inspect with mask_and_scale=False to retain raw uint8 codes; normal CF decoding may mask 0 and promote the array to floating point without changing the stored data.
output_group = zarr.open_group(str(OUTPUT), mode="a")[GROUP]
output_group["lccs_class"].attrs.update(
{
"long_name": "ESA CCI LCCS land-cover class",
"units": "1",
"_FillValue": 0,
"comment": "Categorical class IDs; nearest-neighbor resampling; 0 is no data",
}
)
zarr.consolidate_metadata(str(OUTPUT))
Validation#
Inspect one output chunk without loading the global product. The regression suite also compares a small level-15 conversion against a direct nearest-source lookup.
with xr.open_zarr(OUTPUT, group=GROUP, mask_and_scale=False) as result:
assert result.lccs_class.dtype == np.uint8
assert result.attrs["dggs"]["refinement_level"] == 15
assert result.attrs["dggs"]["indexing_scheme"] == "nested"
assert result.attrs["dggs"]["ellipsoid"]["name"] == "wgs84"
sample = result.isel(cells=slice(0, 4 ** (15 - 8))).compute()
assert np.all(np.diff(sample.cell_ids.values) == 1)
print(sample)
print("Stored class codes in the first chunk:", np.unique(sample.lccs_class))
sample
Plot the sample#
Decode the Zarr DGGS convention with XDGGS for an interactive polygon map, then render the same categorical values with healpix-plot. The static plot uses nearest-neighbor sampling and a discrete normalization, so class IDs are not interpolated.
import xdggs
from xdggs.plotting import explore
grid_metadata = sample.attrs["dggs"]
grid_info = xdggs.HealpixInfo(
level=grid_metadata["refinement_level"],
indexing_scheme=grid_metadata["indexing_scheme"],
ellipsoid=grid_metadata["ellipsoid"]["name"].upper(),
)
decoded = xdggs.decode(
sample, grid_info=grid_info, name=grid_metadata["coordinate"]
)
class_codes = np.unique(decoded.lccs_class.values)
xdggs_map = explore(
decoded.lccs_class,
colorize_params={
"cmap": "tab20",
"vmin": float(class_codes.min()),
"vmax": float(class_codes.max()),
},
)
xdggs_map
import healpix_plot
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
class_values = class_codes.astype(float)
if class_values.size == 1:
class_boundaries = np.array([class_values[0] - 0.5, class_values[0] + 0.5])
else:
class_boundaries = np.concatenate(
([class_values[0] - 0.5], (class_values[:-1] + class_values[1:]) / 2, [class_values[-1] + 0.5])
)
class_cmap = plt.colormaps["tab20"].resampled(class_codes.size)
class_norm = BoundaryNorm(class_boundaries, class_cmap.N)
healpix_grid = healpix_plot.HealpixGrid(
level=grid_metadata["refinement_level"],
indexing_scheme=grid_metadata["indexing_scheme"],
ellipsoid=grid_metadata["ellipsoid"]["name"].upper(),
)
mappable = healpix_plot.plot(
sample.cell_ids.values,
sample.lccs_class.values,
healpix_grid=healpix_grid,
sampling_grid={"shape": (900, 700)},
projection="PlateCarree",
interpolation="nearest",
cmap=class_cmap,
norm=class_norm,
colorbar={"ticks": class_codes, "label": "ESA CCI LCCS class ID", "shrink": 0.8},
title="ESA CCI land cover — WGS84 HEALPix level 15 (NESTED)",
)
plt.show()