How to Calculate Population-Weighted Centroids

Problem statement

A polygon's centroid is the centre of its land, not of its people. For a compact urban tract the two are close. For a large tract with a town at one end, they can be kilometres apart. Distances, travel times and nearest-facility assignments computed from centroids then start from the wrong place.

Measured on Delaware's 258 populated census tracts, the geometric centroid was a median 429 m from the population-weighted centre, and more than 1 km away for 69 tracts. The largest gap was 9.6 km, in a 355 km² tract of 3,722 people. For the 3,100 counties of the contiguous US the median gap was 5.5 km, and in Nye County, Nevada, it was 182.5 km.

The Census Bureau publishes population centres for tracts, block groups and counties, computed from block counts. The same numbers can be reproduced from the TIGER/Line block file in a few lines, which means you can compute them for any zones you like.

Quick answer

import numpy as np
import pandas as pd
import geopandas as gpd

blocks = gpd.read_file("tl_2020_10_tabblock20.zip", ignore_geometry=True,
                       columns=["GEOID20", "POP20", "INTPTLAT20", "INTPTLON20"])
lat = blocks["INTPTLAT20"].astype(float)
lon = blocks["INTPTLON20"].astype(float)
w = blocks["POP20"]
wc = w * np.cos(np.radians(lat))
tract = blocks["GEOID20"].str[:11]

centres = pd.DataFrame({
    "population": w.groupby(tract).sum(),
    "lat": (w * lat).groupby(tract).sum() / w.groupby(tract).sum(),
    "lon": (wc * lon).groupby(tract).sum() / wc.groupby(tract).sum(),
}).query("population > 0")

Measured: every one of the 258 centres agreed with the Census Bureau's published CenPop2020_Mean_TR10.txt to within 0.14 m. The calculation took about 6 ms.

count     258
median    429
75%     1,020
max     9,637

Step-by-step solution

1. Use the smallest populated units that have a point

The 2020 TIGER/Line tabulation block files carry both the population (POP20) and an internal point (INTPTLAT20, INTPTLON20) for every block. Delaware has 20,198 of them. Blocks are small enough that the internal point is a good stand-in for where the block's residents are.

Block groups would work, but each is one point for 600 to 3,000 people. The coarser the input, the more the result drifts back towards the geometric centroid.

2. Apply the Census Bureau's formula

The published centres of population use a weighted mean latitude and a longitude weighted by population × cos(latitude):

lat = Σ(w · lat) / Σ w
lon = Σ(w · cos(lat) · lon) / Σ(w · cos(lat))

The cosine term accounts for degrees of longitude getting shorter away from the equator. Within Delaware it barely matters: dropping it moved tract centres by at most 2.6 m. Across a large country it matters more.

Reproducing the published file checks both the formula and your data. The tract centres matched to 0.14 m, and Delaware's three county centres to 0.02, 0.05 and 0.08 m.

3. Or project the points and take a weighted mean

pts = gpd.GeoDataFrame(blocks, geometry=gpd.points_from_xy(lon, lat), crs=4269).to_crs(5070)

A population-weighted mean of projected x and y in an equal-area projection (EPSG:5070 for the US) came within 1.9 m of the published centres. Use it when the result has to be a geometry anyway, or when the zones are not census units.

4. Drop units with no population

Four Delaware tracts have no residents. Their weights sum to zero and the formula divides by zero. The published file still lists coordinates for them, with a population of 0. Filter with population > 0 and decide explicitly what a zero-population zone should use, usually its geometric centroid or nothing.

5. Measure how far the centroid was from the people

shift = centroid_shift(tracts, projected_mean_centre(pts, "tract", "POP20"))
shift["shift_m"].describe()

count     258
median    429
75%     1,020
max     9,637

The shift grows with size: across tracts, the rank correlation between shift and land area was 0.77. Relative to tract size it is more stable, with a median shift of 0.17 × √(land area).

6. Check whether each centre falls inside its polygon

