“Dasymetric mapping may be defined as a kind of areal interpolation that uses ancillary (additional and related) data to aid in the areal interpolation process. Dasymetric mapping differs from choropleth mapping in that the boundaries of cartographic representation are not arbitrary but reflect the spatial distribution of the variable being mapped (Eicher and Brewer 2001).”
Dasymetric mapping is a cartographic technique developed in the early 20th century Russia (Petrov, 2012) and designed to help improve accuracy in interpolation we can use auxiliary information to mask out areas we know should not be used to distribute the source variable. Some of the first papers to introduce the technique in the spatial analysia and GI-Science literature are Goodchild & Lam (1980), Langford et al. (1991), and Langford & Unwin (1994) with subsequent important contributions from Eicher & Brewer (2001) and Mennis (2003). See Kim & Yao (2010) for a nice overview of the logic behind dasymetric mapping, and Comber & Zeng (2019) for a modern review of commonly-used techniques. For extensions including hybrid approaches, see Jia & Gaughan (2016), Langford (2006), Langford (2007) and Qiu & Cromley (2013).
One standard approach is to use additional data sourced from satellite imagery or remote sensing and use it to mask out uninhabited areas (Fisher & Langford, 1996; Langford et al., 1991; Langford & Unwin, 1994; Ruther et al., 2015). For example we can use raster data like the National Land Cover Database NLCD or OpenLandMap to mask out uninhabited land uses from the source data. To do so, we need to provide a path to the raster and a list of pixel values that are considered developed. The tobler package can accept any kind of raster data that can be read by rasterio, so you can provide your own, or download directly from NLCD linked above. Alternatively, you can read a compressed version of the NLCD hosted in the spatialucr quilt directly over S3.
For this analysis we will convert Census data into a regular hexagonal grid, using NLCD data selecting from the raster image cells classified as ‘developed land’ in various intensities. To begin, we will collect tract-level ACS data for Washington D.C., then generate a set of h3 hexagons that cover the surface of those tracts.
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/pyproj/crs/crs.py:1295: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj = self._crs.to_proj4(version=version)
D.C. Tract Boundaries vs. Regular Hexgrid
5.1 Simple Dasymetric
The easiest way to carry out a dasymetric interpolation is to use tobler’s, masked_area_interpolate function, which takes raster data as an ancillary mask and uses it to modify the input source polygons. This is a simple “binary dasymetric” technique, one of the most common types of dasymetric interpolation in the literature, which means we use the ancillary mask as a binary switch to define which parts of the source are valid for interpolation into the target. In the cartography and spatial analysis literatures, other dasymetric techniques have also been developed that use different classes present in the anciallary data to assign differential weights during the areal interpolation process, but those methods have vary levels of efficacy are not implemented in the tobler package (Mennis, 2003, 2009; Mennis & Hultgren, 2006).
The interpolation itself uses the masked_area_interpolate function, which works almost identically to the area_interpolate function introduced in the previous section, save that this function also accepts the path to a raster file (often a GeoTIFF format, but not necessarily), and a set of pixel values that represent areas considered to be inhabited; area of the census tracts covered by all other cell values will be masked out. In this example we interpolate the total population (which results in a density surface, given that we have consistently-sized zones), and select the NLCD values of 22, 23, and 24, corresponding to low, medium, and high intensity developed land uses.
Code
dc_hexes = masked_area_interpolate( source_df=dc, target_df=hexes, raster="https://spatial-ucr.s3.amazonaws.com/nlcd/landcover/nlcd_landcover_2019.tif", pixel_values=[22, 23, 24], # low, medium, and high intensity development extensive_variables=["n_total_pop"],)ax = dc_hexes.plot("n_total_pop", figsize=(3,3))ax.set_title("D.C. Population in a Regular Hexgrid")ax.axis("off")plt.show()
D.C. Population Interpolated to Hexgrid using Ancillary Information
5.1.1 Comparing Dasymetric with Area-Weighted
To get a sense for the dasymetric technique and how it compares (hopefully improves) the generic area-weighted interpolation approach, we can take the difference of both techniques, then plot a map of the difference laid over an arial image of the study area. This should give a sense for whether the uninhabited areas are being skirted by the interpolation process.
Code
# hexes using regular areal interpolation for comparisondc_hexes_simple = area_interpolate( source_df=dc, target_df=hexes, extensive_variables=["n_total_pop"])# take the difference between dasymetric and area-weighteddiff = dc_hexes.n_total_pop - dc_hexes_simple.n_total_pop# plot a histogram of the differences by hexdiff.plot(kind="hist")plt.show()
Histogram of Difference in Area-Weighted vs Dasymetric Population Estimates
Total Population Difference Between Area-Weighted and Dasymetric Techniques
The blue hexes have lower estimated population when we use the raster, and the red hexes have higher values. In downtown DC, the differences are small. But in certain places, like rock creek park in the northwest, and along the Anacostia river, the values are lower when using the raster as a dasymetric layer (showing it’s doing the correct thing). Instead of population being allocated to this green space (and water), it is instead being distributed to the developed areas nearby. If we plot the difference map on top of a satellite image, it’s easy to see most of the time the blue hexes are parks and forests where we do not want to allocate population, and the darkest red hexes are dense neighborhoods adjacent to the blue hexes. In these cases, we are using the raster data effectively to send data to the populated places instead of the uninhabited places. Although it’s not quite perfectly centered on white with this colormap, the Fisher Jenks classification scheme works nicely here.
An important thing to remember about most of the commonly-used dasymetric approaches (e.g. binary dasymetric) is their focus on the lower part of the distribution. That is, these techniques help ensure that population, etc is not distributed to places that are probably uninhabited, but they still assume a uniform distribution within each of the masked areas. Remotely-sensed land-use data often classifies high intensity land uses like industrial sites and shopping centers in the same category as high-rise apartment buildings, so dasymetric methods may over-allocate variables to developed but non-residential spaces. Again it is important to consider the variable under study, the quality of the source and target datasets, and the ancillary data used as a dasymetric mask
5.2 Advanced Dasymetric
In addition to the simple approach demonstrated above, we can also carry out a more complex estimation leveraging more than one ancillary data source, for example using satellite imagery and the presence of roads (meaning we will use both vector and raster data), albeit using a few more geoprocessing steps. The dasymetric data include impervious surface (from satellite imagery) and proximity to roads, each of which have been shown to improve estimation accuracy (Reibel & Bufalino, 2005), and in this case we require both (although there is no clear best answer to ‘which collection of data is best’ to use as a dasymetric mask (Zandbergen & Ignizio, 2010)). Thus we need to collect raster data, convert it to vector format, collect and buffer road data, then use the intersection of these two layers to define our dasymetric surface before interpolating source into target geometries.
NHGIS uses a two step process, in which they first apply a binary dasymetric filter, then carry out target density weighted interpolation (Schroeder, 2017). Following, we (nearly) replicate the NHGIS binary dasymetric filtering, albeit with slightly different open data and open source tooling, before carrying out a simpler areal interpolation (without the target density weighting, which is not possible with this target data (Schroeder, 2007)).
In NHGIS’s BD model for 2000 block data, the inhabited zone consists of all areas that are at least 5% developed impervious surface (within each 30-meter square cell of NLCD 2001 data) and lie within 300 feet of a residential road center line2 but not in a water body, using road definitions and water polygons from the 2010 TIGER/Line Shapefiles.
To begin, we collect data representing residential streets from OpenStreetMap, then buffer the line features out 300 feet.
Now we use the street buffers to extract any cells between 5 and 100% impervious surface, converting from a raster to a vector. This geodataframe represents the areas that meet both criteria NHGIS suggests. The remaining task is to cut our input census data, “restricting” it to the locations NLCD defines as 5-100% impervious. To do so, we just need to clip the census data using our extracted urban features as a mask. Using the intersection overlay, “clips” the census geometries by the ‘urban features’ leaving the original attributes intact but reshaping the geometries.
If a census tract had 2000 people in it and covered one square mile, but only a quarter mile of the tract was urbanized, the new geodataframe effectively shows these 2000 people occupying the quarter mile area instead of the original tract area. By passing these new data to the areal_interpolate function, we’re still assuming that population density is constant across our new source geometries (a condition that may not be true in reality) but that assumption is much more plausible than when using original census boundaries. Thus we have not changed the values associated with our original source, we have only refined the area of overlap between source and target.
This also means any combination of ancillary data can be combined to create (multiple) dasymetric mask(s). If you have a geodataframe representing water features, it can be used to “erase” overlapping portion of the source (exclusion) or if you have a layer of buildings, it could be used to select only the areas of the source covered by building footprints1 (inclusion) (if you have a crazy architectural marvel that straddles the water, you could use both!). The point here is it is best to consider the kinds of ancillary data available, and how to refine the source polygons before carrying out the area-weighting process; this is the essence of ‘dasymetric interpolation’.
