Add raster-based dasymetric mapping module

CHANGELOG.md

@@ -5,6 +5,7 @@
 ### Version 0.8.0 - 2026-03-03
 
 #### New Features
+- Add raster-based dasymetric mapping: disaggregate, pycnophylactic interpolation, and validation helper (#930)
 - Add enhanced terrain generation features (#929)
 - Add fire module: dNBR, RdNBR, burn severity, fireline intensity, flame length, rate of spread (Rothermel), and KBDI (#927)
 - Add flood prediction tools (#926)

docs/source/reference/dasymetric.rst

@@ -0,0 +1,26 @@
+..  _dasymetric:
+
+***********
+Dasymetric
+***********
+
+Disaggregate
+============
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.dasymetric.disaggregate
+
+Pycnophylactic
+==============
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.dasymetric.pycnophylactic
+
+Validate
+========
+.. autosummary::
+    :toctree: _autosummary
+
+    xrspatial.dasymetric.validate_disaggregation

docs/source/reference/index.rst

@@ -8,6 +8,7 @@ Reference
    :maxdepth: 2
 
    classification
+   dasymetric
    fire
    flood
    focal

examples/user_guide/14_Dasymetric_Mapping.ipynb

@@ -0,0 +1,216 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Dasymetric Mapping\n",
+    "\n",
+    "Dasymetric mapping redistributes aggregate zone-level data (like population per census tract) onto a finer raster grid using ancillary weight information (land cover, nighttime lights, etc.).\n",
+    "\n",
+    "This is the spatial inverse of `zonal.stats`: instead of aggregating pixel values into zone summaries, we spread zone-level totals back across pixels, weighted by auxiliary data.\n",
+    "\n",
+    "`xrspatial.disaggregate` supports three methods:\n",
+    "- **`'binary'`** -- split value equally among nonzero-weight pixels\n",
+    "- **`'weighted'`** (default) -- distribute proportionally to weight values\n",
+    "- **`'limiting_variable'`** -- three-class dasymetric with density caps (numpy-only)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "from xrspatial import disaggregate"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generate Synthetic Data\n",
+    "\n",
+    "We create a small grid of census-tract-like zones, assign population values per zone, and build a land-cover-based weight raster."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 10x10 zone raster with 4 \"census tracts\"\n",
+    "zones_data = np.ones((10, 10), dtype=np.float64)\n",
+    "zones_data[:5, :5] = 1\n",
+    "zones_data[:5, 5:] = 2\n",
+    "zones_data[5:, :5] = 3\n",
+    "zones_data[5:, 5:] = 4\n",
+    "\n",
+    "zones = xr.DataArray(zones_data, dims=['y', 'x'])\n",
+    "\n",
+    "# Population per tract\n",
+    "population = {1: 1000.0, 2: 500.0, 3: 2000.0, 4: 300.0}\n",
+    "\n",
+    "# Ancillary weight raster -- simulates land cover suitability\n",
+    "np.random.seed(42)\n",
+    "weight_data = np.random.rand(10, 10).astype(np.float64)\n",
+    "# Make some areas uninhabitable (zero weight)\n",
+    "weight_data[0:2, 0:2] = 0.0  # water body in tract 1\n",
+    "weight_data[7:9, 7:9] = 0.0  # park in tract 4\n",
+    "\n",
+    "weight = xr.DataArray(weight_data, dims=['y', 'x'])\n",
+    "\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
+    "zones.plot(ax=axes[0], cmap='Set2')\n",
+    "axes[0].set_title('Zones (Census Tracts)')\n",
+    "weight.plot(ax=axes[1], cmap='YlGn')\n",
+    "axes[1].set_title('Weight (Land Cover)')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Binary Method\n",
+    "\n",
+    "The `'binary'` method binarizes the weight raster (nonzero becomes 1, zero stays 0) and splits each zone's value equally among its nonzero pixels."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "result_binary = disaggregate(zones, population, weight, method='binary')\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 5))\n",
+    "result_binary.plot(ax=ax, cmap='YlOrRd')\n",
+    "ax.set_title('Binary Dasymetric Result')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Weighted Method\n",
+    "\n",
+    "The `'weighted'` method (default) distributes each zone's value proportionally to the weight values:\n",
+    "\n",
+    "$$\\text{pixel} = \\text{zone\\_value} \\times \\frac{\\text{pixel\\_weight}}{\\sum_{\\text{zone}} \\text{weights}}$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "result_weighted = disaggregate(zones, population, weight, method='weighted')\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 5))\n",
+    "result_weighted.plot(ax=ax, cmap='YlOrRd')\n",
+    "ax.set_title('Weighted Dasymetric Result')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Conservation Property\n",
+    "\n",
+    "A key property of dasymetric mapping: the sum of output pixel values within each zone should equal the original zone total."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for zid, expected in population.items():\n",
+    "    actual = float(np.nansum(result_weighted.values[zones_data == zid]))\n",
+    "    print(f'Zone {zid}: expected={expected:.1f}, actual={actual:.1f}, '\n",
+    "          f'diff={abs(expected - actual):.2e}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Comparing Methods\n",
+    "\n",
+    "The binary method produces uniform density within each zone (among habitable pixels), while the weighted method produces spatially varying density."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
+    "\n",
+    "result_binary.plot(ax=axes[0], cmap='YlOrRd', vmin=0,\n",
+    "                   vmax=float(np.nanmax(result_weighted.values)))\n",
+    "axes[0].set_title('Binary')\n",
+    "\n",
+    "result_weighted.plot(ax=axes[1], cmap='YlOrRd', vmin=0,\n",
+    "                     vmax=float(np.nanmax(result_weighted.values)))\n",
+    "axes[1].set_title('Weighted')\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Limiting Variable Method\n",
+    "\n",
+    "The `'limiting_variable'` method applies per-class density caps and redistributes overflow iteratively. This is useful when you know certain land cover types have maximum population densities. Currently only available for numpy arrays."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "result_lv = disaggregate(zones, population, weight,\n",
+    "                         method='limiting_variable')\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 5))\n",
+    "result_lv.plot(ax=ax, cmap='YlOrRd')\n",
+    "ax.set_title('Limiting Variable Result')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

xrspatial/__init__.py

@@ -1,6 +1,9 @@
 from xrspatial.aspect import aspect  # noqa
 from xrspatial.bump import bump  # noqa
 from xrspatial.cost_distance import cost_distance  # noqa
+from xrspatial.dasymetric import disaggregate  # noqa
+from xrspatial.dasymetric import pycnophylactic  # noqa
+from xrspatial.dasymetric import validate_disaggregation  # noqa
 from xrspatial.classify import binary  # noqa
 from xrspatial.classify import box_plot  # noqa
 from xrspatial.classify import head_tail_breaks  # noqa

xrspatial/accessor.py

@@ -235,6 +235,16 @@ def regions(self, **kwargs):
         from .zonal import regions
         return regions(self._obj, **kwargs)
 
+    # ---- Dasymetric ----
+
+    def disaggregate(self, values, weight, **kwargs):
+        from .dasymetric import disaggregate
+        return disaggregate(self._obj, values, weight, **kwargs)
+
+    def pycnophylactic(self, values, **kwargs):
+        from .dasymetric import pycnophylactic
+        return pycnophylactic(self._obj, values, **kwargs)
+
     # ---- Terrain generation ----
 
     def generate_terrain(self, **kwargs):