A weighted mean is not guaranteed to land inside its zone. Two of Delaware's 258 tract centres fell outside their own tract, 62 m and 109 m beyond its edge. Of the 3,100 contiguous-US counties, 13 had a population centre outside the county. One is Baltimore County, Maryland, which wraps around the separate city of Baltimore. Twelve counties' geometric centroids also fell outside the county.

If the centre must be a location in the zone, snap it to the nearest populated block inside the zone rather than using the mean.

7. Use the centres where distances start

Replace centroids with population centres as the origins for distance, travel-time and access calculations. For small urban tracts the change is a few hundred metres; for large rural ones it can move the nearest facility.

Code examples

The data used on this page

From the Census Bureau: Delaware's 2020 TIGER/Line blocks, which carry POP20 and an internal point, and its tracts; the published 2020 centres of population for Delaware's tracts and for every county; and the 2023 cartographic county boundaries at 1:500,000.

import os

import requests

G = "https://www2.census.gov/geo"
FILES = {
    "tl_2020_10_tabblock20.zip": f"{G}/tiger/TIGER2020/TABBLOCK20/tl_2020_10_tabblock20.zip",
    "tl_2020_10_tract.zip": f"{G}/tiger/TIGER2020/TRACT/tl_2020_10_tract.zip",
    "CenPop2020_Mean_TR10.txt": f"{G}/docs/reference/cenpop2020/tract/CenPop2020_Mean_TR10.txt",
    "CenPop2020_Mean_CO.txt": f"{G}/docs/reference/cenpop2020/county/CenPop2020_Mean_CO.txt",
    "cb_2023_us_county_500k.zip": f"{G}/tiger/GENZ2023/shp/cb_2023_us_county_500k.zip",
}
for name, url in FILES.items():
    if not os.path.exists(name):
        r = requests.get(url, timeout=600)
        r.raise_for_status()
        open(name, "wb").write(r.content)
print({name: f"{os.path.getsize(name) / 1e6:.2f} MB" for name in FILES})

{'tl_2020_10_tabblock20.zip': '21.47 MB', 'tl_2020_10_tract.zip': '1.08 MB', 'CenPop2020_Mean_TR10.txt': '0.01 MB', 'CenPop2020_Mean_CO.txt': '0.17 MB', 'cb_2023_us_county_500k.zip': '11.63 MB'}

Example 1 — the Census Bureau's mean centre for any grouping

import numpy as np
import pandas as pd


def census_mean_centre(points, group, weight, lat="lat", lon="lon"):
    """Population-weighted mean centre using the Census Bureau's formula."""
    phi = np.radians(points[lat])
    w = points[weight]
    wc = w * np.cos(phi)
    sums = pd.DataFrame({
        "w": w, "w_lat": w * points[lat], "wc": wc, "wc_lon": wc * points[lon],
        group: points[group],
    }).groupby(group).sum()
    out = pd.DataFrame({
        "population": sums["w"],
        "lat": sums["w_lat"] / sums["w"],
        "lon": sums["wc_lon"] / sums["wc"],
    })
    return out[out["population"] > 0]


print(census_mean_centre(blocks.assign(lat=lat, lon=lon, tract=tract), "tract", "POP20").head(2))

             population        lat        lon
tract
10001040100        7315  39.237913 -75.680801
10001040201        5446  39.295250 -75.626417

The published file lists 39.237913, −75.680800 for the first tract. The difference in the sixth decimal place of longitude is 0.09 m.

Example 2 — a projected version that returns geometries

import geopandas as gpd


def projected_mean_centre(points, group, weight, crs=5070):
    """Population-weighted mean of projected coordinates, as a GeoDataFrame."""
    p = points.to_crs(crs)
    frame = pd.DataFrame({"x": p.geometry.x * p[weight], "y": p.geometry.y * p[weight],
                          "w": p[weight], group: p[group]}).groupby(group).sum()
    frame = frame[frame["w"] > 0]
    geometry = gpd.points_from_xy(frame["x"] / frame["w"], frame["y"] / frame["w"])
    return gpd.GeoDataFrame({"population": frame["w"]}, geometry=geometry, crs=crs, index=frame.index)

