SpatiaLite & GeoPackage Automation with Python

Choosing Between EPSG:4326 and a Projected CRS

Store in whatever the data will mostly be used in; measure in a projected system chosen for the extent; interchange in EPSG:4326. Those are three separate…

Store in whatever the data will mostly be used in; measure in a projected system chosen for the extent; interchange in EPSG:4326. Those are three separate decisions, and treating them as one is why so many containers end up computing areas in degrees. The short version: 4326 is the right answer for handing a file to someone else and the wrong answer for anything with units.

This page belongs to the Datum Transformations & Projection Accuracy guide.

Why This Matters

EPSG:4326 is a geographic reference system: coordinates are angles on an ellipsoid, not positions on a plane. Every function that assumes a plane — length, area, buffer, distance — therefore produces a number whose unit is “degrees”, which is not a unit of length and varies with latitude in one axis and not the other.

Nothing raises. ST_Area returns a float, ST_Buffer returns a polygon, ST_Distance returns a number that is even roughly monotonic with real distance. The result is a pipeline that computes plausible nonsense and a report that is wrong by a factor which changes as you move north.

Prerequisites

  • A clear statement of what the container is for — analysis, distribution, or both
  • The geographic extent the data covers
  • pyproj 3.4+ if you want the library to suggest a projected system
  • Familiarity with the registry, from How to Add a Custom CRS to gpkg_spatial_ref_sys

Primary Method

python
# Ask PROJ which projected systems suit an extent, ranked by area of use
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info, query_crs_info


def suggest_projected(min_lon, min_lat, max_lon, max_lat):
    aoi = AreaOfInterest(min_lon, min_lat, max_lon, max_lat)

    utm = query_utm_crs_info(datum_name="WGS 84", area_of_interest=aoi)
    for info in utm:
        print(f"UTM      EPSG:{info.code:<6} {info.name}")

    projected = query_crs_info(
        pj_types="PROJECTED_CRS", area_of_interest=aoi, contains=True
    )
    for info in projected[:5]:
        print(f"national EPSG:{info.code:<6} {info.name}")

For most work the answer is either the national grid where one exists — it is what the surrounding data uses and what local consumers expect — or the UTM zone covering the extent, which is a reasonable default anywhere and stays accurate over a few hundred kilometres.

Three jobs, three answersStorage should use whatever the data is mostly used in, so the common case needs no transformation. Measurement must use a projected system appropriate to the extent, because area and distance have no meaning in degrees. Interchange should use EPSG 4326, because every consumer can read it without a definition. Treating the three as one decision is what produces areas computed in degrees.storagewhat it is mostly used inavoids the common transformindex works without oneoften the national gridmeasurementprojected, suited to the extentmetres, not degreesor a geodesic functionnever 4326interchangeEPSG:4326every consumer reads itneeds no extra definitionthe safe handoverthe three columns need not agree, and collapsing them into one is the usual mistakea container can store one system, measure in another and export in a third
Storing in a projected system and exporting to 4326 at the boundary costs one transformation and removes the whole class of unit errors.

Step-by-Step Walkthrough

1. Demonstrate the problem before choosing

The fastest way to settle an argument about this is to compute the same area both ways.

python
from shapely.geometry import Polygon
from pyproj import Transformer

# A square roughly 1 km on a side, at 51°N
poly_4326 = Polygon([
    (-0.10, 51.50), (-0.10, 51.509), (-0.086, 51.509), (-0.086, 51.50), (-0.10, 51.50)
])
print(f"area in 4326: {poly_4326.area:.10f}")        # square degrees — meaningless

to_27700 = Transformer.from_crs("EPSG:4326", "EPSG:27700", always_xy=True)
from shapely.ops import transform
poly_27700 = transform(to_27700.transform, poly_4326)
print(f"area in 27700: {poly_27700.area:,.1f} m²")   # ~1,000,000

The first number is not “the area in different units” — it is a quantity with no physical interpretation, because a degree of longitude and a degree of latitude are different lengths and the ratio changes as you move.

2. Decide storage from the dominant use

If the container is read mostly by a field application that renders a Web Mercator basemap, storing in EPSG:3857 removes a per-feature transform on every frame. If it is read mostly by an analysis pipeline computing areas, store in the projected system that pipeline measures in. If it is read mostly by external parties, store in 4326.

The decision is about which transformation you are willing to pay for repeatedly, because you will pay for the others once at the boundary.

3. Never measure in a geographic system

Where measurement must happen against 4326-stored data, two options are correct and one is not.

python
# Correct: transform, then measure
projected = transform(to_27700.transform, geom)
area_m2 = projected.area

# Also correct: measure geodesically, no projection involved
from pyproj import Geod
geod = Geod(ellps="WGS84")
area_m2_geodesic, perimeter_m = geod.geometry_area_perimeter(geom)
area_m2_geodesic = abs(area_m2_geodesic)

# Wrong: geom.area on a 4326 geometry

The geodesic route is the more accurate of the two over large extents, because it does not inherit a projection’s distortion. Over a survey-sized area the two agree closely, and the projected route is faster.

Three ways to compute an area from 4326 dataCalling the planar area function directly on geographic coordinates produces square degrees, which has no physical meaning and is wrong by a factor that changes with latitude. Transforming to a projected system first produces metres and is accurate over a survey-sized extent. Using a geodesic area function produces metres without any projection and stays accurate over large extents, at a higher computational cost.planar area on 4326 — wrongsquare degrees · no physical meaning · error changes with latitudetransform, then measure — correct and fastmetres · accurate over a survey-sized extent · one transform per geometrygeodesic area — correct at any sizemetres · no projection distortion · slower, and worth it over large extents
Only the top row is a bug. The lower two are a speed-versus-extent trade, and both give an answer in metres.

