Context
zonal.stats aggregates raster values within zones, but there's nothing for the inverse: redistributing aggregate zone-level data (population per census tract, say) onto a finer-resolution raster using ancillary weights (land cover, nighttime lights, building footprints, etc).
This is usually called dasymetric mapping. The existing Python options (like tobler) are vector-based and don't work with dask or cupy.
Proposed module: xrspatial/dasymetric.py
Core function: disaggregate(zones, values, weight, method='binary')
zones -- integer raster of zone membership per pixel (e.g. rasterized census tracts)
values -- dict or Series mapping zone ID to the value to redistribute (e.g. {tract_id: population})
weight -- continuous raster of ancillary weights
method -- 'binary', 'weighted', or 'limiting_variable'
Returns a float raster where each pixel gets its share of the zone total, proportional to its weight.
Methods
'binary' -- weight treated as a mask (0/1), value split equally among non-zero pixels
'weighted' -- pixel_value = zone_value * (pixel_weight / zone_weight_sum)
'limiting_variable' -- three-class dasymetric: ancillary classes define density caps, iterative redistribution of residuals
Algorithm for 'weighted':
- For each zone z:
zone_sum[z] = sum(weight[zones == z])
- For each pixel:
result[i,j] = values[zones[i,j]] * weight[i,j] / zone_sum[zones[i,j]]
Grouped normalization + per-zone scalar multiply. The zonal sum can reuse zonal.stats internals or be computed directly.
Backend notes
- numpy/cupy: vectorized with
np.unique + fancy indexing
- dask: zones span chunk boundaries, so compute zonal sums eagerly (one float per zone, small), then distribute lazily per chunk via
map_blocks
- dask+cupy: same pattern with cupy chunk functions
Possible follow-ups (separate PRs)
pycnophylactic(zones, values) -- Tobler's smooth interpolation that preserves zone totals (iterative Laplacian smoothing with mass correction)
validate_disaggregation(result, zones, values) -- check that zone totals are preserved within tolerance
References
- Mennis (2003). Generating surface models of population using dasymetric mapping. The Professional Geographer, 55(1), 31-42.
- Tobler (1979). Smooth pycnophylactic interpolation for geographical regions. JASA, 74(367), 519-530.
Context
zonal.statsaggregates raster values within zones, but there's nothing for the inverse: redistributing aggregate zone-level data (population per census tract, say) onto a finer-resolution raster using ancillary weights (land cover, nighttime lights, building footprints, etc).This is usually called dasymetric mapping. The existing Python options (like
tobler) are vector-based and don't work with dask or cupy.Proposed module:
xrspatial/dasymetric.pyCore function:
disaggregate(zones, values, weight, method='binary')zones-- integer raster of zone membership per pixel (e.g. rasterized census tracts)values-- dict or Series mapping zone ID to the value to redistribute (e.g.{tract_id: population})weight-- continuous raster of ancillary weightsmethod--'binary','weighted', or'limiting_variable'Returns a float raster where each pixel gets its share of the zone total, proportional to its weight.
Methods
'binary'-- weight treated as a mask (0/1), value split equally among non-zero pixels'weighted'--pixel_value = zone_value * (pixel_weight / zone_weight_sum)'limiting_variable'-- three-class dasymetric: ancillary classes define density caps, iterative redistribution of residualsAlgorithm for
'weighted':zone_sum[z] = sum(weight[zones == z])result[i,j] = values[zones[i,j]] * weight[i,j] / zone_sum[zones[i,j]]Grouped normalization + per-zone scalar multiply. The zonal sum can reuse
zonal.statsinternals or be computed directly.Backend notes
np.unique+ fancy indexingmap_blocksPossible follow-ups (separate PRs)
pycnophylactic(zones, values)-- Tobler's smooth interpolation that preserves zone totals (iterative Laplacian smoothing with mass correction)validate_disaggregation(result, zones, values)-- check that zone totals are preserved within toleranceReferences