The grouping column can be anything present on the points: a tract prefix, a school district from a spatial join, or an H3 cell. That is the advantage over the published files, which exist only for census units.

Example 3 — the shift from centroid to population centre

from pyproj import Geod


def centroid_shift(polygons, centres, crs=5070):
    """Geodesic distance in metres from each polygon's centroid to its population centre."""
    geod = Geod(ellps="WGS84")
    centroids = polygons.to_crs(crs).centroid.to_crs(4326)
    centres = centres.to_crs(4326).reindex(centroids.index)
    _, _, dist = geod.inv(centroids.x.to_numpy(), centroids.y.to_numpy(),
                          centres.geometry.x.to_numpy(), centres.geometry.y.to_numpy())
    inside = polygons.to_crs(4326).geometry.contains(centres.geometry)
    return pd.DataFrame({"shift_m": dist, "centre_inside": inside}, index=centroids.index)

The centroid is computed in a projected CRS and converted afterwards. For contiguous-US counties, a centroid computed directly on longitude and latitude was a median 12 m and at most 1,203 m from the projected one. The difference is small but avoidable.

Example 4 — checking every Delaware figure on this page

from scipy.stats import spearmanr

geod = Geod(ellps="WGS84")
pub = pd.read_csv("CenPop2020_Mean_TR10.txt", dtype=str, encoding="utf-8-sig")
pub.index = pub.STATEFP + pub.COUNTYFP + pub.TRACTCE
pub = pub[pub.POPULATION.astype(int) > 0]
pts = gpd.GeoDataFrame(blocks.assign(lat=lat, lon=lon, tract=tract, county=blocks.GEOID20.str[:5]),
                       geometry=gpd.points_from_xy(lon, lat), crs=4269)


def gap_m(a_lat, a_lon, b):
    return geod.inv(a_lon.to_numpy(), a_lat.to_numpy(),
                    b.LONGITUDE.astype(float).reindex(a_lat.index).to_numpy(),
                    b.LATITUDE.astype(float).reindex(a_lat.index).to_numpy())[2]


ours = census_mean_centre(pts, "tract", "POP20")
print(f"{len(ours)} populated tracts; largest gap to the published centres {gap_m(ours.lat, ours.lon, pub).max():.2f} m")
flat = (w * lon).groupby(tract).sum() / w.groupby(tract).sum()
print(f"without the cosine term: tract centres move at most "
      f"{geod.inv(flat[ours.index].to_numpy(), ours.lat.to_numpy(), ours.lon.to_numpy(), ours.lat.to_numpy())[2].max():.1f} m")
co_pub = pd.read_csv("CenPop2020_Mean_CO.txt", dtype=str, encoding="utf-8-sig")
co_pub.index = co_pub.STATEFP + co_pub.COUNTYFP
co = census_mean_centre(pts, "county", "POP20")
print("Delaware counties:", ", ".join(f"{x:.2f} m" for x in gap_m(co.lat, co.lon, co_pub)))
proj = projected_mean_centre(pts, "tract", "POP20").to_crs(4326)
print(f"projected mean: at most {gap_m(proj.geometry.y, proj.geometry.x, pub).max():.1f} m from the published centres")
print("tracts with no residents:", int((w.groupby(tract).sum() == 0).sum()))

tracts = gpd.read_file("tl_2020_10_tract.zip").set_index("GEOID")
shift = centroid_shift(tracts.loc[proj.index], projected_mean_centre(pts, "tract", "POP20"))
print(shift["shift_m"].describe(percentiles=[0.25, 0.5, 0.75]).round(0).to_string())
print("more than 1 km:", int((shift.shift_m > 1000).sum()))
area = tracts.loc[shift.index, "ALAND"].astype(float)                   # land area, m²
print(f"rank correlation with land area {spearmanr(shift.shift_m, area)[0]:.2f}; "
      f"median shift / sqrt(area) {(shift.shift_m / np.sqrt(area)).median():.2f}")
