Implementing FAR checks with Shapely and GeoPandas

Implementing FAR checks with Shapely and GeoPandas requires calculating the ratio of a building’s total floor area to its underlying parcel area, then comparing that ratio against jurisdiction-specific zoning thresholds. The core workflow loads parcel and building footprint datasets, enforces a metric-projected coordinate reference system (CRS), computes planar areas, aggregates floor areas via spatial joins, and applies vectorized compliance logic. Because FAR is a dimensionless ratio, any deviation from accurate 2D area measurement cascades into false compliance flags. The pipeline below replaces manual desktop GIS workflows with a reproducible, open-source stack that scales across municipal datasets.

Prerequisites & CRS Validation #

FAR validation sits at the intersection of spatial geometry and tabular zoning rules. Before executing any area calculations, verify that all input shapefiles or GeoJSON files use a projected CRS (e.g., UTM, State Plane, or local metric grid). Geographic CRS like EPSG:4326 returns measurements in square degrees, which breaks compliance math and triggers silent calculation errors. Transform datasets using .to_crs() and validate projection status with .crs.is_projected. Consult the GeoPandas projections guide for authoritative transformation workflows.

Additionally, run make_valid on all geometries to repair self-intersections, duplicate vertices, or invalid ring orientations that silently corrupt area outputs. The Shapely validation reference documents how topology errors propagate into spatial operations.

Production-Ready Pipeline #

import geopandas as gpd
import numpy as np
from shapely.geometry import box
from shapely.validation import make_valid

# 1. Load datasets (replace paths with your actual sources)
# parcels = gpd.read_file("zoning_parcels.gpkg")
# buildings = gpd.read_file("building_footprints.gpkg")

# Mock data for immediate reproducibility
parcels = gpd.GeoDataFrame({
    "parcel_id": ["P1", "P2"],
    "max_far": [2.5, 1.8],
    "geometry": [
        box(0, 0, 100, 100),    # 10,000 m² parcel
        box(150, 0, 350, 200),  # 40,000 m² parcel
    ]
}, crs="EPSG:32618")  # UTM Zone 18N (meters)

buildings = gpd.GeoDataFrame({
    "bldg_id": ["B1", "B2", "B3"],
    "stories": [4, 2, 3],
    "geometry": [
        box(10, 10, 60, 60),    # inside parcel P1
        box(160, 10, 260, 110), # inside parcel P2
        box(170, 120, 300, 190) # inside parcel P2
    ]
}, crs="EPSG:32618")

# 2. Enforce valid geometries & verify projected CRS
parcels["geometry"] = parcels.geometry.apply(make_valid)
buildings["geometry"] = buildings.geometry.apply(make_valid)

if not (parcels.crs.is_projected and buildings.crs.is_projected):
    raise ValueError("Both datasets must use a projected CRS for accurate planar area calculations.")

# 3. Calculate footprint areas & spatial join
buildings["footprint_area_m2"] = buildings.geometry.area

# Join buildings to parcels; 'intersects' captures edge overlaps safely
bldg_parcel_joined = gpd.sjoin(
    buildings, parcels, how="inner", predicate="intersects"
)

# 4. Compute total floor area & FAR
# Simple multiplier: footprint × stories (adjust for jurisdictional inclusions/exclusions)
bldg_parcel_joined["total_floor_area_m2"] = (
    bldg_parcel_joined["footprint_area_m2"] * bldg_parcel_joined["stories"]
)

# Parcel area: look up from the parcels GeoDataFrame by parcel_id
# (sjoin uses suffixes _left/_right for non-geometry columns; geometry stays as 'geometry')
parcel_areas = parcels.set_index("parcel_id")["geometry"].area
bldg_parcel_joined["parcel_area_m2"] = bldg_parcel_joined["parcel_id"].map(parcel_areas)
bldg_parcel_joined["calculated_far"] = (
    bldg_parcel_joined["total_floor_area_m2"] / bldg_parcel_joined["parcel_area_m2"]
)

