Automated State Plane Coordinate System Selection Tool

Selecting the correct State Plane Coordinate System (SPCS) zone is essential for ensuring spatial accuracy and precision in geospatial workflows. However, manually identifying the correct zone among the 124 NAD 83 State Plane zones can be tedious and error-prone, especially when relying on printed maps or reference tables. To address this, I developed an automated Python tool that simplifies this process.

Tool Features:

• Versatile Input: Accepts addresses, city/state combinations, full ZIP codes, or latitude/longitude coordinates.

• Automated Geocoding: Utilizes Nominatim/OpenStreetMap for converting addresses into geographic coordinates.

• Zone Identification: Determines the appropriate SPCS zone using spatial intersection within GeoPandas.

• EPSG Code Retrieval: Fetches exact EPSG codes through a local SQLite database created from the EPSG registry.

• Interactive Visualization: Generates clear, visual maps displaying the precise location within the selected SPCS zone.

Technologies Used:

• Python 3.12

• GeoPandas, Shapely, PyProj, GeoPy

• Matplotlib for map visualization

Data Sources

• US Census Bureau TIGER (2024)

• EPSG registry table that links every county GEOID to its NAD 83 State Plane code

• Geocoder, Nominatim OpenStreetMap

This automated approach significantly reduces the time required for identifying correct SPCS zones from minutes down to seconds and ensures more reliable results by eliminating human error.

StatePlaneCoordinateSystemTool.py

"""

Interactive helper for the NAD 83 State‑Plane Coordinate System (SPCS)

---------------------------------------------------------------------

• Accepts a city/state, street address, or lon/lat pair

• Returns the county GEOID and matching SPCS EPSG code

• Generates a PNG map: red SPCS zone bbox, gray county outlines, black query point

A small SQLite cache speeds up repeat geocoder calls.

Dependencies: GeoPandas · Shapely ≥ 2 · PyProj · Matplotlib · GeoPy

"""

import time

import sqlite3

from pathlib import Path

import re

import geopandas as gpd

import pandas as pd

from shapely.geometry import Point

from geopy.geocoders import Nominatim

import matplotlib.pyplot as plt

from pyproj import CRS

from shapely.geometry import box

import os, sys

import webbrowser

# ---------------------------------------------------------------

# Set-Up

# ---------------------------------------------------------------

Path("cache").mkdir(exist_ok=True)

# single-table cache: stores geocoder results for ~12 h

CACHE = sqlite3.connect("cache/geocode_cache.db")

CUR = CACHE.cursor()

CUR.execute(

    """

    CREATE TABLE IF NOT EXISTS geocode (

        query TEXT PRIMARY KEY,

        lon   REAL,

        lat   REAL,

        county_fips TEXT,

        expires INTEGER

    )

    """

)

CACHE.commit()

COUNTIES_4326 = gpd.read_file("data/tl_2024_us_county.shp").to_crs(4326) # Counties in WGS-84    

SPCS = pd.read_csv("data/state_plane_reference.csv")

geocoder = Nominatim(user_agent="spcs_selector")

# ---------------------------------------------------------------

# Helpers

# ---------------------------------------------------------------

def get_cached(q):

    row = CUR.execute(

        "SELECT lon, lat, county_fips, expires FROM geocode WHERE query = ?",

        (q,),

    ).fetchone()

    if row and row[3] > int(time.time()):

        return row

    return None

def save_cache(q, lon, lat, county):

    CUR.execute(

        "INSERT OR REPLACE INTO geocode VALUES (?,?,?,?,?)",

        (q, lon, lat, county, int(time.time()) + 43200),

    )

    CACHE.commit()

def as_float_pair(text):

    """Return (lon, lat) if the string holds two numbers, else None."""

    parts = text.replace(",", " ").split()

    if len(parts) != 2:

        return None

    try:

        return float(parts[0]), float(parts[1])

    except ValueError:

        return None

def candidate_county(lon, lat):

    pt = Point(lon, lat)

    hit = COUNTIES_4326[COUNTIES_4326.contains(pt)]

    return hit.iloc[0]["GEOID"] if not hit.empty else None

def parse_lonlat(text):

    """Try lon-lat, then lat-lon; return (lon,lat,county) or None."""

    pair = as_float_pair(text)

    if not pair:

        return None

    a, b = pair

    # first interpret as lon, lat

    county = candidate_county(a, b)

    if county:

        return a, b, county

    # then interpret as lat, lon

    county = candidate_county(b, a)

    if county:

        return b, a, county

    return None  # neither order works

def geocode_query(q):

    hit = get_cached(q)

    if hit:

        return hit[0], hit[1], hit[2]

    triple = parse_lonlat(q)

    if triple:

        lon, lat, county = triple

    else:

        loc = geocoder.geocode(q, exactly_one=True, timeout=10)

        if loc is None:

            raise ValueError("Geocoding failed")

        lon, lat = loc.longitude, loc.latitude

        county = candidate_county(lon, lat)

        if county is None:

            raise ValueError("Point is outside US counties")

    save_cache(q, lon, lat, county)

    return lon, lat, county

