{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "esacci-00",
   "metadata": {},
   "source": [
    "# ESA CCI land cover → WGS84 ellipsoidal HEALPix\n",
    "\n",
    "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.\n",
    "\n",
    "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/`.\n",
    "\n",
    "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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-01",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "import rioxarray\n",
    "import shapely\n",
    "import xarray as xr\n",
    "import zarr\n",
    "\n",
    "from healpix_convert.cache import create_staging_cache\n",
    "from healpix_convert.convert import convert_group_to_healpix, prepare_healpix_dataset\n",
    "from healpix_convert.settings.common import ConvertSettings\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "esacci-02",
   "metadata": {},
   "source": [
    "## Configuration\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-03",
   "metadata": {},
   "outputs": [],
   "source": [
    "ROOT = Path.cwd().parent if Path.cwd().name == \"notebooks\" else Path.cwd()\n",
    "SOURCE = ROOT / \"ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif\"\n",
    "FULL_GLOBE = False\n",
    "# The region controls output chunk coverage, not an exact polygon mask.\n",
    "REGION = shapely.box(2.0, 48.0, 2.02, 48.02)\n",
    "SUFFIX = \"global\" if FULL_GLOBE else \"trial\"\n",
    "STAGED = ROOT / f\"esacci-lc-2015-{SUFFIX}-source.zarr\"\n",
    "OUTPUT = ROOT / f\"esacci-lc-2015-{SUFFIX}-healpix-l15.zarr\"\n",
    "GROUP = \"measurements/land_cover\"\n",
    "\n",
    "settings = ConvertSettings.from_dict(\n",
    "    {\n",
    "        \"group_settings\": {\n",
    "            GROUP: {\n",
    "                \"healpix\": {\n",
    "                    \"refinement_level\": 15,\n",
    "                    \"indexing_scheme\": \"nested\",\n",
    "                    \"ellipsoid\": {\"name\": \"wgs84\"},\n",
    "                },\n",
    "                \"chunk\": {\n",
    "                    \"method\": \"healpix_cell_dense\",\n",
    "                    \"healpix\": {\n",
    "                        \"refinement_level\": 8,\n",
    "                        \"indexing_scheme\": \"nested\",\n",
    "                        \"ellipsoid\": {\"name\": \"wgs84\"},\n",
    "                    },\n",
    "                    \"chunk_buffer_width\": 1000.0,  # metres; halo beyond each chunk\n",
    "                },\n",
    "                \"resampler\": {\"name\": \"nearest\"},\n",
    "            }\n",
    "        }\n",
    "    }\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "esacci-04",
   "metadata": {},
   "source": [
    "## Step 1 — Stage the TIFF as a geographic Zarr group\n",
    "\n",
    "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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-05",
   "metadata": {},
   "outputs": [],
   "source": [
    "if STAGED.exists() or OUTPUT.exists():\n",
    "    raise FileExistsError(\"Choose fresh STAGED and OUTPUT paths before running\")\n",
    "\n",
    "raster = rioxarray.open_rasterio(\n",
    "    SOURCE, chunks={\"x\": 2048, \"y\": 2048}, masked=False, mask_and_scale=False\n",
    ")\n",
    "assert raster.sizes[\"band\"] == 1\n",
    "assert raster.rio.crs.to_epsg() == 4326\n",
    "assert raster.dtype == np.uint8\n",
    "assert raster.rio.nodata == 0\n",
    "if not FULL_GLOBE:\n",
    "    west, south, east, north = REGION.bounds\n",
    "    raster = raster.rio.clip_box(west - 1, south - 1, east + 1, north + 1)\n",
    "bounds = raster.rio.bounds()\n",
    "land_cover = raster.squeeze(\"band\", drop=True).drop_vars(\"spatial_ref\")\n",
    "land_cover = land_cover.rename({\"x\": \"longitude\", \"y\": \"latitude\"})\n",
    "land_cover.attrs = {\"long_name\": \"ESA CCI LCCS land-cover class\", \"units\": \"1\"}\n",
    "land_cover.encoding = {\"_FillValue\": None}\n",
    "land_cover.longitude.attrs = {\"standard_name\": \"longitude\", \"units\": \"degrees_east\"}\n",
    "land_cover.latitude.attrs = {\"standard_name\": \"latitude\", \"units\": \"degrees_north\"}\n",
    "ds = land_cover.to_dataset(name=\"lccs_class\")\n",
    "footprint = shapely.box(*bounds)\n",
    "root = xr.Dataset(\n",
    "    attrs={\n",
    "        \"stac_discovery\": {\n",
    "            \"type\": \"Feature\",\n",
    "            \"stac_version\": \"1.0.0\",\n",
    "            \"stac_extensions\": [],\n",
    "            \"id\": SOURCE.stem,\n",
    "            \"bbox\": list(bounds),\n",
    "            \"geometry\": shapely.geometry.mapping(footprint),\n",
    "            \"properties\": {\n",
    "                \"datetime\": None,\n",
    "                \"start_datetime\": \"2015-01-01T00:00:00Z\",\n",
    "                \"end_datetime\": \"2015-12-31T23:59:59Z\",\n",
    "                \"proj:epsg\": 4326,\n",
    "            },\n",
    "            \"links\": [],\n",
    "            \"assets\": {},\n",
    "        }\n",
    "    }\n",
    ")\n",
    "xr.DataTree.from_dict({\"/\": root, GROUP: ds}).to_zarr(STAGED, mode=\"w-\")\n",
    "raster.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "esacci-06",
   "metadata": {},
   "source": [
    "## Step 2 — Prepare the output and convert all spatial chunks\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the exact world extent to avoid floating-point TIFF bounds near +/-180.\n",
    "extent = shapely.box(-180, -90, 180, 90) if FULL_GLOBE else REGION\n",
    "cache = create_staging_cache([str(STAGED)], settings, output_extent=extent)\n",
    "prepared = prepare_healpix_dataset(cache, settings, str(OUTPUT))\n",
    "prepared.close()\n",
    "convert_group_to_healpix(\n",
    "    GROUP,\n",
    "    cache=cache,\n",
    "    settings=settings,\n",
    "    output_path=str(OUTPUT),\n",
    "    load_input_data=False,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "esacci-08",
   "metadata": {},
   "source": [
    "## Step 3 — Add categorical metadata and consolidate\n",
    "\n",
    "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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-09",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_group = zarr.open_group(str(OUTPUT), mode=\"a\")[GROUP]\n",
    "output_group[\"lccs_class\"].attrs.update(\n",
    "    {\n",
    "        \"long_name\": \"ESA CCI LCCS land-cover class\",\n",
    "        \"units\": \"1\",\n",
    "        \"_FillValue\": 0,\n",
    "        \"comment\": \"Categorical class IDs; nearest-neighbor resampling; 0 is no data\",\n",
    "    }\n",
    ")\n",
    "zarr.consolidate_metadata(str(OUTPUT))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "esacci-10",
   "metadata": {},
   "source": [
    "## Validation\n",
    "\n",
    "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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-11",
   "metadata": {},
   "outputs": [],
   "source": [
    "with xr.open_zarr(OUTPUT, group=GROUP, mask_and_scale=False) as result:\n",
    "    assert result.lccs_class.dtype == np.uint8\n",
    "    assert result.attrs[\"dggs\"][\"refinement_level\"] == 15\n",
    "    assert result.attrs[\"dggs\"][\"indexing_scheme\"] == \"nested\"\n",
    "    assert result.attrs[\"dggs\"][\"ellipsoid\"][\"name\"] == \"wgs84\"\n",
    "    sample = result.isel(cells=slice(0, 4 ** (15 - 8))).compute()\n",
    "    assert np.all(np.diff(sample.cell_ids.values) == 1)\n",
    "    print(sample)\n",
    "    print(\"Stored class codes in the first chunk:\", np.unique(sample.lccs_class))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2083d31a-98a2-4347-b679-7ecb12490694",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "esacci-visualization",
   "metadata": {},
   "source": [
    "## Plot the sample\n",
    "\n",
    "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.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-xdggs-plot",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xdggs\n",
    "from xdggs.plotting import explore\n",
    "\n",
    "grid_metadata = sample.attrs[\"dggs\"]\n",
    "grid_info = xdggs.HealpixInfo(\n",
    "    level=grid_metadata[\"refinement_level\"],\n",
    "    indexing_scheme=grid_metadata[\"indexing_scheme\"],\n",
    "    ellipsoid=grid_metadata[\"ellipsoid\"][\"name\"].upper(),\n",
    ")\n",
    "decoded = xdggs.decode(\n",
    "    sample, grid_info=grid_info, name=grid_metadata[\"coordinate\"]\n",
    ")\n",
    "class_codes = np.unique(decoded.lccs_class.values)\n",
    "xdggs_map = explore(\n",
    "    decoded.lccs_class,\n",
    "    colorize_params={\n",
    "        \"cmap\": \"tab20\",\n",
    "        \"vmin\": float(class_codes.min()),\n",
    "        \"vmax\": float(class_codes.max()),\n",
    "    },\n",
    ")\n",
    "xdggs_map\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "esacci-healpix-plot",
   "metadata": {},
   "outputs": [],
   "source": [
    "import healpix_plot\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import BoundaryNorm\n",
    "\n",
    "class_values = class_codes.astype(float)\n",
    "if class_values.size == 1:\n",
    "    class_boundaries = np.array([class_values[0] - 0.5, class_values[0] + 0.5])\n",
    "else:\n",
    "    class_boundaries = np.concatenate(\n",
    "        ([class_values[0] - 0.5], (class_values[:-1] + class_values[1:]) / 2, [class_values[-1] + 0.5])\n",
    "    )\n",
    "class_cmap = plt.colormaps[\"tab20\"].resampled(class_codes.size)\n",
    "class_norm = BoundaryNorm(class_boundaries, class_cmap.N)\n",
    "healpix_grid = healpix_plot.HealpixGrid(\n",
    "    level=grid_metadata[\"refinement_level\"],\n",
    "    indexing_scheme=grid_metadata[\"indexing_scheme\"],\n",
    "    ellipsoid=grid_metadata[\"ellipsoid\"][\"name\"].upper(),\n",
    ")\n",
    "\n",
    "mappable = healpix_plot.plot(\n",
    "    sample.cell_ids.values,\n",
    "    sample.lccs_class.values,\n",
    "    healpix_grid=healpix_grid,\n",
    "    sampling_grid={\"shape\": (900, 700)},\n",
    "    projection=\"PlateCarree\",\n",
    "    interpolation=\"nearest\",\n",
    "    cmap=class_cmap,\n",
    "    norm=class_norm,\n",
    "    colorbar={\"ticks\": class_codes, \"label\": \"ESA CCI LCCS class ID\", \"shrink\": 0.8},\n",
    "    title=\"ESA CCI land cover — WGS84 HEALPix level 15 (NESTED)\",\n",
    ")\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}