# 5. Apply compliance logic
bldg_parcel_joined["compliant"] = bldg_parcel_joined["calculated_far"] <= bldg_parcel_joined["max_far"]
bldg_parcel_joined["far_variance"] = (
    bldg_parcel_joined["calculated_far"] - bldg_parcel_joined["max_far"]
)

# Output results
print(bldg_parcel_joined[["bldg_id", "parcel_id", "calculated_far", "max_far", "compliant"]])

Step-by-Step Execution Logic #

1. Geometry Sanitization: make_valid resolves topological defects before area computation. Skipping this step often produces negative or inflated areas in legacy municipal datasets.
2. Spatial Join Strategy: gpd.sjoin(predicate="intersects") maps each building footprint to its host parcel(s). After the join, parcel areas are retrieved from the original parcels GeoDataFrame by parcel_id — not from a geometry_right column, which does not exist in the sjoin output. For buildings crossing parcel boundaries, group by parcel_id and distribute floor area proportionally using geometry.intersection() to avoid double-counting.
3. Vectorized Area Calculation: .geometry.area operates on the entire column at once, leveraging NumPy-backed C extensions. Avoid row-wise apply() for area math; it degrades performance on datasets exceeding 10k features.
4. FAR Computation: The ratio divides aggregated floor area by parcel area. Store both raw values and variance metrics (calculated_far - max_far) to support downstream audit trails and exception reporting.

Integrating Zoning Rules & Compliance Flags #

Raw FAR calculations rarely map 1:1 to municipal zoning codes. Jurisdictions frequently exclude basements, mechanical penthouses, or public plazas from floor area totals. When designing automated compliance workflows, the Height & FAR Compliance Logic layer must isolate planar measurements before applying multipliers, exemptions, or overlay district modifiers.

For enterprise deployments, embed these calculations into a broader Rule Engine Design for Zoning & Setback Automation architecture. Decouple spatial geometry operations from policy evaluation: store zoning thresholds in a version-controlled lookup table, apply them via pd.merge(), and output compliance flags as boolean masks. This separation enables planners to update code without redeploying the spatial pipeline, and it simplifies auditability for regulatory reviews.

Getting the Denominator Right #

Most FAR bugs are denominator bugs. The numerator is a sum of floor areas and is at least obviously a sum; the site area looks like parcel.area and is frequently not.

Codes distinguish gross site area — the parcel as recorded — from net site area, which excludes some combination of public rights-of-way inside the parcel boundary, recorded easements, private access strips, and sometimes land under water or above a slope threshold. A parcel with a twenty-foot right-of-way running through it can have a net area ten per cent below its gross, which moves the ratio by the same ten per cent in the direction that matters.

def site_area(parcel, exclusions, basis="net"):
    """Gross is the parcel; net subtracts the excluded areas the code names.

    exclusions: a GeoDataFrame of right-of-way, easement and other excluded
    polygons, already clipped to the parcel and in the working CRS.
    """
    if basis == "gross":
        return parcel.area
    excluded = exclusions.geometry.unary_union
    return parcel.difference(excluded).area if excluded else parcel.area

Two practical cautions. First, clip the exclusions to the parcel before subtracting: a right-of-way polygon extending beyond the parcel will not affect difference, but it will affect any area you compute from the exclusion layer directly for reporting. Second, dissolve overlapping exclusions before subtracting, or an easement crossing a right-of-way will be counted once and reported twice, which produces a report whose components do not sum to its total.

Reporting Components, Not Just the Ratio #

A FAR verdict that reports only a ratio cannot be checked, and a verdict that cannot be checked is one that gets argued about. The components are already computed; persisting them costs nothing.

A useful verdict payload carries the counted floor area with its per-floor breakdown, the excluded areas by category, the site area on the stated basis, the resulting ratio, the threshold, and the citation. With those, a planner can confirm the arithmetic in under a minute and can see immediately whether a disagreement is about measurement or about definition — which are resolved by very different conversations.

