Computing Dwelling-Unit Density per Acre in Python

Density caps decide how many homes a parcel may legally hold, yet dwelling-unit counts and land areas usually live in separate layers measured in mismatched units. This guide walks through aggregating dwelling units onto parcels or a regular grid, computing dwelling units per acre inside a metric CRS, and flagging every cell that breaches its zoning density ceiling. The result is a reproducible DU/acre surface that planners can defend line by line during entitlement review.

Prerequisites

Step-by-step

Step 1: Project to a metric CRS and compute area

Acreage must come from a linear CRS, never from degrees. Reproject the parcels to a state plane or UTM system, derive area in square meters, then convert to acres with the exact factor. One acre equals 4046.8564224 square meters.

import geopandas as gpd

SQM_PER_ACRE = 4046.8564224

# State Plane California Zone 3 (feet) reprojected to meters via UTM 10N
parcels = gpd.read_file("parcels.gpkg").to_crs("EPSG:32610")

# Area is only valid in a projected CRS; degrees would distort it badly
parcels["area_m2"] = parcels.geometry.area
parcels["area_acre"] = parcels["area_m2"] / SQM_PER_ACRE

Step 2: Aggregate dwelling units to the analysis unit

If dwelling counts already sit on the parcels, this step is a passthrough. When units arrive as building points, spatially join them to parcels and sum. The join uses GeoPandas’ built-in spatial index, so it stays fast even on large point sets.

import pandas as pd

units = gpd.read_file("dwelling_points.gpkg").to_crs("EPSG:32610")

# Assign each dwelling point to the parcel that contains it
joined = gpd.sjoin(units, parcels[["parcel_id", "geometry"]],
                   how="inner", predicate="within")

# Sum dwelling units per parcel; parcels with none get zero
du_by_parcel = joined.groupby("parcel_id")["units"].sum()
parcels = parcels.merge(du_by_parcel.rename("dwelling_units"),
                        on="parcel_id", how="left")
parcels["dwelling_units"] = parcels["dwelling_units"].fillna(0)

Step 3: Optionally bin into a regular density grid

Parcel-level density is spiky; a fishnet grid smooths it for regional comparison. Build a square grid across the study extent, then area-weight each parcel’s units into the cells it overlaps so partial coverage is counted fairly.

import numpy as np
from shapely.geometry import box

CELL = 200  # 200 m cells in the metric CRS
xmin, ymin, xmax, ymax = parcels.total_bounds
cols = np.arange(xmin, xmax + CELL, CELL)
rows = np.arange(ymin, ymax + CELL, CELL)
grid = gpd.GeoDataFrame(
    {"geometry": [box(x, y, x + CELL, y + CELL) for x in cols for y in rows]},
    crs=parcels.crs,
)
grid["cell_id"] = range(len(grid))

# Area-weight dwelling units into the cells each parcel touches
parts = gpd.overlay(parcels, grid, how="intersection")
parts["frac"] = parts.geometry.area / parts["area_m2"]
parts["du_share"] = parts["dwelling_units"] * parts["frac"]
grid = grid.merge(parts.groupby("cell_id")["du_share"].sum(),
                  on="cell_id", how="left")
grid["du_share"] = grid["du_share"].fillna(0)

Step 4: Compute DU per acre and flag against the cap

Divide dwelling units by acreage on whichever unit you chose, then attach the applicable zoning cap and compare. A representative-point join keeps each parcel tied to a single governing district.

# Density on the parcel unit
parcels["du_per_acre"] = parcels["dwelling_units"] / parcels["area_acre"]

# Attach the zoning cap via the parcel's representative point
zoning = gpd.read_file("zoning.gpkg").to_crs("EPSG:32610")
reps = parcels.copy()
reps["geometry"] = reps.geometry.representative_point()
capped = gpd.sjoin(reps, zoning[["max_du_acre", "geometry"]],
                   how="left", predicate="within")
parcels["max_du_acre"] = capped["max_du_acre"].values
parcels["over_cap"] = parcels["du_per_acre"] > parcels["max_du_acre"]

Verification

Spot-check that densities are physically plausible and that the grid conserves total dwelling units. The area-weighted shares distributed across cells must sum back to the original parcel total within floating-point tolerance.

total_parcel_du = parcels["dwelling_units"].sum()
total_grid_du = grid["du_share"].sum()
assert abs(total_parcel_du - total_grid_du) < 1.0  # units conserved
print("over-cap parcels:", int(parcels["over_cap"].sum()))
print("max density seen:", round(parcels["du_per_acre"].max(), 2), "DU/acre")

This DU/acre surface feeds naturally into the broader automated density calculation grids topic area, and the point-to-parcel aggregation reuses the tuning from optimizing spatial joins for 100k parcel datasets.

Common Pitfalls

  • Computing acreage in EPSG:4326. Degree-based area is nonsense at parcel scale; a lot near 45° latitude can be off by more than a third. Always project first.
  • Double-counting units on shared boundaries. A within join keeps each dwelling point in exactly one parcel, whereas intersects can assign boundary points twice. Use within or representative_point deduplication.
  • Ignoring net-versus-gross density. Some ordinances exclude rights-of-way and easements from the denominator. If the cap is a net-density figure, subtract those areas before dividing, or the flag will misfire.

Frequently Asked Questions

Should I measure density on parcels or on a grid?

Use parcels when the ordinance regulates each lot individually, since that is the legal unit of enforcement. Use a grid when you need a smooth regional picture for comprehensive-plan analysis or heatmaps. Many teams compute both: parcel density for compliance flags and grid density for visualization and trend reporting.

Why convert to acres instead of leaving density in units per square meter?

Zoning codes in the United States almost universally express caps as dwelling units per acre, so matching that unit avoids conversion errors during review. Compute area in square meters for accuracy, then divide by the exact square-meters-per-acre constant so the reported figure lines up with the ordinance text.

How do I handle a parcel that straddles two zoning districts?

Attach the cap using the parcel’s representative point so a single governing district is chosen deterministically, then escalate genuinely split parcels to manual review. Splitting a parcel’s density between two caps is rarely how the code reads, so a documented single-district assignment plus a review flag is the defensible default.

What CRS should I pick for a multi-county study area?

Choose one projected CRS that covers the whole extent with acceptable distortion, such as a state plane zone or a single UTM zone. Consult the EPSG registry to pick the zone whose central meridian sits nearest your study area, and reproject every input to it before any area math.

Conclusion

Dwelling-unit density becomes auditable once area, unit counts, and caps all live in the same metric frame. By projecting first, aggregating with spatial joins, optionally binning into a conserving grid, and comparing against district ceilings, you produce a DU/acre layer that flags overbuilt parcels and holds up under scrutiny.