Skip to content

Commit 27615ba

Browse files
brendancolbrendancol
committed
Add raster-based dasymetric mapping module (#930)
Adds disaggregate(), pycnophylactic(), and validate_disaggregation() for redistributing zone-level values onto finer raster grids. disaggregate supports three methods (binary, weighted, limiting_variable) across all four backends (numpy, cupy, dask+numpy, dask+cupy). pycnophylactic implements Tobler's smooth interpolation with iterative Laplacian smoothing and mass correction (numpy/cupy). validate_disaggregation checks that zone totals are preserved within tolerance across all backends. 51 tests covering known values, conservation, edge cases, cross-backend parity, and input validation.
1 parent 25309d8 commit 27615ba

File tree

8 files changed

+1584
-0
lines changed

8 files changed

+1584
-0
lines changed

‎CHANGELOG.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
### Version 0.8.0 - 2026-03-03
66

77
#### New Features
8+
- Add raster-based dasymetric mapping: disaggregate, pycnophylactic interpolation, and validation helper (#930)
89
- Add enhanced terrain generation features (#929)
910
- Add fire module: dNBR, RdNBR, burn severity, fireline intensity, flame length, rate of spread (Rothermel), and KBDI (#927)
1011
- Add flood prediction tools (#926)

‎docs/source/reference/dasymetric.rst‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
.. _dasymetric:
2+
3+
***********
4+
Dasymetric
5+
***********
6+
7+
Disaggregate
8+
============
9+
.. autosummary::
10+
:toctree: _autosummary
11+
12+
xrspatial.dasymetric.disaggregate
13+
14+
Pycnophylactic
15+
==============
16+
.. autosummary::
17+
:toctree: _autosummary
18+
19+
xrspatial.dasymetric.pycnophylactic
20+
21+
Validate
22+
========
23+
.. autosummary::
24+
:toctree: _autosummary
25+
26+
xrspatial.dasymetric.validate_disaggregation

‎docs/source/reference/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ Reference
88
:maxdepth: 2
99

1010
classification
11+
dasymetric
1112
fire
1213
flood
1314
focal

‎examples/user_guide/14_Dasymetric_Mapping.ipynb‎

Lines changed: 216 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,216 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Dasymetric Mapping\n",
8+
"\n",
9+
"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",
10+
"\n",
11+
"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",
12+
"\n",
13+
"`xrspatial.disaggregate` supports three methods:\n",
14+
"- **`'binary'`** -- split value equally among nonzero-weight pixels\n",
15+
"- **`'weighted'`** (default) -- distribute proportionally to weight values\n",
16+
"- **`'limiting_variable'`** -- three-class dasymetric with density caps (numpy-only)"
17+
]
18+
},
19+
{
20+
"cell_type": "code",
21+
"execution_count": null,
22+
"metadata": {},
23+
"outputs": [],
24+
"source": [
25+
"import numpy as np\n",
26+
"import xarray as xr\n",
27+
"import matplotlib.pyplot as plt\n",
28+
"\n",
29+
"from xrspatial import disaggregate"
30+
]
31+
},
32+
{
33+
"cell_type": "markdown",
34+
"metadata": {},
35+
"source": [
36+
"## Generate Synthetic Data\n",
37+
"\n",
38+
"We create a small grid of census-tract-like zones, assign population values per zone, and build a land-cover-based weight raster."
39+
]
40+
},
41+
{
42+
"cell_type": "code",
43+
"execution_count": null,
44+
"metadata": {},
45+
"outputs": [],
46+
"source": [
47+
"# 10x10 zone raster with 4 \"census tracts\"\n",
48+
"zones_data = np.ones((10, 10), dtype=np.float64)\n",
49+
"zones_data[:5, :5] = 1\n",
50+
"zones_data[:5, 5:] = 2\n",
51+
"zones_data[5:, :5] = 3\n",
52+
"zones_data[5:, 5:] = 4\n",
53+
"\n",
54+
"zones = xr.DataArray(zones_data, dims=['y', 'x'])\n",
55+
"\n",
56+
"# Population per tract\n",
57+
"population = {1: 1000.0, 2: 500.0, 3: 2000.0, 4: 300.0}\n",
58+
"\n",
59+
"# Ancillary weight raster -- simulates land cover suitability\n",
60+
"np.random.seed(42)\n",
61+
"weight_data = np.random.rand(10, 10).astype(np.float64)\n",
62+
"# Make some areas uninhabitable (zero weight)\n",
63+
"weight_data[0:2, 0:2] = 0.0 # water body in tract 1\n",
64+
"weight_data[7:9, 7:9] = 0.0 # park in tract 4\n",
65+
"\n",
66+
"weight = xr.DataArray(weight_data, dims=['y', 'x'])\n",
67+
"\n",
68+
"fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
69+
"zones.plot(ax=axes[0], cmap='Set2')\n",
70+
"axes[0].set_title('Zones (Census Tracts)')\n",
71+
"weight.plot(ax=axes[1], cmap='YlGn')\n",
72+
"axes[1].set_title('Weight (Land Cover)')\n",
73+
"plt.tight_layout()\n",
74+
"plt.show()"
75+
]
76+
},
77+
{
78+
"cell_type": "markdown",
79+
"metadata": {},
80+
"source": [
81+
"## Binary Method\n",
82+
"\n",
83+
"The `'binary'` method binarizes the weight raster (nonzero becomes 1, zero stays 0) and splits each zone's value equally among its nonzero pixels."
84+
]
85+
},
86+
{
87+
"cell_type": "code",
88+
"execution_count": null,
89+
"metadata": {},
90+
"outputs": [],
91+
"source": [
92+
"result_binary = disaggregate(zones, population, weight, method='binary')\n",
93+
"\n",
94+
"fig, ax = plt.subplots(figsize=(6, 5))\n",
95+
"result_binary.plot(ax=ax, cmap='YlOrRd')\n",
96+
"ax.set_title('Binary Dasymetric Result')\n",
97+
"plt.tight_layout()\n",
98+
"plt.show()"
99+
]
100+
},
101+
{
102+
"cell_type": "markdown",
103+
"metadata": {},
104+
"source": [
105+
"## Weighted Method\n",
106+
"\n",
107+
"The `'weighted'` method (default) distributes each zone's value proportionally to the weight values:\n",
108+
"\n",
109+
"$$\\text{pixel} = \\text{zone\\_value} \\times \\frac{\\text{pixel\\_weight}}{\\sum_{\\text{zone}} \\text{weights}}$$"
110+
]
111+
},
112+
{
113+
"cell_type": "code",
114+
"execution_count": null,
115+
"metadata": {},
116+
"outputs": [],
117+
"source": [
118+
"result_weighted = disaggregate(zones, population, weight, method='weighted')\n",
119+
"\n",
120+
"fig, ax = plt.subplots(figsize=(6, 5))\n",
121+
"result_weighted.plot(ax=ax, cmap='YlOrRd')\n",
122+
"ax.set_title('Weighted Dasymetric Result')\n",
123+
"plt.tight_layout()\n",
124+
"plt.show()"
125+
]
126+
},
127+
{
128+
"cell_type": "markdown",
129+
"metadata": {},
130+
"source": [
131+
"## Conservation Property\n",
132+
"\n",
133+
"A key property of dasymetric mapping: the sum of output pixel values within each zone should equal the original zone total."
134+
]
135+
},
136+
{
137+
"cell_type": "code",
138+
"execution_count": null,
139+
"metadata": {},
140+
"outputs": [],
141+
"source": [
142+
"for zid, expected in population.items():\n",
143+
" actual = float(np.nansum(result_weighted.values[zones_data == zid]))\n",
144+
" print(f'Zone {zid}: expected={expected:.1f}, actual={actual:.1f}, '\n",
145+
" f'diff={abs(expected - actual):.2e}')"
146+
]
147+
},
148+
{
149+
"cell_type": "markdown",
150+
"metadata": {},
151+
"source": [
152+
"## Comparing Methods\n",
153+
"\n",
154+
"The binary method produces uniform density within each zone (among habitable pixels), while the weighted method produces spatially varying density."
155+
]
156+
},
157+
{
158+
"cell_type": "code",
159+
"execution_count": null,
160+
"metadata": {},
161+
"outputs": [],
162+
"source": [
163+
"fig, axes = plt.subplots(1, 2, figsize=(12, 5))\n",
164+
"\n",
165+
"result_binary.plot(ax=axes[0], cmap='YlOrRd', vmin=0,\n",
166+
" vmax=float(np.nanmax(result_weighted.values)))\n",
167+
"axes[0].set_title('Binary')\n",
168+
"\n",
169+
"result_weighted.plot(ax=axes[1], cmap='YlOrRd', vmin=0,\n",
170+
" vmax=float(np.nanmax(result_weighted.values)))\n",
171+
"axes[1].set_title('Weighted')\n",
172+
"\n",
173+
"plt.tight_layout()\n",
174+
"plt.show()"
175+
]
176+
},
177+
{
178+
"cell_type": "markdown",
179+
"metadata": {},
180+
"source": [
181+
"## Limiting Variable Method\n",
182+
"\n",
183+
"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."
184+
]
185+
},
186+
{
187+
"cell_type": "code",
188+
"execution_count": null,
189+
"metadata": {},
190+
"outputs": [],
191+
"source": [
192+
"result_lv = disaggregate(zones, population, weight,\n",
193+
" method='limiting_variable')\n",
194+
"\n",
195+
"fig, ax = plt.subplots(figsize=(6, 5))\n",
196+
"result_lv.plot(ax=ax, cmap='YlOrRd')\n",
197+
"ax.set_title('Limiting Variable Result')\n",
198+
"plt.tight_layout()\n",
199+
"plt.show()"
200+
]
201+
}
202+
],
203+
"metadata": {
204+
"kernelspec": {
205+
"display_name": "Python 3",
206+
"language": "python",
207+
"name": "python3"
208+
},
209+
"language_info": {
210+
"name": "python",
211+
"version": "3.11.0"
212+
}
213+
},
214+
"nbformat": 4,
215+
"nbformat_minor": 4
216+
}

‎xrspatial/__init__.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,9 @@
11
from xrspatial.aspect import aspect # noqa
22
from xrspatial.bump import bump # noqa
33
from xrspatial.cost_distance import cost_distance # noqa
4+
from xrspatial.dasymetric import disaggregate # noqa
5+
from xrspatial.dasymetric import pycnophylactic # noqa
6+
from xrspatial.dasymetric import validate_disaggregation # noqa
47
from xrspatial.classify import binary # noqa
58
from xrspatial.classify import box_plot # noqa
69
from xrspatial.classify import head_tail_breaks # noqa

‎xrspatial/accessor.py‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -235,6 +235,16 @@ def regions(self, **kwargs):
235235
from .zonal import regions
236236
return regions(self._obj, **kwargs)
237237

238+
# ---- Dasymetric ----
239+
240+
def disaggregate(self, values, weight, **kwargs):
241+
from .dasymetric import disaggregate
242+
return disaggregate(self._obj, values, weight, **kwargs)
243+
244+
def pycnophylactic(self, values, **kwargs):
245+
from .dasymetric import pycnophylactic
246+
return pycnophylactic(self._obj, values, **kwargs)
247+
238248
# ---- Terrain generation ----
239249

240250
def generate_terrain(self, **kwargs):

0 commit comments

Comments
 (0)