11 min read 2402 words

How to align EPSG:4326 and EPSG:3857 for solar site mapping

Scenario / symptom: gpd.overlay(parcels, raster_bbox, how="intersection") returns an empty or near-zero-area GeoDataFrame (intersection.area.sum() prints 0.0 or a nonsensical value), even though the parcel and the orthomosaic clearly cover the same field. This failure lands in the CRS alignment stage — the moment EPSG:4326 (WGS84) parcel boundaries are overlaid against an EPSG:3857 (Web Mercator) satellite tile without an explicit transformation between them. It is a special case of the CRS drift problem covered by the parent workflow: degrees are compared against metres, the geometries never intersect, and the pipeline produces a confident-but-wrong answer instead of raising an exception.

Solar development pipelines routinely ingest heterogeneous spatial assets: parcel boundaries, interconnection queue data, and environmental constraint layers arrive in EPSG:4326, while satellite orthomosaics, web-mapped terrain models, and utility distribution overlays default to EPSG:3857. When these layers are overlaid without explicit transformation, site boundaries shift by hundreds of metres, irradiance footprints misalign with regulatory setbacks, and capacity-factor models return invalid geometries. The fix is to make the coordinate frame explicit, transform once into a single working CRS, and gate every output with a latitude-aware area check.

Root-cause analysis

The misalignment is not a single bug — it is three compounding causes that each pass silently on their own:

  1. Implicit CRS assumption. geopandas.read_file() and rasterio.open() silently inherit or guess a coordinate frame when metadata is absent. A shapefile written without a .prj sidecar loads with parcels.crs is None, so no transformation is attempted and the raw degree values flow straight into the overlay.
  2. Planar vs. angular arithmetic. Overlay operations (overlay, clip, intersection) treat coordinates as planar numbers. EPSG:4326 stores position in decimal degrees (≈ ±180 / ±90), while EPSG:3857 stores metres (millions). The two number ranges never coincide, so the spatial predicate returns no intersection.
  3. Web Mercator scale drift. EPSG:3857 preserves shape at the equator but inflates area with latitude. The areal scale factor for a conformal Mercator frame is

so at 45°N a footprint computed directly in EPSG:3857 is inflated by roughly 2× in area. A solar array measured in that frame yields inaccurate MW estimates and can violate interconnection filing tolerances even after the layers visually line up.

Establishing a consistent baseline means tagging every dataset explicitly on ingestion, transforming into one metric working CRS, then validating area against a known-good projection before any geometry feeds a siting or regulatory boundary overlay.

How three silent CRS causes converge on a misaligned overlay, and the fix that resolves them Left: three warning boxes — implicit CRS assumption, planar versus angular arithmetic, and Web Mercator scale drift — each arrow into a central warning node, the misaligned overlay that returns an empty intersection. The overlay arrows into a success node on the right: a single working CRS reached via an explicit CRS tag, one pyproj.Transformer call, and a make_valid repair. Three silent causes converge on one misaligned overlay — and the contract that fixes it Implicit CRS assumption parcels.crs is None Planar vs angular arithmetic degrees compared to metres Web Mercator scale drift sec²φ area inflation Misaligned overlay empty intersection area.sum() = 0.0 One working CRS 1 · explicit CRS tag 2 · single pyproj.Transformer 3 · make_valid repair Web Mercator areal scale factor sec²φ versus latitude The areal scale factor equals one over the square of the cosine of latitude. It stays near 1.0 through the low latitudes, passes 1.33 at 30 degrees north, reaches exactly 2.0 (a doubling of measured area) at 45 degrees north, and accelerates to 4.0 at 60 degrees north — so a footprint measured directly in EPSG:3857 is inflated more and more with latitude and must be verified in an equal-area or UTM frame. 1× 1× 2× 3× 4× 0° 15° 30° 45° 60°N EPSG:3857 area within 5% — safe to report directly 1.0× at equator 1.33× at 30°N 2.0× at 45°N 4.0× at 60°N Site latitude φ Areal scale factor sec²φ

Pre-flight validation

Surface the root cause before the overlay runs. The check below inspects both inputs, flags an undefined parcel CRS, detects the degrees-vs-metres mismatch by comparing coordinate magnitudes, and refuses to proceed until both layers share a frame.

