How to Normalise Census Counts into Rates and Densities

Problem statement

Turning a census count into a rate or a density is a single division. Each half of that division has a way to go wrong without an error, and all four showed up in the 2023 ACS tract tables:

  • Sentinel values in the numerator. 1,547 tract median incomes are the code โˆ’666666666. The mean of the column as downloaded is โˆ’11,997,152 dollars; with the codes removed it is 83,552.
  • Zero denominators. 1,040 tracts have a poverty universe of zero, so their rate is undefined, not zero.
  • The wrong population. Dividing people below poverty by total population instead of the poverty universe moved 1,378 tract rates by more than five percentage points, and the two versions agreed on only 52 of the top 100 tracts.
  • The wrong area. Dividing by geometry.area in longitude and latitude gives people per square degree, which shifted 2,079 tracts by more than five percentile points โ€” despite a rank correlation of 0.997 that makes it look harmless.

Quick answer

import pandas as pd

SENTINELS = [-222222222, -333333333, -555555555, -666666666, -888888888, -999999999]

def clean(series):
    return series.where(~series.isin(SENTINELS))

tracts["poverty_rate"] = clean(tracts["below"]) / clean(tracts["universe"]).where(tracts["universe"] > 0)
tracts["density_km2"] = tracts["population"] / (tracts["ALAND"] / 1e6)    # ALAND is square metres

Three rules sit behind those two lines: remove sentinels before any arithmetic, divide by the population the table was measured on, and divide by land area in square metres rather than by the shape of the polygon.

Triage table of four silent errors when normalising census counts โ€” sentinel values, zero denominators, the wrong population and the wrong area โ€” with the fix for each.
None of the four raises an exception; all four change what the map shows.

Step-by-step solution

1. Remove sentinel values first

ACS tables encode "not available" with large negative numbers: โˆ’666666666 when an estimate could not be computed, โˆ’222222222 and โˆ’333333333 in margins, and others. They are valid integers, so every aggregate happily includes them:

income = tracts["B19013_E001"]
print(round(income.mean()), round(income.where(income > 0).mean()))
-11997152 83552

Replace the codes with missing values when the table is loaded, and never later โ€” a rate or density built on a sentinel is just a very large negative number.

2. Divide by the table's own universe

Every ACS table is measured on a universe, named in the table shell and published as its first row. For poverty status, B17001_E001 is "Population for whom poverty status is determined", which excludes people in institutions, college dormitories and military barracks.

Across tracts the universe is usually close to total population, with a median of 99.7%. In 623 tracts it is less than half, and those are the tracts where the denominator choice decides the answer: one Miami-Dade tract of 3,737 residents has a poverty universe of 3.

3. Leave zero denominators undefined

Division by zero in pandas does not raise. Zero over zero gives NaN; a positive count over zero gives inf, which sorts to the top of a ranking and breaks colour classes. Mask the denominator first:

rate = numerator / denominator.where(denominator > 0)

1,040 tracts have a poverty universe of zero and 849 have no residents at all. Show them as "no data", which is what they are.

4. Use total population only for totals

"Per capita" is the right normalisation for totals that belong to everyone, such as aggregate household income or total vehicles. For a characteristic that is only measured on part of the population โ€” poverty, school enrolment, employment โ€” use that part. Measured on tract poverty, the wrong choice moved 1,378 rates by more than five points.

5. Build densities from land area

The Census Bureau publishes land and water area for every geography as ALAND and AWATER, in square metres, in both boundary files and gazetteer files. Dividing by ALAND / 1e6 gives people per square kilometre with no projection involved.

Polygon area is a worse substitute than it looks. Across 84,532 populated tracts, the equal-area polygon area was close to land area at the median โ€” a ratio of 1.004 โ€” but 2,832 tracts had polygons at least 25% larger than their land, and 472 more than twice as large. Those are the coastal and lakeside tracts, where the polygon includes water.

6. If you must use geometry, project to an equal-area CRS

For custom polygons with no ALAND, reproject first:

