Data scientists encounter GIS when spatial data shows up in a project: a dataset with latitude/longitude columns, a shapefile from a government portal, a GeoJSON file someone dropped in the repo. The first few operations usually work. Then something breaks silently — a spatial join returns an empty DataFrame, a buffer produces geometry the size of a continent, or exported shapefiles arrive at the lab with every field name truncated.
These aren’t random bugs. They stem from a small set of spatial concepts that GIS practitioners learn early and that data science curricula skip entirely. This post covers the ones that cost the most time.
CRS: The Concept That Causes Most Silent Failures
Coordinate Reference System (CRS) is the framework that maps coordinates to actual locations on Earth. Different CRS values produce different coordinate numbers for the same physical point — and mixing them causes incorrect results with no error message.
The EPSG codes that matter
- EPSG:4326 — WGS84, geographic. Degrees, lon/lat. What GPS outputs, what GeoJSON requires by spec. Units are degrees, not meters.
- EPSG:3857 — Web Mercator. What Google Maps, Leaflet, and OpenStreetMap tile servers use. Units are meters but the projection severely distorts areas above ~60° latitude. Greenland appears the same size as Africa; Africa is 14× larger.
- EPSG:326xx — UTM zones.
EPSG:32618= UTM Zone 18N (eastern US),EPSG:32610= Zone 10N (Pacific coast). Meters, low distortion within a 6°-wide strip. The right CRS for measuring distances and areas anywhere in the world. - EPSG:4269 — NAD83. Looks visually identical to WGS84 but coordinates differ by up to 2 meters. Enough to break any precision work.
set_crs() vs to_crs() — the single most common mistake
gdf = gpd.read_file("data.shp") # stored in EPSG:4326
gdf = gdf.set_crs(epsg=3857) # WRONG: relabels, does not transform
result = gdf.buffer(1000) # "1000 degrees" ≈ 111,000 km radius
set_crs() only rewrites metadata. It tells geopandas “these coordinates are in 3857” without moving them. to_crs() mathematically transforms every coordinate. No exception is raised in either case — the buffer produces nonsense silently.
The correct workflow for any spatial measurement operation:
gdf = gpd.read_file("data.shp") # EPSG:4326
gdf_utm = gdf.to_crs("EPSG:32618") # transform to UTM Zone 18N
buffered = gdf_utm.buffer(1000) # real 1 km buffer in meters
result = buffered.to_crs("EPSG:4326") # back to geographic if needed
Axis order: PROJ 6+ broke a lot of code
EPSG:4326 officially defines axis order as (latitude, longitude). But most software, including geopandas, stores coordinates as (x, y) = (longitude, latitude). Raw pyproj respects the official axis order unless you override it:
from pyproj import Transformer
# Without always_xy=True, PROJ 6+ swaps axes per the official CRS definition
t = Transformer.from_crs("EPSG:4326", "EPSG:32618")
x, y = t.transform(40.748, -73.985) # lat, lon input → axes may swap silently
# Correct: always_xy=True forces (longitude, latitude) input regardless of CRS definition
t = Transformer.from_crs("EPSG:4326", "EPSG:32618", always_xy=True)
x, y = t.transform(-73.985, 40.748) # lon, lat input → x, y output in meters
Code written against PROJ 4 that passed (lat, lon) without always_xy=True silently produces wrong output in PROJ 6+.
Spatial joins with mismatched CRS
geopandas.sjoin() does not raise an exception when the two GeoDataFrames are in different CRS. It returns an empty DataFrame or geometrically incorrect results:
# This fails silently — returns 0 rows even if every point is inside a polygon
result = gpd.sjoin(points_4326, polygons_3857, predicate="within")
# Check CRS before every join
print(points.crs, polygons.crs) # both must match
polygons = polygons.to_crs(points.crs) # align before joining
result = gpd.sjoin(points, polygons, predicate="within")
Also watch the predicate direction:
gpd.sjoin(points, polygons, predicate="contains") # "does each point contain a polygon?" → always False
gpd.sjoin(points, polygons, predicate="within") # correct: "is each point within a polygon?"
Projection for Spatial Analysis
Never buffer in degrees
At 40°N, one degree of longitude is ~83 km. At 60°N it is ~55 km. At the equator it is ~111 km. A buffer expressed in degrees produces different physical sizes at every latitude, and there’s no way to interpret the output number as a real-world distance.
# WRONG — "1000 meter" buffer in degree-based CRS
gdf_wgs84.buffer(1000) # 1000° radius ≈ 111,000 km at equator
# CORRECT
gdf_utm = gdf_wgs84.to_crs("EPSG:32618") # US East Coast
gdf_utm.buffer(1000) # 1 km in real meters
Web Mercator for analysis is always wrong
EPSG:3857 area distortion scales as 1/cos²(lat). At 60° latitude areas are inflated 4×. Comparing areas, computing centroids, or running kernel density estimation in Web Mercator produces quantitatively incorrect results even if the map looks right. Use it for display tiles only; measure in UTM or a local projected CRS.
UTM zone boundaries
A dataset spanning multiple UTM zones cannot be accurately analyzed in a single zone. A pipeline from UTM Zone 43N into Zone 44N accumulates meters of positional error at the zone boundary when everything is forced into one zone. Use a national equal-area projection (EPSG:5070 for the contiguous US, EPSG:3338 for Alaska) for cross-zone analysis.
The workflow rule: Store in EPSG:4326. Visualize on the web in EPSG:3857. Measure in the appropriate UTM zone or local projected CRS.
Shapefile Gotchas
Shapefiles are the most-sent GIS format and the most trap-filled. The format dates to 1994 and its constraints show it.
10-character field name limit — silently truncated
gdf["population_density"] = gdf["population"] / gdf["area_km2"]
gdf.to_file("output.shp")
# Field name silently becomes "population_" or "populati_1"
# No warning, no error
When two truncated names collide, ESRI’s convention appends a serial number (populat_1, populat_2). The original names are gone. This breaks downstream code that references the field by name. The 10-character field truncation problem is documented in the GDAL shapefile driver spec — there is no workaround except using GeoPackage or GeoParquet.
Size limits — both geometry and attributes
The .shp geometry file and the .dbf attribute file each have a 2 GB size limit independently. A shapefile with many attributes can exhaust the DBF limit while the geometry file is tiny. The total row count cap is ~2 billion rows, but most software struggles long before that.
Additional constraints:
- Max 255 attribute columns
- Text fields max 254 characters
- No datetime with timezone
- No boolean type (stored as integer 0/1)
- Field names case-insensitive in DBF —
Pop_Totalandpop_totalmay not survive as distinct columns
Missing companion files
A shapefile is 4 files minimum: .shp (geometry), .shx (spatial index), .dbf (attributes), .prj (CRS). When the .prj file is missing, geopandas reads the file with gdf.crs is None. Calling to_crs() on a CRS-less GeoDataFrame raises a CRSError that traces to the missing file, not to the data itself. Always check print(gdf.crs) immediately after read_file().
GeoJSON Gotchas
Coordinate order is [longitude, latitude]
This is the most common GeoJSON bug:
// WRONG — places a point somewhere off the Antarctic coast
{"type": "Point", "coordinates": [40.7484, -73.9857]}
// CORRECT — Empire State Building
{"type": "Point", "coordinates": [-73.9857, 40.7484]}
RFC 7946 is explicit: coordinates are [longitude, latitude]. But many APIs and mapping libraries use lat/lon order internally, and the conversion is wrong half the time because both orderings are plausible numbers. A point placed at [lat, lon] instead of [lon, lat] lands on a real place — just the wrong one — with no validation error.
The crs key is illegal in RFC 7946
The original GeoJSON spec (2008) allowed a "crs" member specifying the projection. RFC 7946 (2016) removed it entirely. GeoJSON must be WGS84. If you receive a GeoJSON file with a "crs" key, it was written before 2016 and the projection metadata in that key is not reliable — the actual coordinate values may or may not be WGS84. Treat it as unknown projection and verify against a known point before using.
Antimeridian crossing
A polygon spanning the 180th meridian (the antimeridian) is broken for Shapely:
# A polygon crossing the antimeridian from 170°E to -170°E (190°W)
# Shapely interprets this as spanning the entire globe westward
poly = Polygon([(170, 0), (-170, 10), (-170, -10), (170, -10)])
poly.area # ~340 × 20 = ~6800 square degrees — the whole globe minus a sliver
Operations like contains, intersection, and area on antimeridian-crossing geometries return wrong results in Shapely without raising any exception. Split along the antimeridian line before processing. The gis-file-conversion-silent-failures post covers this and four other categories of conversion failures that produce valid-looking but wrong output.
Polygon winding order
RFC 7946 requires exterior rings to be counterclockwise (right-hand rule) and interior rings (holes) to be clockwise. Most renderers tolerate violations, but some topology validators and area algorithms are sensitive. If you’re generating GeoJSON programmatically, enforce winding order:
from shapely.geometry import mapping
from shapely.ops import orient
oriented = orient(polygon, sign=1.0) # 1.0 = counterclockwise exterior
geojson = mapping(oriented)
Topology: Why Shapefiles and GeoJSON Silently Store Wrong Geometry
Neither shapefiles nor GeoJSON have any concept of topology. Adjacent polygons each store a full independent copy of their shared boundary. Edit one vertex in polygon A, and the neighboring polygon B still holds the old coordinate — producing a sliver gap or overlap that is invisible in normal rendering but causes errors in spatial operations.
# This returns non-empty geometry that should be empty
polygons.geometry.unary_union.difference(polygons.geometry.unary_union)
# If there are slivers, this returns small artifact polygons
PostGIS ST_IsValid() and QGIS’s built-in Topology Checker detect these. For automated validation in Python, use shapely.validation.explain_validity(geom) on each feature.
TopoJSON solves this by encoding shared edges as arcs referenced by index. Simplifying one arc updates all polygons that reference it. Simplifying GeoJSON polygons independently reintroduces slivers that TopoJSON was designed to eliminate. For web display where polygon boundaries must align exactly, TopoJSON is more reliable than GeoJSON for topologically sensitive datasets.
Python Library Pitfalls
Deprecated pyproj syntax that silently works wrong
# PROJ 4 syntax — still accepted but broken in PROJ 6+
Proj(init='EPSG:4326') # raises: "Invalid projection: +init=..."
Proj('+proj=utm +zone=18 +datum=WGS84') # deprecated +proj string
# Correct
Proj('EPSG:4326') # EPSG code directly
CRS.from_epsg(4326) # explicit constructor
pyproj.Proj ignores datum shifts
Proj handles unit conversions within a single datum but silently ignores datum shifts. The difference between NAD27 and WGS84 is ~27 meters in the eastern US — meaningful for utility mapping or parcel boundary work but invisible if you use Proj where Transformer.from_crs() is required.
# Misses NAD27 → WGS84 datum shift of up to 27 meters
Proj("EPSG:26718").transform(x, y) # wrong
# Correct: Transformer handles datum shift
t = Transformer.from_crs("EPSG:26718", "EPSG:32618", always_xy=True)
x_new, y_new = t.transform(x, y) # applies full geodetic transformation
PROJ string storage is lossy
crs = CRS.from_epsg(4326)
crs.to_proj4() # "+proj=longlat +datum=WGS84 +no_defs +type=crs"
# loses authority information, geographic metadata, datum ensemble info
When saving CRS to a database or file header, use crs.to_wkt() or store the EPSG code directly. PROJ4 strings were the standard until 2018 and remain in widespread use, but they are a lossy representation that causes ambiguity in datum-sensitive applications.
Null Island in your data
Many ETL pipelines store geocoding failures as (0.0, 0.0) — a real coordinate in the Gulf of Guinea. Spatial joins and aggregations on global datasets silently include these phantom points:
# Filter Null Island before any spatial operation
gdf = gdf[~((gdf.geometry.x == 0) & (gdf.geometry.y == 0))]
Format Choices for Analysis vs Storage
For working datasets that live on your machine, GeoJSON is a poor choice: it stores coordinates as text (slow to parse, no spatial indexing), requires full file scan for any query, and has no column-level typing. Use GeoPackage (.gpkg) for working data — single file, SQLite-backed, spatial index built in, field names up to 63 characters. Use GeoParquet for large datasets that go into notebooks or Pandas pipelines. The GeoParquet vs GeoPackage vs Shapefile decision tree covers when each format wins.
When you need to hand off to systems that only accept shapefiles, convert at the last step and accept the field name truncation — don’t work in shapefiles internally. If you’re converting between formats regularly, especially from GeoJSON to shapefile, GeoConvert handles the field name mapping, CRS preservation, and geometry validation automatically, which matters when you’re processing files you didn’t create.
Checklist Before Any Spatial Operation
print(gdf.crs)— is it what you expect? Is it None?- If measuring (buffer, area, distance):
to_crs("EPSG:326xx")first - Before
sjoin(): verify both DataFrames share the same CRS - After
to_file("out.shp"): read back and verify field names survived - Check for Null Island:
(x==0) & (y==0)rows - For polygons from external sources: run
shapely.validation.explain_validity()on anything used in spatial predicates