python
import geopandas as gpd
import rasterio
import pyproj


def preflight_crs_check(parcel_path: str, raster_path: str) -> dict:
    """Surface CRS mismatches before any overlay executes.

    Returns a diagnostic dict; raises ValueError on an unrecoverable mismatch.
    """
    parcels = gpd.read_file(parcel_path)
    with rasterio.open(raster_path) as src:
        raster_crs = src.crs
        left, bottom, right, top = src.bounds

    report = {
        "parcel_crs": str(parcels.crs),
        "raster_crs": str(raster_crs),
        "parcel_bounds": tuple(round(v, 3) for v in parcels.total_bounds),
        "raster_bounds": (round(left, 1), round(bottom, 1), round(right, 1), round(top, 1)),
    }

    # 1. Undefined parcel CRS — the most common silent failure
    if parcels.crs is None:
        raise ValueError(
            "Parcel CRS is undefined. Assign EPSG:4326 (or the true frame) "
            "explicitly before any geometric operation."
        )

    # 2. Degrees-vs-metres magnitude mismatch (4326 stays within +/-180/+/-90)
    px_max = max(abs(v) for v in parcels.total_bounds)
    rx_max = max(abs(v) for v in (left, bottom, right, top))
    parcel_is_degrees = px_max <= 360
    raster_is_metres = rx_max > 360
    if parcel_is_degrees and raster_is_metres and parcels.crs != raster_crs:
        report["diagnosis"] = (
            "MISMATCH: parcels in angular degrees, raster in projected metres. "
            "Overlay will return empty geometry until both share one CRS."
        )
    else:
        report["diagnosis"] = "OK: coordinate magnitudes are compatible."

    report["needs_transform"] = parcels.crs != raster_crs
    return report

Running preflight_crs_check against a 4326 parcel and a 3857 ortho prints the MISMATCH diagnosis instead of letting the overlay quietly return 0.0.

Fix implementation

Align the layers by declaring the parcel CRS explicitly, transforming geometry with a single pyproj.Transformer (vectorised across every vertex via shapely.ops.transform), repairing topology, then performing the intersection in one shared working frame. EPSG:3857 is acceptable as the display working frame, but area is verified separately against an equal-area or UTM frame in the next step.

python
import geopandas as gpd
import rasterio
import pyproj
from rasterio.windows import Window
from shapely.geometry import box
from shapely.ops import transform as shapely_transform
from shapely.validation import make_valid
import logging

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")


def align_spatial_assets(parcel_path: str, raster_path: str, target_crs: str = "EPSG:3857"):
    # 1. Explicit CRS assignment on ingestion
    parcels = gpd.read_file(parcel_path)
    if parcels.crs is None:
        parcels.set_crs("EPSG:4326", inplace=True)
        logging.warning("Parcel CRS was undefined. Defaulted to EPSG:4326.")

    # 2. Raster metadata extraction with windowed bounds for memory efficiency
    with rasterio.open(raster_path) as src:
        src_crs = src.crs
        raster_window = Window(0, 0, src.width, src.height)
        raster_bounds = src.window_bounds(raster_window)

    # 3. Fast, memory-safe coordinate transformation using pyproj.Transformer.
    # shapely.ops.transform applies the transformer to every coordinate of an
    # arbitrary geometry (Point, Polygon, MultiPolygon, ...), so parcel
    # boundaries are preserved instead of collapsed to a single (x, y) tuple.
    transformer = pyproj.Transformer.from_crs(
        parcels.crs, target_crs, always_xy=True
    )
    parcels_aligned = parcels.copy()
    parcels_aligned.geometry = parcels.geometry.apply(
        lambda geom: shapely_transform(transformer.transform, geom) if geom else geom
    )
    parcels_aligned.set_crs(target_crs, inplace=True)

    # 4. Raster bounding box alignment — build a Shapely box from the bounds
    # captured inside the `with` block so the source CRS metadata is retained.
    raster_bbox = gpd.GeoDataFrame(
        geometry=[box(*raster_bounds)],
        crs=src_crs
    ).to_crs(target_crs)

    # 5. Spatial intersection with topology repair
    parcels_aligned.geometry = parcels_aligned.geometry.apply(make_valid)
    intersection = gpd.overlay(parcels_aligned, raster_bbox, how="intersection")

    logging.info(f"Alignment complete. Valid geometries: {intersection.geometry.is_valid.sum()}")
    return intersection