area_km2 = gdf.to_crs("EPSG:6933").area / 1e6       # or EPSG:5070 for the conterminous US

Mercator-style and longitudeโ€“latitude coordinates are not equal area, and any density built on them is distorted by latitude.

7. Never divide by square degrees

The area of a square degree shrinks towards the poles. Measured from the tract file, the median tract in Florida covers 10,912 kmยฒ per square degree of its polygon, in Minnesota 8,763 and in Alaska 5,991. A density in people per square degree therefore favours northern areas.

The damage is easy to underestimate. Ranking tracts by people per square degree gave a rank correlation of 0.997 with the true density, and still moved 2,079 tracts by more than five percentile points โ€” enough to change their colour class on a map.

8. Carry the uncertainty with the rate

A rate from survey counts has a margin of error of its own, and the median tract poverty rate has one of 5.6 percentage points. Compute it alongside the rate rather than as an afterthought; the margins guide has the formula and its fallback.

Table comparing tract poverty rates computed with the poverty universe and with total population: tracts whose rate moves, top-100 overlap and tracts whose universe is under half their population.
The two denominators agree for most tracts and disagree exactly where rankings and extremes are decided.

Code examples

Example 1 โ€” a rate that cannot produce 0 or infinity by accident

SENTINELS = [-222222222, -333333333, -555555555, -666666666, -888888888, -999999999]


def add_rate(frame, numerator, denominator, name):
    """numerator / denominator, with sentinels and zero denominators as missing, never 0 or inf."""
    num = frame[numerator].where(~frame[numerator].isin(SENTINELS))
    den = frame[denominator].where(~frame[denominator].isin(SENTINELS))
    rate = num / den.where(den > 0)
    out = frame.assign(**{name: rate})
    print(f"{name}: {rate.notna().sum():,} defined, {int((den == 0).sum()):,} with a zero denominator, "
          f"{int(num.isna().sum() + den.isna().sum()):,} sentinels")
    return out
poverty_rate: 84,341 defined, 1,040 with a zero denominator, 0 sentinels

The printed line is the audit trail. A rate column with more undefined values than zero denominators plus sentinels has lost rows somewhere else.

Example 2 โ€” test the denominator choice

def compare_denominators(frame, numerator, universe, population, threshold=0.05):
    """How much a rate changes when total population replaces the proper universe."""
    ok = frame[frame[universe] > 0]
    by_universe = ok[numerator] / ok[universe]
    by_population = ok[numerator] / ok[population]
    moved = ((by_universe - by_population).abs() > threshold).sum()
    coverage = ok[universe] / ok[population]
    top_u = set(by_universe.nlargest(100).index)
    top_p = set(by_population.nlargest(100).index)
    print(f"universe / population: median {coverage.median():.3f}; below 0.5 in {int((coverage < 0.5).sum()):,} areas")
    print(f"rates moving more than {threshold:.0%} points: {moved:,}; top-100 lists share {len(top_u & top_p)}")
    return coverage
universe / population: median 0.997; below 0.5 in 623 areas
rates moving more than 5% points: 1,378; top-100 lists share 52

The returned coverage ratio is worth keeping as a column. Areas where the universe is a small share of the population are the ones to annotate on the map.

Example 3 โ€” check a density before mapping it

import warnings

from scipy.stats import spearmanr


def density_check(gdf, population, land="ALAND", equal_area="EPSG:6933"):
    """Density from land area against density from polygon area, in degrees and in an equal-area CRS."""
    g = gdf[(gdf[land] > 0) & (gdf[population] > 0)]
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        degrees = g.geometry.area
    polygon_m2 = g.to_crs(equal_area).area
    ratio = polygon_m2 / g[land]
    true_density = g[population] / (g[land] / 1e6)
    degree_density = g[population] / degrees
    moved = ((true_density.rank(pct=True) - degree_density.rank(pct=True)).abs() > 0.05).sum()
    print(f"polygon / land area: median {ratio.median():.3f}, over 1.25 in {int((ratio > 1.25).sum()):,}, "
          f"over 2 in {int((ratio > 2).sum()):,}")
    print(f"density from square degrees: rank correlation {spearmanr(degree_density, true_density)[0]:.3f}, "
          f"{moved:,} of {len(g):,} areas move more than 5 percentile points")
    return ratio
