Support string/categorical columns in rasterize with QGIS-visible labels

docs/source/user_guide/rasterize.ipynb

@@ -3,7 +3,7 @@
   {
    "cell_type": "markdown",
    "id": "jb1yj58wq4q",
-   "source": "## Rasterize\n\n`xrspatial.rasterize` converts vector geometries (polygons, lines, points) into a 2D `xr.DataArray`. No GDAL dependency required.\n\nThis guide covers:\n- [Basic rasterization](#Basic-rasterization) -- polygons, lines, and points\n- [Merge modes](#Merge-modes) -- controlling how overlapping geometries combine\n- [Custom merge functions](#Custom-merge-functions) -- user-defined numba-jitted merge logic\n- [Dask parallel rasterization](#Dask-parallel-rasterization) -- tile-based output chunking",
+   "source": "## Rasterize\n\n`xrspatial.rasterize` converts vector geometries (polygons, lines, points) into a 2D `xr.DataArray`. No GDAL dependency required.\n\nThis guide covers:\n- [Basic rasterization](#Basic-rasterization) -- polygons, lines, and points\n- [Categorical columns](#Categorical-columns) -- burn string/categorical labels, readable in QGIS\n- [Merge modes](#Merge-modes) -- controlling how overlapping geometries combine\n- [Custom merge functions](#Custom-merge-functions) -- user-defined numba-jitted merge logic\n- [Dask parallel rasterization](#Dask-parallel-rasterization) -- tile-based output chunking",
    "metadata": {}
   }
  ],

examples/user_guide/28_Rasterize.ipynb

@@ -18,20 +18,21 @@
     "### What you'll build\n",
     "\n",
     "1. Rasterize land-use zones with the `.xrs` accessor\n",
-    "2. Handle overlapping polygons and interior holes\n",
-    "3. Burn lines and points into a raster\n",
-    "4. Compare merge modes for overlapping features\n",
-    "5. Write a custom numba merge function\n",
-    "6. Use multi-column properties for density mapping\n",
-    "7. Run Dask parallel rasterization\n",
-    "8. Combine rasterization with zonal statistics\n",
-    "9. Compare default vs. `all_touched` rasterization\n",
-    "10. Use the standalone `rasterize()` function with geometry pairs\n",
+    "2. Burn a string or categorical column into a labeled raster\n",
+    "3. Handle overlapping polygons and interior holes\n",
+    "4. Burn lines and points into a raster\n",
+    "5. Compare merge modes for overlapping features\n",
+    "6. Write a custom numba merge function\n",
+    "7. Use multi-column properties for density mapping\n",
+    "8. Run Dask parallel rasterization\n",
+    "9. Combine rasterization with zonal statistics\n",
+    "10. Compare default vs. `all_touched` rasterization\n",
+    "11. Use the standalone `rasterize()` function with geometry pairs\n",
     "\n",
     "![Rasterize preview](images/rasterize_preview.png)\n",
     "\n",
     "**Jump to a section:**\n",