4. Keep the spatial index in the storage system

A query whose geometry is in a different reference system than the layer cannot use the R-tree, because the index holds bounding boxes in the stored system. Transform the query geometry once, before the query, rather than transforming every row inside it.

python
# Transform the search window into the layer's system, once
window_4326 = box(-0.11, 51.49, -0.08, 51.52)
window = transform(to_27700.transform, window_4326)
minx, miny, maxx, maxy = window.bounds

rows = conn.execute("""
    SELECT p.fid FROM parcels p
    JOIN rtree_parcels_geom r ON r.id = p.fid
    WHERE r.maxx >= ? AND r.minx <= ? AND r.maxy >= ? AND r.miny <= ?
""", (minx, maxx, miny, maxy)).fetchall()

5. Convert at the boundary, and record that you did

Exporting to 4326 for handover is a transformation like any other, which means it takes a route and deserves the same provenance record the parent guide describes. A recipient who knows the operation used can reproduce your coordinates; one who knows only the target system cannot.

Verification

sql
-- Every layer's declared system, in one place
SELECT c.table_name, c.srs_id, s.srs_name, s.organization
FROM gpkg_contents c
JOIN gpkg_spatial_ref_sys s ON s.srs_id = c.srs_id
WHERE c.data_type = 'features';
python
# An area computed in the storage system must match a geodesic area
from pyproj import Geod

geod = Geod(ellps="WGS84")
planar = projected_geom.area
geodesic = abs(geod.geometry_area_perimeter(geom_4326)[0])

assert abs(planar - geodesic) / geodesic < 0.005, (
    f"projected area differs from geodesic by "
    f"{abs(planar - geodesic) / geodesic:.2%} — wrong projection for this extent"
)

That assertion is a genuinely useful gate. A projected system appropriate to the extent agrees with the geodesic figure to a fraction of a per cent; one chosen badly — a UTM zone two zones away, a national grid used outside its area of use — diverges enough to fail immediately.

Picking a projected system from the extentThree extents and the systems that suit them. A site of a few kilometres is served by any local grid or the containing UTM zone, with distortion far below survey accuracy. A region of a few hundred kilometres wants a national grid designed for it, or the UTM zone it sits in. A territory spanning several UTM zones has no good single planar system, and either needs an equal-area projection or should be measured geodesically.a site — kilometres acrossany local grid, or the containing UTM zone · distortion far below survey accuracya region — hundreds of kilometresthe national grid designed for it, or the UTM zone it sits ina territory — several UTM zonesno good single planar system · use an equal-area projection, or measure geodesically
The extent decides this, not the data's origin — a national grid used outside its area of use is as wrong as a distant UTM zone.

Alternative Approaches or Edge Cases

Web Mercator (EPSG:3857) for storage. It is projected, so planar functions return numbers with units — but the units are badly distorted away from the equator, and area error reaches a factor of two at high latitude. It is the right storage system when the dominant use is rendering web tiles and the wrong one whenever anything is measured.

Data spanning several UTM zones. No single UTM zone is appropriate, and forcing one produces increasing distortion at the edges. Options are a national grid with a wider area of use, an equal-area projection if area is what matters, or measuring geodesically and storing in 4326.

Mixed-system containers. A GeoPackage may hold layers in different systems, which is legitimate and occasionally right — a basemap in Web Mercator beside survey data in a national grid. The cost lands on any query that joins across them, and the fix is to transform one side once rather than to redesign the container.

Troubleshooting

Areas are out by a factor of roughly 10,000

Cause: Area computed in square degrees and then treated as square metres, or the reverse. The ratio at mid-latitudes is close to that magnitude. Fix: Check what system the geometry was in when .area was called — the number is not convertible after the fact, because the conversion factor varies across the extent.

Distances are right east-to-west and wrong north-to-south

Cause: A degree of latitude is nearly constant and a degree of longitude shrinks with latitude, so a planar distance in 4326 is wrong differently in each axis. Fix: Any of the correct routes above. This asymmetry is a reliable fingerprint for geographic coordinates being measured as planar.

A projected layer renders in the wrong place on a web map

Cause: The map library assumes 4326 or 3857 and the layer is in something else. Fix: Transform on export for web consumption, or store in 3857 if rendering is the dominant use — and keep the measurement path in a suitable projected system regardless.

Frequently Asked Questions

Is EPSG:4326 ever right for storage?

Yes — when the container’s dominant use is handing it to someone else. An interchange artefact stored in 4326 needs no reference-system definition beyond the mandatory row, opens correctly in every tool, and asks nothing of the recipient. The rule that matters is not “never store 4326” but “never measure in it”, and those are separable.

How much accuracy does UTM lose at the edge of a zone?

A UTM zone is six degrees wide and its scale factor varies from about 0.9996 at the central meridian to about 1.0010 at the zone edge, so distances are off by roughly one part in a thousand at worst — a metre per kilometre. That is fine for most survey work and not fine for engineering; where it matters, a national grid designed for the area, or a local system with a suitable central meridian, is the answer.

Should the tile pyramid and the feature layers share a system?

It saves the client a transform on every frame, which on a field device is measurable. The pyramid usually commits to Web Mercator because that is what mobile map libraries expect, so matching means storing features in 3857 too — acceptable when nothing is measured from them, and a poor trade when something is. Storing features in the measurement system and letting the client transform is the more common compromise.