polygon / land area: median 1.004, over 1.25 in 2,832, over 2 in 472
density from square degrees: rank correlation 0.997, 2,079 of 84,532 areas move more than 5 percentile points

The warning suppressed inside is GeoPandas saying exactly what goes wrong โ€” area in a geographic CRS is likely incorrect. The function computes it on purpose, to measure how wrong.

Explanation

Why sentinels survive arithmetic

The ACS publishes every value in a numeric column, so "cannot be computed" has to be a number too. Codes such as โˆ’666666666 were chosen to be impossible as real estimates, not to be safe in calculations. A pandas mean, sum or division treats them as ordinary numbers, and because they are enormous, a single one outweighs thousands of real values: 1,547 of them turned a mean income of 83,552 into โˆ’11,997,152.

Why the universe matters more in some places than others

For most tracts, the population whose poverty status is determined is nearly everyone, and the two denominators give almost the same rate. The difference concentrates where institutions and group quarters are โ€” prisons, universities, military bases โ€” which are also the tracts with the most unusual rates. That is why a small median difference coexists with a top-100 list that changes by half.

Why polygon area is not land area

Census polygons follow legal and statistical boundaries, which often run through lakes, bays and rivers. The cartographic boundary files clip to the shoreline at a generalised scale, but inland water and detailed coastlines remain. ALAND is computed from the full-resolution geography with water removed, which makes it the right denominator whatever file you draw with.

Why a 0.997 correlation still changes the map

A rank correlation summarises the whole ordering, and most of it survives any monotonic-ish distortion. A choropleth depends on local order near class breaks. Moving 2,079 tracts by five or more percentile points is enough to move many of them into another quintile or decile, concentrated at the northern and southern extremes of the country where the square-degree error is largest.

Bar chart of the median area in square kilometres of one square degree of tract polygon in Florida, Texas, New York, Minnesota and Alaska.
A density per square degree is a density per a different amount of ground in every state.

Edge cases or notes

  • ALAND is in square metres. Divide by 1,000,000 for kmยฒ, or by 2,589,988 for square miles.
  • Water-only tracts have ALAND of zero. Mask them before dividing, exactly like zero populations.
  • Medians and averages are already normalised. Do not divide median household income by population.
  • Percentages of a percentage are not rates. Check whether a published column is already a share before dividing again.
  • Density at block level uses ALAND20 in 2020 block files; the suffix marks the vintage.
  • Rates in small areas are unstable even with the right denominator; see the small numbers guide before ranking them.
  • UK data carries the same choices. Nomis tables state their population base โ€” all usual residents, households, people aged 16 and over โ€” which is the denominator to use.

FAQ

How do I calculate a rate from census counts in pandas?

Divide the count by its table's universe, after replacing sentinel values with NaN and masking zero denominators: numerator / denominator.where(denominator > 0). That left 1,040 tract poverty rates undefined rather than zero or infinite.

Should I divide by total population?

Only for totals that apply to everyone. For characteristics measured on part of the population, such as poverty status, total population moved 1,378 tract rates by more than five percentage points.

How do I calculate population density for census tracts?

Divide population by ALAND divided by 1,000,000 to get people per square kilometre. ALAND excludes water and needs no projection.

Can I use the polygon area instead of ALAND?

Only after projecting to an equal-area CRS, and even then water inflates coastal tracts: 2,832 populated tracts had polygons at least 25% larger than their land area.

Why is my average income negative?

The column contains sentinel codes such as โˆ’666666666. Across tracts, 1,547 of them turned a mean of 83,552 dollars into โˆ’11,997,152.

What should a zero denominator produce?

A missing value. Zero over zero is undefined, and a count over zero is infinite; neither belongs on a map or in a ranking.