-    "[Basic rasterization](#Basic-rasterization) | [Overlapping polygons](#Overlapping-polygons) | [Lines and points](#Lines-and-points) | [Merge modes](#Merge-modes) | [Custom merge](#Custom-merge) | [Multi-column properties](#Multi-column-properties) | [Dask parallel](#Dask-parallel) | [Zonal statistics](#Zonal-statistics) | [All touched](#All-touched) | [Standalone function](#Standalone-function)"
+    "[Basic rasterization](#Basic-rasterization) | [Categorical columns](#Categorical-columns) | [Overlapping polygons](#Overlapping-polygons) | [Lines and points](#Lines-and-points) | [Merge modes](#Merge-modes) | [Custom merge](#Custom-merge) | [Multi-column properties](#Multi-column-properties) | [Dask parallel](#Dask-parallel) | [Zonal statistics](#Zonal-statistics) | [All touched](#All-touched) | [Standalone function](#Standalone-function)"
    ]
   },
   {
@@ -120,6 +121,58 @@
     "plt.tight_layout()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "cat3482md",
+   "metadata": {},
+   "source": [
+    "## Categorical columns\n",
+    "\n",
+    "Pass a string or categorical column straight to `column=` and `rasterize` label-encodes it for you. Each distinct label becomes an integer code, the output is an `int32` band with a `-1` nodata value, and the label map rides along in `result.attrs['category_names']` (the list index is the pixel code). A matching `attrs['category_colors']` holds one RGBA per class.\n",
+    "\n",
+    "The plot below burns the land-use `label` column directly, with no manual encoding step."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "cat3482code",
+   "metadata": {},
+   "source": [
+    "# Burn the string 'label' column directly -- no manual encoding\n",
+    "landcover = template.xrs.rasterize(gdf, column='label')\n",
+    "\n",
+    "names = landcover.attrs['category_names']\n",
+    "print('category_names:', names)\n",
+    "print('dtype:', landcover.dtype, '| nodata:', landcover.attrs['nodata'])\n",
+    "\n",
+    "# Colors the encoder assigned, scaled to 0-1 for matplotlib\n",
+    "colors = [tuple(c / 255 for c in rgba)\n",
+    "          for rgba in landcover.attrs['category_colors']]\n",
+    "cmap = ListedColormap(colors)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10, 4))\n",
+    "landcover.where(landcover >= 0).plot.imshow(\n",
+    "    ax=ax, cmap=cmap, add_colorbar=False, vmin=0, vmax=len(names) - 1)\n",
+    "ax.legend(handles=[Patch(facecolor=colors[i], label=name)\n",
+    "                   for i, name in enumerate(names)],\n",
+    "          loc='upper right', fontsize=10, framealpha=0.9)\n",
+    "ax.set_title('Rasterized by string label')\n",
+    "ax.set_axis_off()\n",
+    "plt.tight_layout()"
+   ],
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cat3482alert",
+   "metadata": {},
+   "source": [
+    "<div class=\"alert alert-block alert-info\">\n",
+    "<b>Labels travel to QGIS.</b> Writing this result with <code>to_geotiff(landcover, 'landcover.tif')</code> also writes a <code>landcover.tif.aux.xml</code> sidecar holding the category names and colors. GDAL reads the sidecar, so the file opens in QGIS showing the class names instead of bare codes. Keep the sidecar next to the <code>.tif</code> when you move the file. <code>open_geotiff</code> restores <code>category_names</code> / <code>category_colors</code> back onto <code>attrs</code>.\n",
+    "</div>"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "wtyzoml2vc",
@@ -596,4 +649,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}

xrspatial/geotiff/__init__.py

@@ -910,6 +910,17 @@ def open_geotiff(source: str | BinaryIO, *,
 
     source = _coerce_path(source)
 
+    def _attach_category_attrs(da):
+        # Categorical labels/colors live in a PAM ``<file>.aux.xml`` sidecar
+        # (GDAL ignores them embedded in the TIFF). Merge them back onto the
+        # result so a rasterize -> to_geotiff -> open_geotiff round-trip
+        # preserves ``category_names`` / ``category_colors``. Local string
+        # sources only; a missing/malformed sidecar yields {} and is a no-op.
+        if isinstance(source, str):
+            from ._pam import read_pam_sidecar
+            da.attrs.update(read_pam_sidecar(source))
+        return da
+
     # Resolve the rioxarray-compatible renames. ``masked`` / ``default_name``
     # are the canonical names; ``mask_nodata`` / ``name`` are deprecated
     # aliases kept for back-compat. Mirrors the sentinel-based deprecation in
@@ -1115,7 +1126,8 @@ def open_geotiff(source: str | BinaryIO, *,
         vrt_kwargs = {}
         if missing_sources_passed:
             vrt_kwargs['missing_sources'] = missing_sources
-        return _read_vrt(source, dtype=dtype, window=window, band=band,
+        return _attach_category_attrs(_read_vrt(
+                         source, dtype=dtype, window=window, band=band,
                          name=default_name, chunks=chunks, gpu=gpu,
                          max_pixels=max_pixels,
                          allow_rotated=allow_rotated,
@@ -1126,7 +1138,7 @@ def open_geotiff(source: str | BinaryIO, *,
                          allow_internal_only_jpeg=allow_internal_only_jpeg,
                          band_nodata=band_nodata,
                          mask_nodata=masked,
-                         **vrt_kwargs)
+                         **vrt_kwargs))
 
     # File-like buffer rejections for ``gpu=True`` / ``chunks=...`` already
     # fired inside ``_validate_dispatch_kwargs`` above; the non-VRT branches
@@ -1138,7 +1150,8 @@ def open_geotiff(source: str | BinaryIO, *,
         gpu_kwargs = {}
         if on_gpu_failure is not _ON_GPU_FAILURE_SENTINEL:
             gpu_kwargs['on_gpu_failure'] = on_gpu_failure
-        return _read_geotiff_gpu(source, dtype=dtype,
+        return _attach_category_attrs(_read_geotiff_gpu(
+                                 source, dtype=dtype,
                                  overview_level=overview_level,
                                  window=window, band=band,
                                  name=default_name, chunks=chunks,
@@ -1153,11 +1166,12 @@ def open_geotiff(source: str | BinaryIO, *,
                                     allow_internal_only_jpeg),
                                  mask_nodata=masked,
                                  mask_and_scale=unpack,
-                                 **gpu_kwargs)
+                                 **gpu_kwargs))
 
     # Dask path (CPU)
     if chunks is not None:
-        return _read_geotiff_dask(source, dtype=dtype, chunks=chunks,
+        return _attach_category_attrs(_read_geotiff_dask(
+                                  source, dtype=dtype, chunks=chunks,
                                   overview_level=overview_level,
                                   window=window, band=band,
                                   max_pixels=max_pixels, name=default_name,
@@ -1171,7 +1185,7 @@ def open_geotiff(source: str | BinaryIO, *,
                                      allow_internal_only_jpeg),
                                   mask_nodata=masked,
                                   mask_and_scale=unpack,
-                                  parse_coordinates=parse_coordinates)
+                                  parse_coordinates=parse_coordinates))
 
     kwargs = {}
     if max_pixels is not None:
@@ -1212,7 +1226,7 @@ def open_geotiff(source: str | BinaryIO, *,
         getattr(geo_info, '_mask_nodata', nodata)
         if nodata is not None else None
     )
-    return _finalize_eager_read(
+    return _attach_category_attrs(_finalize_eager_read(
         arr,
         geo_info=geo_info,
         nodata=nodata,
@@ -1226,7 +1240,7 @@ def open_geotiff(source: str | BinaryIO, *,
         mask_and_scale=unpack,
         parse_coordinates=parse_coordinates,
         band=band,
-    )
+    ))
 
 
 def plot_geotiff(da: xr.DataArray, **kwargs):

xrspatial/geotiff/_attrs.py

@@ -116,6 +116,12 @@
 - ``extra_samples``: TIFF ExtraSamples tag.
 - ``colormap``: raw uint16 RGB triples from the TIFF ColorMap tag (320),
   attached to single-band paletted images.
+- ``category_names``: ordered list of class label strings (index == pixel
+  value) for a categorical raster. Written to / read from a PAM
+  ``<file>.aux.xml`` sidecar (``<CategoryNames>`` plus a thematic
+  ``<GDALRasterAttributeTable>``); see :mod:`xrspatial.geotiff._pam`.
+- ``category_colors``: list of ``(r, g, b, a)`` int tuples (0-255), one per
+  category, emitted as the RAT's Red/Green/Blue/Alpha columns.
 
 Removed in contract v2 (issue #2016):