big = shift.shift_m.idxmax()
print(f"largest: {big}, {tracts.to_crs(5070).area[big] / 1e6:.0f} km², {int(ours.population[big]):,} people")
centres = projected_mean_centre(pts, "tract", "POP20")
outside = tracts.loc[shift.index[~shift.centre_inside]].to_crs(5070)
print("centres outside their tract:", {i: f"{outside.geometry[i].distance(centres.geometry[i]):.0f} m beyond the edge"
                                         for i in outside.index})

258 populated tracts; largest gap to the published centres 0.14 m
without the cosine term: tract centres move at most 2.6 m
Delaware counties: 0.02 m, 0.05 m, 0.08 m
projected mean: at most 1.9 m from the published centres
tracts with no residents: 4
count     258.0
mean      770.0
std       965.0
min         4.0
25%       172.0
50%       429.0
75%      1020.0
max      9637.0
more than 1 km: 69
rank correlation with land area 0.77; median shift / sqrt(area) 0.17
largest: 10001043202, 355 km², 3,722 people
centres outside their tract: {'10003016309': '62 m beyond the edge', '10005051313': '109 m beyond the edge'}

Example 5 — the same comparison for every contiguous-US county

The published county centres are national, so the same check runs on the Census Bureau's 2023 cartographic county boundaries. Connecticut's current planning regions have no 2020 centre and drop out.

counties = gpd.read_file("cb_2023_us_county_500k.zip")
counties = counties[~counties.STATEFP.isin(["02", "15", "60", "66", "69", "72", "78"])].set_index("GEOID")
co_pub = pd.read_csv("CenPop2020_Mean_CO.txt", dtype=str, encoding="utf-8-sig")
co_pub.index = co_pub.STATEFP + co_pub.COUNTYFP
unmatched = counties.index.difference(co_pub.index)
print(f"{len(counties):,} contiguous-US counties; {len(unmatched)} with no published centre:",
      sorted(counties.loc[unmatched, "STATEFP"].unique()))
counties = counties.drop(unmatched)
centres = gpd.GeoDataFrame(geometry=gpd.points_from_xy(co_pub.LONGITUDE.astype(float), co_pub.LATITUDE.astype(float)),
                           index=co_pub.index, crs=4269).loc[counties.index]
shift = centroid_shift(counties, centres)
names = counties.NAME + ", " + counties.STUSPS
print(f"{len(shift):,} counties: median gap {shift.shift_m.median() / 1000:.1f} km; "
      f"more than 10 km: {(shift.shift_m > 10_000).sum()}")
print((shift.shift_m.nlargest(3) / 1000).round(1).rename(names).to_string())
out = shift.index[~shift.centre_inside]
print(f"centre outside its county: {len(out)}, including", ", ".join(names[out][:5]))
cent = counties.to_crs(5070).centroid
print("geometric centroid outside its county:", int((~counties.to_crs(5070).contains(cent)).sum()))
deg = counties.centroid                                       # geographic CRS: warns, and is off
dd = Geod(ellps="WGS84").inv(deg.x.to_numpy(), deg.y.to_numpy(),
                             cent.to_crs(4269).x.to_numpy(), cent.to_crs(4269).y.to_numpy())[2]
print(f"centroid in degrees vs projected: median {np.median(dd):.0f} m, max {dd.max():,.0f} m")

3,109 contiguous-US counties; 9 with no published centre: ['09']
3,100 counties: median gap 5.5 km; more than 10 km: 808
GEOID
Nye, NV               182.5
San Bernardino, CA    130.5
Washoe, NV            125.1
centre outside its county: 13, including Baltimore, MD, Island, WA, Wise, VA, Roanoke, VA, Lee, IA
geometric centroid outside its county: 12
centroid in degrees vs projected: median 12 m, max 1,203 m

Explanation

Why a centroid describes land, not people

The geometric centroid weights every square metre equally, and census zones are drawn to hold similar numbers of people, not similar areas. A rural tract has to be large to reach its population, and its people gather in villages and along roads. The centroid lands in whatever the tract has most of, which is usually fields.