verdict = {
    "parcel_id": apn,
    "floor_area_counted_sqft": round(counted, 1),
    "floor_area_excluded": {"basement": 1180.0, "parking": 2400.0},
    "site_area_basis": "net",
    "site_area_sqft": round(net_site, 1),
    "site_area_excluded": {"right_of_way": 2100.0},
    "far": round(counted / net_site, 3),
    "far_limit": 1.25,
    "citation": "§ 17.40.060(A)",
    "outcome": "violation" if counted / net_site > 1.25 else "compliant",
}

Validating the Whole Pipeline on One Parcel #

Before running a county, run one parcel whose answer you know independently — ideally one a planner has assessed by hand — and reconcile every intermediate number, not just the final ratio.

The reconciliation catches the errors that a ratio comparison hides. If the ratio matches but the floor area is ten per cent high and the site area is ten per cent high with it, two errors are cancelling and will stop cancelling on the next parcel. If the site area matches the assessor’s recorded acreage exactly, that is evidence the frame and units are right. And if the counted floor area differs from the planner’s figure, the per-floor breakdown localises the disagreement to a specific floor and usually to a specific definition.

Keep that parcel as a fixture. It costs one test and it is the single most effective guard against the class of change — a new exclusion rule, a library upgrade, a frame migration — that alters every FAR on the site without failing anything else.

Frequently Asked Questions #

Should FAR be computed from a 3D model when one exists? #

Where a model exists and is trustworthy, yes, because it removes the floor-count assumption entirely. But treat the model’s provenance the way you would any other source: a designer’s submitted model is authoritative for a proposal, while a citywide model derived from imagery is a proxy and should carry a proxy’s tolerance.

How are partial floors and mezzanines handled? #

By the code’s definition, which almost always exists and is almost always specific — commonly a minimum head height or a percentage-of-floor-below threshold. Encode it as an inclusion rule with its citation rather than as a rounding convention, so the choice is visible in the verdict.

What if the parcel has no recorded exclusions but the code specifies net area? #

Then net equals gross for that parcel, and the verdict should say so explicitly rather than silently reporting a gross figure under a net label. The distinction matters when the exclusions layer is later corrected, because it identifies which parcels need re-evaluating.

Does the working CRS affect the ratio? #

Not the ratio itself, since both areas scale together, but it certainly affects the reported areas — and those are what get compared against the assessor’s record. Compute in the authoritative projected frame and report in the unit the code uses, as covered in best practices for CRS standardization.

Related #

Part of: Height and FAR compliance logic

• Computing building height from lidar-derived surfaces — the other half of the envelope.
• Computing dwelling unit density per acre in Python — the same denominator question for density.
• Spatial threshold configuration — recording the basis and inclusion list.
• Apportioning split-zoned parcels by area — when one parcel has two ratios.

Performance & Validation Checklist #

• Handle Multi-Parcel Buildings: If a single footprint intersects multiple parcels, use gpd.overlay(how="intersection") to split the building geometry, then allocate floor area by intersection proportion.
• Avoid Floating-Point Drift: Round calculated_far to 3 decimal places before threshold comparison. Use np.isclose() or tolerance bands (<= max_far + 0.01) when municipal codes allow minor rounding variances.
• Index Spatial Data: For datasets >50k features, build spatial indexes with .sindex before joins. GeoPandas automatically leverages R-tree indexes during sjoin, but explicit .sindex calls improve cache locality in iterative workflows.
• Cross-Validate Outputs: Compare pipeline results against a manual QGIS/ArcGIS sample. Verify that parcel_area_m2 matches official assessor records within ±0.5%. Discrepancies usually trace to unprojected inputs or invalid ring orientations.
• Export for Reporting: Write compliant/non-compliant subsets to separate GeoParquet files. Parquet preserves CRS metadata, supports columnar compression, and integrates cleanly with BI dashboards used by planning departments.

By standardizing CRS enforcement, vectorizing spatial joins, and decoupling geometry math from policy rules, teams can scale FAR validation across entire municipalities without proprietary licensing overhead.