always_xy=True pins (lon, lat) ordering and eliminates the legacy axis-order bug in pyproj 6+. Using a single pyproj.Transformer for batch conversion is faster and far less memory-intensive than repeated per-row GeoDataFrame.to_crs() calls, and make_valid repairs the self-intersections that surface when degree-precision rings are reprojected to metres.

Fallback routing & performance tuning

For national-scale screening, CI/CD runs, and memory-constrained cloud nodes, layer these strategies on top of the core fix:

  • Equal-area fallback for any MW figure. EPSG:3857 is for display only. Before reporting array area, reproject to parcels.estimate_utm_crs() or an Albers equal-area frame (e.g. EPSG:5070 for CONUS). If EPSG:3857 distortion exceeds 5% at the site latitude, route the area calculation through the equal-area frame automatically and log the trigger.
  • Windowed raster I/O. Never load a multi-gigabyte orthomosaic in full. Use rasterio.windows.from_bounds() to read only the tile overlapping the parcel extent — this cuts peak memory by 60–80% on large terrain stacks.
  • Spatial index pre-filter. Call parcels_aligned.sindex.query(raster_geom) before overlay() to drop non-overlapping features and avoid O(n²) geometric comparisons during batch site screening, mirroring the index discipline used in proximity and distance calculations.
  • Pin the PROJ stack. Pin pyproj to an exact minor version in requirements.txt (e.g. pyproj==3.7.2) so the bundled PROJ datum database is identical across CI/CD and production, keeping transformations deterministic.
  • Isolate transform failures. Wrap the overlay in a try/except for pyproj.exceptions.ProjError and shapely.errors.TopologicalError; on failure, revert to the source CRS, apply a conservative 50 m buffer to absorb coordinate drift, and queue the asset for manual GIS review rather than dropping it silently.

Downstream validation

Gate every alignment output in CI/CD with an assertion that fails the build on CRS, emptiness, validity, or latitude-distortion regressions — the same audit posture used in grid capacity buffer analysis.

python
from datetime import datetime, timezone
import math


def audit_alignment(intersection_gdf, expected_crs: str = "EPSG:3857",
                    site_latitude_deg: float | None = None,
                    max_mercator_inflation: float = 1.05) -> dict:
    """Assert alignment integrity. Raises AssertionError on any CI/CD-blocking issue."""
    assert intersection_gdf.crs is not None, "Output CRS is undefined."
    assert str(intersection_gdf.crs) == expected_crs, (
        f"CRS drift: expected {expected_crs}, got {intersection_gdf.crs}"
    )
    assert not intersection_gdf.empty, "Empty overlay — inputs did not intersect (CRS mismatch?)."
    invalid = (~intersection_gdf.geometry.is_valid).sum()
    assert invalid == 0, f"{invalid} invalid geometries remain after make_valid."

    audit = {
        "source_crs": "EPSG:4326",
        "target_crs": str(intersection_gdf.crs),
        "transformer_method": "pyproj.Transformer(always_xy=True)",
        "feature_count": len(intersection_gdf),
        "all_valid": bool(invalid == 0),
        "timestamp": datetime.now(timezone.utc).isoformat(),
    }

    # Warn (don't trust 3857 area) when Mercator inflation is material at this latitude
    if site_latitude_deg is not None:
        k_area = 1.0 / (math.cos(math.radians(site_latitude_deg)) ** 2)
        audit["mercator_area_inflation"] = round(k_area, 3)
        assert k_area <= max_mercator_inflation or expected_crs != "EPSG:3857", (
            f"Mercator area inflation {k_area:.2f}x at {site_latitude_deg} deg — "
            "report area from an equal-area/UTM frame, not EPSG:3857."
        )
    return audit

Attaching the returned audit dictionary to each deliverable satisfies ISO 19115 lineage requirements and lets a permitting authority or independent engineer reproduce exactly how a footprint was derived.

Related articles

Back to Coordinate Reference Systems for Energy Projects