Why the gap scales with the zone

In a large zone the people have more room to be off-centre, which is why the tract shift correlated with land area at 0.77. At county scale, the largest gaps are desert counties with one city at the edge. Nye County, Nevada, is 182.5 km; San Bernardino County, California, 130.5 km; Washoe County, Nevada, 125.1 km. 808 of the 3,100 counties had a gap of more than 10 km.

Why the Census formula uses cos(latitude)

A degree of longitude spans about 111 km at the equator and less further north: about 86 km at Delaware's latitude. Averaging raw longitudes gives a northern resident's degree as much weight as a southern resident's. The cosine weight corrects that without projecting. The published file uses it, which is why the formula with the term reproduced the file to 0.14 m.

Why a mean centre can fall outside

A mean is a balance point, not a location in the set. A county shaped like a ring, or one whose people live at two far ends, balances in a place where nobody lives. Baltimore County's residents surround the city of Baltimore, so their mean lands inside the city, which is a separate county-equivalent.

Why blocks and not block groups

A block group's internal point stands in for up to 3,000 people. Using it reintroduces the problem this calculation is meant to remove, only at a smaller scale. The block file is the finest free population surface the Census Bureau publishes, and for Delaware it is 21 MB zipped.

Edge cases or notes

• Convert the internal point strings. INTPTLAT20 is text like +39.2379130; astype(float) handles the sign.
• Zero-population zones divide by zero. Four Delaware tracts have none; filter them before dividing.
• The published file name is per state. CenPop2020_Mean_TR10.txt is Delaware's tracts, not a 2010 file.
• Connecticut counties changed in 2022. The 2020 centres of population use the old counties, so 9 current planning regions had no match.
• A centre can be outside its zone. 13 contiguous-US counties and 2 Delaware tracts had one.
• Project before computing a centroid. GeoPandas warns on geographic CRSs, and the county centroids moved by up to 1.2 km.
• Use block counts from the same census as the zones. 2010 blocks will not sum to 2020 tracts.
• A median centre is a different statistic. It minimises total distance and does not have to agree with the mean.

Internal links

• How to get polygon centroids in GeoPandas — the geometric version this guide replaces
• How to redistribute population with dasymetric mapping — where people are, when blocks are not enough
• How to count the population within reach of each facility — a calculation that should start from these centres
• How to measure distance to the nearest facility for every home — the distances that move with the origin
• How to measure distance accurately in Python — geodesic versus projected distance
• Census geographies explained: blocks, tracts, output areas and why they nest — why blocks sum to tracts
• How to download census boundaries in Python — getting the block file
• Choosing a map projection for display — which projections preserve area
• Urban extent explained — where a city's population ends

FAQ

What is a population-weighted centroid?

The average location of a zone's residents, computed from small populated units such as census blocks. It is where the people are centred, as opposed to the geometric centroid, which is the centre of the land.

Does the Census Bureau publish them?

Yes. The 2020 centres of population are published for tracts, block groups and counties under www2.census.gov/geo/docs/reference/cenpop2020/. Computing them from blocks reproduced the tract file to within 0.14 m.

How far is the centroid from the population centre?

For Delaware tracts, a median of 429 m and a maximum of 9.6 km. For contiguous-US counties, a median of 5.5 km and a maximum of 182.5 km, in Nye County, Nevada.

Should I project the coordinates first?

Either way works if done correctly. A weighted mean of projected coordinates in EPSG:5070 was within 1.9 m of the published centres, and the Census formula with a cosine weight on longitude was within 0.14 m.

Can the population centre fall outside the polygon?

Yes. It happened for 2 of 258 Delaware tracts and 13 of 3,100 counties, including Baltimore County, whose residents surround a city that is a separate county-equivalent.

Can I use block groups instead of blocks?

You can, but each block group's point represents up to 3,000 people, so the result is pulled back towards the geometric centroid. Use blocks when they are available.