def epsg_for_county(fips5):

    row = SPCS.loc[SPCS["fips"].astype(str).str.zfill(5) == fips5, "nad83_epsg"]

    return int(row.iloc[0]) if not row.empty else None

# ---------------------------------------------------------------

# Main Loop

# ---------------------------------------------------------------

def visualize_spcs(lon, lat, epsg, county_fips, counties_gdf,

                   query,save_dir="images"):

    """

    Draw the State-Plane zone (bbox), the query point, zoom to the whole

    state, and show a small info panel with CRS details.

    """

    Path(save_dir).mkdir(exist_ok=True)

    # Sanitize the query text into a filesystem-safe slug

    safe = re.sub(r'[^A-Za-z0-9]+', '_', query).strip('_')

# Build map layers ----------------------------------------------------------

    # 1. State extent -------------------------------------------------------------

    state_code = county_fips[:2]

    state_counties = counties_gdf[counties_gdf["GEOID"].str.startswith(state_code)]

    sxmin, symin, sxmax, symax = state_counties.total_bounds

    xbuf = (sxmax - sxmin) * 0.05

    ybuf = (symax - symin) * 0.05

    # 2. EPSG area-of-use bbox (quick visual) --------------------------------------

    w, s, e, n = CRS.from_epsg(epsg).area_of_use.bounds

    zone_poly = box(w, s, e, n)

    zone_gdf  = gpd.GeoDataFrame({"epsg":[epsg]}, geometry=[zone_poly], crs=4326)

    # 3. Query point ---------------------------------------------------------------

    pt_gdf = gpd.GeoDataFrame(geometry=[Point(lon, lat)], crs=4326)

    # 4. title text: include CRS name, EPSG, key parameters ----------------------

    crs   = CRS.from_epsg(epsg)

    meth  = crs.coordinate_operation.method_name

    units = crs.coordinate_system.axis_list[0].unit_name

    params = {p.name: p.value for p in crs.coordinate_operation.params}

    if "Longitude of natural origin" in params:        # Transverse Mercator

        title = (

            f"CRS    : {crs.name}\n"

            f"EPSG   : {epsg}\n"

            f"Method : {meth}\n"

            f"Lon₀   : {params['Longitude of natural origin']:.4f}°\n"

            f"Scale₀ : {params['Scale factor at natural origin']}\n"

            f"Units  : {units}"

        )

    elif "Standard parallel 1" in params:              # Lambert / Albers

        title = (

            f"CRS    : {crs.name}\n"

            f"EPSG   : {epsg}\n"

            f"Method : {meth}\n"

            f"Std ∥  : {params['Standard parallel 1']:.4f}°, "

            f"{params['Standard parallel 2']:.4f}°\n"

            f"Lon₀   : "

            f"{params.get('Longitude of false origin', params.get('Longitude of natural origin')):.4f}°\n"

            f"Units  : {units}"

        )

    else:

        title = (

            f"CRS: {crs.name}\n"

            f"EPSG   : {epsg}\n"

            f"Method: {meth}\n"

            f"Units: {units}"

        )

    # 5. Plot ---------------------------------------------------------------

    fig, ax = plt.subplots(figsize=(6, 6))

    state_counties.boundary.plot(ax=ax, linewidth=0.3, color="#7f7f7f")

    zone_gdf.plot(ax=ax, facecolor="#2A4A70", edgecolor="#2A4A70",

                  alpha=0.25, linewidth=1.2)

    pt_gdf.plot(ax=ax, color="black", markersize=25)

    ax.set_xlim(sxmin - xbuf, sxmax + xbuf)

    ax.set_ylim(symin - ybuf, symax + ybuf)

    ax.set_title(title, pad=12, loc="left", fontdict={"fontsize": 10}    )

    ax.set_axis_off()

    plt.tight_layout()

    out = Path(save_dir) / f"{safe}_EPSG{epsg}.png"

    # save + open PNG

    plt.savefig(out, dpi=300)

    plt.close(fig)

    print(f"✓ map saved → {out}")

    if sys.platform.startswith("win"):

        os.startfile(out)

    else:

        webbrowser.open(out.as_uri())

if __name__ == "__main__":

    # --- Simple Command Line Interface ---

    while True:

        query = input("\nEnter a City and State, address, or lon/lat (Hit enter without input to exit): ").strip()

        if not query:

            break

        try:

            lon, lat, county = geocode_query(query)

            epsg = epsg_for_county(county)

            if epsg:

                print(f"GEOID {county}   EPSG {epsg}")

                visualize_spcs(lon, lat, epsg, county, COUNTIES_4326, query)

            else:

                print("County found but no EPSG code in lookup table")

        except Exception as err:

            print("Error:", err)