Migrating to gometry¶
gometry provides geometry, CRS, geodesy, grid, and spatial-index APIs for workflows commonly split across Shapely + pyproj + h3-py + s2sphere + rtree.
Use the cheatsheet to look up equivalent names, and this guide to adapt your workflow. See CRS, units & measurement, Spatial indexing & joins, and Grids & geocodes.
Coming from Shapely¶
gometry keeps Shapely's role as the in-process geometry type system — the same
seven Simple Features families, the same
predicates and constructive operations, the same WKB/WKT/GeoJSON interop and
__geo_interface__. Scalar predicates become free functions (gm.contains(a, b), vectorized over
arrays), metrics read the CRS, and the kernels are GEOS-free. Behavior differences
include validity, CRS-aware metrics, index refinement, and point-set overlay
presentation: when union, symmetric_difference, or difference with a
higher-dimensional geometry on the left retains that operand, its stored
line/polygon coordinates and Z/M pass through. See Overlay (set
operations) for the exact
cases and when to call normalize().
make_valid → repair¶
Use repair in place of make_valid.
validate()
returns a structured report with the reason and location, where Shapely splits the
job across is_valid, explain_validity, and make_valid.
from shapely import make_valid
from shapely.validation import explain_validity
if not geom.is_valid:
print(explain_validity(geom))
fixed = make_valid(geom)
import gometry as gm
bowtie = gm.from_wkt("POLYGON((0 0, 1 1, 1 0, 0 1, 0 0))")
print("valid:", bowtie.is_valid)
report = bowtie.validate()
print("reason:", report.reason)
fixed = bowtie.repair()
print("repaired:", fixed.geometry_type, "| now valid:", fixed.is_valid)
valid: False
reason: exterior ring has a self-intersection
repaired: MultiPolygon | now valid: True
CRS-aware .area, .distance, and .buffer¶
On degree-valued lon/lat data Shapely's bare
.area returns degrees², .distance returns degrees, and .buffer(100) buffers
by 100 degrees. gometry has one geom.area (and one gm.distance /
Geometry.buffer) whose answer is decided by the geometry's CRS and is native by
default — geographic shape area and length measure straight longitude/latitude
edges on the ellipsoid, point distances follow ellipsoidal geodesics, a projected
CRS uses its native linear units, and CRS-free geometry uses raw coordinate
units.
poly.area
a.distance(b)
poly.buffer(100)
import gometry as gm
area = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)
a = gm.Point(21.0, 52.0, crs=4326)
b = gm.Point(22.0, 52.0, crs=4326)
print("ellipsoidal area (m²):", round(area.area))
print("geodesic distance (m):", round(gm.distance(a, b)))
print("projected area (m²):", round(area.to_crs(area.estimate_local_crs()).area))
ellipsoidal area (m²): 30564870978
geodesic distance (m): 68677
projected area (m²): 30535837688
STRtree and prepare → SpatialIndex¶
Shapely already refines when you pass predicate= to STRtree.query.
SpatialIndex covers many geometries, with query, join,
and nearest, explicit candidates vs query(..., predicate=...), frame checks,
and missing-row semantics. There is no prepare: a geometry
probed repeatedly indexes itself after its first few probes, so a plain loop of
predicate calls gets the prepared speed.
from shapely import STRtree, prepare
tree = STRtree(geoms)
candidates = tree.query(area)
hits = tree.query(area, predicate="contains")
prepare(poly)
poly.contains(pt)
import gometry as gm
pts = gm.points([20.25, 21.0, 30.0], [51.25, 52.0, 52.0], crs=4326)
area = gm.Polygon(
[(20.0, 51.0), (22.0, 51.0), (22.0, 53.0), (20.0, 53.0), (20.0, 51.0)],
holes=[[(20.5, 51.5), (21.5, 51.5), (21.5, 52.5), (20.5, 52.5), (20.5, 51.5)]],
crs=4326,
)
idx = gm.SpatialIndex(pts)
print("candidates:", idx.candidates(area).tolist())
print("refined:", idx.query(area, predicate="contains").tolist())
print("contains first point:", gm.contains(area, pts[0]))
candidates: [0, 1]
refined: [0]
contains first point: True
Build through the geometry classes
(gm.Point, gm.Polygon, …) with explicit crs=, pack arrays with the plural
builders (gm.points, gm.polygons) or gm.GeometryArray([...]), and ingest
GeoJSON / __geo_interface__ with gm.from_geojson. Readers are verb-prefixed
free functions (gm.from_wkb, gm.from_wkt), and encoders are methods on the
geometry (geom.to_wkb(include_srid=True), geom.to_wkt(), geom.to_geojson(),
geom.to_arrow() for GeoArrow columnar interchange).
Vertex-subset ops (simplify, convex_hull, triangulation) preserve Z/M where
Shapely drops it. Hulls copy the first input vertex at a duplicated XY.
Coming from pyproj¶
pyproj owns two jobs: CRS transformation (the CRS / Transformer objects
over libPROJ) and geodesy (the Geod object). gometry keeps
PROJ as the authority backend, bundled behind the API, and
exposes transformations through to_crs or the stateless
gm.crs_transform. Point-navigation free functions
provide geodesy on geographic CRSs. Coordinate boundaries are always X/Y.
Transformer → to_crs and gm.crs_transform¶
For geometry, transform via the geometry. For raw coordinate values, use the
stateless gm.crs_transform — there is no transformer
object to construct, cache, or reuse.
t = Transformer.from_crs(4326, 3857, always_xy=True)
x2, y2 = t.transform(21.0, 52.0)
minx2, miny2, maxx2, maxy2 = t.transform_bounds(20, 51, 22, 53)
import gometry as gm
x2, y2 = gm.crs_transform(4326, 3857, 21.0, 52.0)
print("point:", round(x2, 1), round(y2, 1))
print("bounds:", [round(v) for v in gm.crs_transform_bounds(4326, 3857, (20, 51, 22, 53))])
point: 2337709.3 6800125.5
bounds: [2226390, 6621294, 2449029, 6982998]
Scalar inputs return a tuple; bulk inputs return one row-major (N, 2) / (N, 3)
NumPy matrix. transform_bounds takes the bounds as one (minx, miny, maxx, maxy)
tuple and accepts a densify= count for curved CRS edges. To pick a working
projection without hand-choosing a UTM zone, ask
geom.estimate_local_crs(). When every
coordinate has Z and the operation consumes it (geocentric, vertical, compound
CRS), to_crs transforms Z through the same pipeline; otherwise Z and M are
preserved unchanged — never invented, never silently dropped.
Geod → point navigation¶
pyproj's Geod (Karney/GeographicLib
WGS 84 ellipsoidal geodesy) maps onto gm.bearing,
point.destination, gm.point_between, and
point-to-point distances. Geographic shape area, length, and default LRS instead
measure the ellipsoid along straight longitude/latitude chart edges; unit='planar'
selects Cartesian measurement in stored coordinates. There is no Geod object
to construct. For a
non-WGS 84 ellipsoid, give the points that geographic CRS (for example
gm.Point(lon, lat, crs=4230)); every geodesic operation runs on the CRS's own
ellipsoid.
from pyproj import Geod
g = Geod(ellps="WGS84")
dist = g.inv(lon1, lat1, lon2, lat2)[2]
az = g.inv(lon1, lat1, lon2, lat2)[0]
lon2, lat2, _ = g.fwd(lon1, lat1, az, dist)
length = g.line_length(lons, lats)
area, perim = g.geometry_area_perimeter(poly)
pts = g.npts(lon1, lat1, lon2, lat2, n)
import gometry as gm
a = gm.Point(21.0, 52.0, crs=4326)
b = gm.Point(22.0, 52.0, crs=4326)
area = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)
print("distance (m):", round(gm.distance(a, b))) # geographic → geodesic
print("bearing (deg):", round(gm.bearing(a, b), 1))
print("destination:", a.destination(90.0, 1000.0))
print("midpoint:", gm.point_between(a, b, 0.5, normalized=True))
print("area (m²):", round(area.area))
print("perimeter (m):", round(area.length))
distance (m): 68677
bearing (deg): 89.6
destination: POINT (21.014560699969515 51.999999100100496)
midpoint: POINT (21.5 52.001061140753436)
area (m²): 30564870978
perimeter (m): 719739
The gm.CRS surface includes is_* classification, standards export, catalog
discovery, and coordinate operations.
Coming from h3-py & s2sphere¶
h3-py and s2sphere both turn geometry into compact cell IDs for aggregation, joins, partitioning, and coarse pre-filtering. gometry unifies H3 and S2 under typed cells and CellArray, and makes the cell-selection rule explicit. The returned cells are candidate keys; keep the source geometry for exact predicates.
polyfill → cover with an explicit rule¶
h3-py's polygon_to_cells (formerly polyfill) has an implicit containment rule.
gometry names the rule and returns a typed CellArray.
import h3
cells = h3.polygon_to_cells(h3.LatLngPoly(shell), 8) # containment rule is implicit
# Exact membership is not provided.
import gometry as gm
area = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)
cov = gm.h3_cover(area, resolution=5)
print("cells:", len(cov))
# A cover is a candidate set; retain the source for exact predicates.
print("covers:", gm.covers(area, gm.Point(21.0, 52.0, crs=4326)))
cells: 159
covers: True
cell_rule is "within" (cells the area fully owns), "center" (h3-py's
binning behavior), "overlap" (default — a complete-coverage superset, safe
candidate keys), or "bbox" (loosest). Use the free predicates
gm.covers(...) / gm.contains(...) / gm.intersects(...) against the source
geometry when exact geometry answers are needed. gm.s2_cover(geom, level=...)
selects a fixed-level cover. Omit level and pass target_cells or a level
range to adapt refinement through the same operation.
Cell-set algebra¶
Every grid family exposes its own prefixed union / intersection /
difference functions (h3_*, s2_*, geohash_*, tile_*). CellArrays
also carry compact() / uncompact(...) for multi-resolution
rollups. Family algebra preserves hierarchy-aware behavior that plain Python
sets do not.
import gometry as gm
a = gm.s2_cover(gm.box(0.0, 0.0, 2.0, 2.0, crs=4326), level=8)
b = gm.s2_cover(gm.box(1.0, 1.0, 3.0, 3.0, crs=4326), level=8)
a_present = [cell for cell in a if cell is not None]
b_present = [cell for cell in b if cell is not None]
shared = gm.s2_intersection(a_present, b_present)
rollup = shared.compact()
print("shared S2 cells:", len(shared), "compacted:", len(rollup))
h3_a = gm.h3_cover(gm.box(0.0, 0.0, 2.0, 2.0, crs=4326), resolution=5)
h3_b = gm.h3_cover(gm.box(1.0, 1.0, 3.0, 3.0, crs=4326), resolution=5)
h3_a_present = [cell for cell in h3_a if cell is not None]
h3_b_present = [cell for cell in h3_b if cell is not None]
combined = gm.h3_union(h3_a_present, h3_b_present)
print("combined H3 cells:", len(combined))
shared S2 cells: 13 compacted: 13
combined H3 cells: 115
H3Cell exposes adjacency and local indexing (neighbors, is_neighbor,
local_ij, base_cell), with resolution metadata in the H3 function family
(gm.h3_pentagons(resolution), gm.h3_base_cells()); .polygon returns one
cell's boundary polygon and cells.polygon a whole GeometryArray.
The cell classes provide point-to-cell construction, boundaries, and
compact/uncompact operations.
Coming from rtree & STRtree¶
rtree, Shapely's STRtree, and GeoPandas sjoin all accelerate spatial
predicates with a bounding-box prefilter. Shapely 2.x and GeoPandas already
refine exactly when you pass a predicate (STRtree.query(..., predicate=...),
sjoin(..., predicate=...)); rtree's intersection path is bbox-only
unless you refine yourself. gometry exposes the corresponding stages and adds
frame and row-identity checks:
- Explicit stage names —
candidates(bbox) vsquery(..., predicate=...)(exact). - Strict frame checks — CRS and epoch must match; mixed frames raise rather than silently comparing across datums.
- Missing query rows keep their input positions as missing result groups.
- Prepared caches on every geometry, built once it is probed repeatedly;
there is no separate
preparestep. - CRS-aware nearest / dwithin queries when the index frame is geographic:
point pairs use ellipsoidal geodesics, while shape queries use straight
chart-edge support with ellipsoidal separation (
unit=selects the CRS metric or planar coordinates).
gm.SpatialIndex and gm.join are the
entry points; see Spatial indexing & joins for examples.
Nearest¶
rtree (Index.nearest), Shapely (STRtree.nearest / query_nearest), and
gometry all expose nearest-neighbour APIs.
# rtree — bbox-based nearest by envelope
nearest = list(idx.nearest(point.bounds, num_results=5))
# Shapely 2.x — geometry nearest (+ distances via return_distance=)
from shapely import STRtree
tree = STRtree(geoms)
idxs = tree.query_nearest(query, return_distance=True)
import gometry as gm
pts = gm.points([21.0, 30.0, 21.2], [52.0, 52.0, 52.1], crs=4326)
target = gm.Point(21.0, 52.0, crs=4326)
idx = gm.SpatialIndex(pts)
print("nearest 2 (planar):", idx.nearest(target, k=2, unit="planar").tolist())
nearest 2 (planar): [0, 2]
idx.nearest(query, k=..., unit=...) returns the k closest indices
(num_results → k). unit="planar" selects projected or CRS-free coordinate
units; on a geographic CRS, point-pair ordering uses geodesics and shape-pair
ordering uses ellipsoidal separation over chart-edge support.
The index is mutable (insert/remove) while still sharing one frame; Shapely's
STRtree is bulk-built and immutable after construction.
Candidate selectivity¶
idx.candidates(...) exposes the bounding-box prefilter, while
idx.query(..., predicate=...) performs exact refinement. Compare their lengths
to see how much work the envelope filter leaves for the predicate:
import gometry as gm
pts = gm.points([20.25, 21.0, 30.0], [51.25, 52.0, 52.0], crs=4326)
area = gm.Polygon(
[(20.0, 51.0), (22.0, 51.0), (22.0, 53.0), (20.0, 53.0), (20.0, 51.0)],
holes=[[(20.5, 51.5), (21.5, 51.5), (21.5, 52.5), (20.5, 52.5), (20.5, 51.5)]],
crs=4326,
)
idx = gm.SpatialIndex(pts)
candidate_count = len(idx.candidates(area))
refined_count = len(idx.query(area, predicate="contains"))
print("candidate count:", candidate_count)
print("refined count:", refined_count)
candidate count: 2
refined count: 1
Spatial joins (sjoin)¶
GeoPandas sjoin
does apply the binary predicate= exactly (it is not a candidates-only API).
gm.join is the geometry-array analogue: prefilter → exact refine,
returning matched index pairs (not a joined DataFrame). gm.join provides:
- pure geometry inputs and index-pair outputs (bring your own table join);
- CRS/epoch frame checks on both sides;
- candidate and refined counts are directly observable with
len(idx.candidates(query))andlen(idx.query(query, predicate=...)); - no pandas/GeoPandas dependency for the join itself.
import geopandas
result = geopandas.sjoin(points, polygons, predicate="within")
import gometry as gm
pts = gm.points([21.0, 30.0], [52.0, 52.0], crs=4326)
polys = gm.GeometryArray([gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)])
left, right = gm.join(pts, polys, predicate="within")
print("pairs:", list(zip(left.tolist(), right.tolist(), strict=True)))
pairs: [(0, 0)]
gm.SpatialIndex(geometries) bulk-loads a balanced in-memory
R-tree over geometry envelopes — suited to
static, query-heavy workloads — while idx.insert(geometry) / idx.remove(...)
handle dynamic sets without a full rebuild. The search structure is in memory.
Pickle preserves
the live geometries and stable handles, rebuilding the search structure on load.
For portable storage, save WKB/EWKB and construct a new index.
What gometry covers¶
Data-modeling layers (GeoJSON object models and dataframe integrations such as GeoPandas) remain adjacent; gometry is the geometry engine they use.
| Package | What it does | gometry surface for the same workflow |
|---|---|---|
shapely |
planar geometry, predicates, overlay | the whole geometry surface — gm.Point, geom.buffer(...), gm.intersection, gm.contains, … |
pyproj |
CRS transforms, projection | geom.to_crs, gm.CRS, gm.crs_transform |
geographiclib |
geodesic & rhumb solutions | Use path='rhumb' with gm.distance for point distance, gm.bearing, point.destination, gm.point_between, or geom.segmentize(..., path='rhumb') for constant-bearing routes. path='geodesic' on segmentize samples the shortest ellipsoidal route. |
pyclipper |
polygon buffer / offset / clip | Geometry.buffer, Geometry.offset_curve, gm.intersection with gm.box, the overlay ops |
pygeohash family |
geohash encode/decode | gm.GeohashCell, cell.polygon, gm.geohash_cover(geom, ...) |
haversine |
point-to-point distance | gm.distance(a, b) on a geographic CRS (ellipsoidal geodesic, not the haversine sphere) |
polyline |
Google encoded polyline | gm.from_polyline, Geometry.to_polyline |
mercantile |
XYZ tiles / quadkeys | gm.Tile, gm.tile_cover(geom, ...) |
s2sphere |
S2 cells & coverings | gm.s2_*, gm.s2_cover(geom, ...) |
h3 |
H3 hexagonal cells | gm.h3_*, gm.h3_cover(geom, ...) |
openlocationcode |
plus codes | gm.pluscode_encode, gm.pluscode_polygon, gm.pluscode_shorten/recover |
utm |
lat/lon ↔ UTM | geom.estimate_local_crs() + geom.to_crs(...) |
rtree / Shapely STRtree |
bounding-box index | gm.SpatialIndex, gm.join |