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.

Choosing the Area Basis Before You Divide

The division is the easy part. Which area goes underneath it is the decision that determines whether the answer matches the one a planner would compute by hand.

Gross and net area, and what separates themGross parcel area, the exclusions a code may name, and the resulting net area that a density figure divides by.Included in the denominatorEffect on the densityGross parcel areaEverything inside the recorded boundaryLowest reported densityLess internal rights-of-wayExcluded where the code says soRaises density; common on subdivisionsLess recorded easementsAccess, utility, drainageRaises density on encumbered parcelsLess undevelopable landWater, steep slope, where namedLargest effect on constrained sites
Same units, two denominators, two defensible densities. The code says which.

Gross parcel area is the recorded parcel, and it is what most quick calculations use. Net area subtracts what the code says to subtract — rights-of-way inside the parcel, recorded easements, sometimes land under water or above a slope threshold — and on parcels with any of those it produces a materially higher density from the same unit count. A subdivision proposal evaluated gross will look compliant where the same proposal evaluated net does not, and the code will have said which one applies.

def du_per_acre(units: int, parcel, exclusions=None, basis="net"):
    """Density with the basis stated. Returns the components, not just the ratio."""
    gross_sqft = parcel.area
    excluded_sqft = 0.0
    if basis == "net" and exclusions is not None and not exclusions.empty:
        excluded_sqft = exclusions.geometry.union_all().intersection(parcel).area
    site_sqft = gross_sqft - excluded_sqft
    if site_sqft <= 0:
        return {"outcome": "indeterminate", "reason": "no net site area after exclusions"}
    acres = site_sqft / 43560.0
    return {"units": units, "basis": basis, "gross_sqft": round(gross_sqft, 1),
            "excluded_sqft": round(excluded_sqft, 1), "acres": round(acres, 4),
            "du_per_acre": round(units / acres, 3)}

Returning the components rather than a bare number is what makes the result checkable, and it costs nothing since every value is already in hand.

Counting Units Without Double-Counting

The numerator carries its own ambiguity, and it is the source of most disagreements between a pipeline’s figure and a planner’s.

Per parcel, or across the assemblage?A development spanning several parcels under one application can be measured per parcel or across the combined site, and the two give different verdicts.Does the applicationcombine several parcelsinto one development site?single parcelPer-parcel density is the whole answerthe ordinary caseassemblageCompute both, per parcel and combinedrequires the application id, not just the parcel fabricRecord which basis governed and why
The pipeline should compute both; only the code can say which governs.

Three questions settle it. What counts as a unit — a dwelling, an accessory dwelling unit, a bedroom in a code that regulates by occupancy? Are existing units included alongside proposed ones, or is the cap on the net addition? And where a development spans several parcels under one application, is the density computed per parcel or across the assemblage?

That last one produces the most surprising results. A project on three parcels with all its units on one of them will fail a per-parcel density test and pass an assemblage test, and both are defensible readings depending on what the code says about combined lots. The pipeline’s job is to compute both and let the rule select, which requires knowing that the three parcels belong to one application — an attribute that comes from the permit system rather than from the parcel fabric.

Verifying Against a Known Case

The fastest way to confirm a density calculation is to reproduce a figure someone has already computed by hand, and to reconcile every component rather than the final number.

Reconciling against an approved applicationUnit count, gross area, excluded area and acres are each compared against a hand-computed figure before the density itself is compared.Unit countfrom the applicationGross areavs assessor acreageExcluded areaby categoryDensityonly now compared
A matching density from two cancelling errors will stop matching on the next parcel. Reconcile the components.

Take an approved application with a stated density, and compare the unit count, the gross area, the excluded area and the resulting acres before comparing the density itself. A matching density built from a unit count that is one too high and an area that is correspondingly large is two errors cancelling, and it will stop cancelling on the next parcel. Where the site area matches the assessor’s recorded acreage, that also confirms the working frame and its units are right, which is a second check for free.

Frequently Asked Questions

Should density be rounded before comparison against the limit?

No. Compare the full-precision value against the threshold with the stated tolerance, and round only for display. Rounding first turns 14.04 into 14.0 and can convert a marginal violation into a pass, which is a rounding convention making a regulatory decision.

What if the parcel has no units yet?

Then its density is zero and it is compliant with a maximum-density rule, which is correct but rarely the question being asked. Density checks are normally run against a proposal, so the unit count should come from the application rather than from the current record — and where it comes from should be recorded with the verdict.

How do accessory dwelling units affect the count?

That depends entirely on the code, and codes have changed on this repeatedly in recent years. Many now exempt them from density calculations up to a limit. Because it is both jurisdiction-specific and recently amended, it belongs in the rule record with an effective date rather than in the counting code.

Part of: Automated density calculation grids

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.