Arrays & performance¶

Many gometry operations on one geometry also accept a whole column and keep the per-row traversal in native code.

Scalar and array calls share names and semantics, with results shaped to match their inputs:

Call shape	Input	Output
`geom.area`	scalar `Geometry`	`float`
`arr.area`	`GeometryArray`	`float64` ndarray
`geom.buffer(d)`	scalar	`Polygon` / `MultiPolygon`
`arr.buffer(d)`	array	`GeometryArray[Polygon \| MultiPolygon]`
`gm.contains(poly, pt)`	scalar pair	`bool`
`gm.contains(poly, pts)`	scalar × array	`bool_` ndarray
`idx.query(scalar)`	one geometry	`int64` ndarray
`idx.query(arr)`	`GeometryArray`	`Groups` of `numpy.typing.NDArray[numpy.int64]`
`geom.coords`	geometry or array	flattened `Coordinates` view; `np.asarray(view, dtype=np.float64)` exports a matrix and `view.row_index` gives source rows

GeometryArray¶

gm.GeometryArray builds a packed, typed container of geometries that every vectorized function consumes. It behaves like a sequence (indexing, iteration, len) while keeping geometry storage in Rust-owned objects rather than a Python list. It carries one CRS for the column. coordinate_axes reports the layout of each row; common_coordinate_axes reports one shared layout for present rows, or None when those layouts differ:

import gometry as gm

arr = gm.GeometryArray([
    gm.box(20.0, 51.0, 22.0, 53.0, crs=4326),
    gm.box(29.0, 51.0, 32.0, 53.0, crs=4326),
])
print("len:", len(arr))
print("crs:", arr.crs)
print("coordinate_axes:", arr.coordinate_axes)
print("common_coordinate_axes:", arr.common_coordinate_axes)
print("element types:", arr.geometry_type)
print("combined bounds:", arr.total_bounds)
print("areas (ellipsoidal m²):", repr(arr.area))

len: 2
crs: EPSG:4326
coordinate_axes: ['XY', 'XY']
common_coordinate_axes: XY
element types: ['Polygon', 'Polygon']
combined bounds: (20.0, 51.0, 32.0, 53.0)
areas (ellipsoidal m²): array([3.05648710e+10, 4.58473065e+10])

GeometryArray.total_bounds is the combined (minx, miny, maxx, maxy) tuple, or None for an empty/fully-empty array. For geographic arrays, longitude bounds are periodic in the CRS's native frame, so minx > maxx denotes an interval crossing the longitude seam.

Packed constructors¶

The packed constructors take parallel coordinate arrays and build a GeometryArray directly. gm.points accepts optional z= and m= arrays, so you can construct XY, XYZ, XYM, or XYZM arrays in one call — the bulk analogue of gm.Point(..., z=..., m=...):

import gometry as gm

pts = gm.points([21.0, 22.0, 23.0], [52.0, 52.5, 53.0], crs=4326)
print("2D per row:", pts.coordinate_axes, "shared:", pts.common_coordinate_axes, "len", len(pts))

pts_z = gm.points([21.0, 22.0], [52.0, 52.5], z=[100.0, 200.0], crs=4979)
print("XYZ per row:", pts_z.coordinate_axes, "shared:", pts_z.common_coordinate_axes)

2D per row: ['XY', 'XY', 'XY'] shared: XY len 3
XYZ per row: ['XYZ', 'XYZ'] shared: XYZ

The sibling column factories gm.line_strings, gm.polygons, and the rest skip the per-row Python object loop the same way; IO parsers (from_wkb, from_wkt, from_geojson) also return packed arrays directly.

Per-row axes and missing rows¶

coordinate_axes is row-aligned and uses None for a missing row. common_coordinate_axes reports one shared layout for formats that require it:

import gometry as gm

mixed = gm.GeometryArray([
    gm.from_wkt("POINT Z (1 2 3)"),
    None,
    gm.from_wkt("POINT M (1 2 9)"),
])
print("per row:", mixed.coordinate_axes)
print("shared:", mixed.common_coordinate_axes)
print("missing:", mixed.is_missing.tolist())

per row: ['XYZ', None, 'XYM']
shared: None
missing: [False, True, False]

Indexing, slicing, and element types¶

Indexing with an integer returns a single typed Geometry; slicing returns a new GeometryArray. Iteration yields typed elements:

import gometry as gm

arr = gm.GeometryArray([gm.Point(i, i) for i in range(5)])
print("element:", type(arr[0]).__name__, arr[0].x)
print("slice:", type(arr[1:4]).__name__, len(arr[1:4]))
print("strided:", [g.x for g in arr[::2]])

element: Point 0.0
slice: GeometryArray 3
strided: [0.0, 2.0, 4.0]

The container is generic over its element type: gm.points(...) and gm.GeometryArray([pt, ...]) are typed GeometryArray[Point], so a type checker knows arr[0] is a Point and for g in arr yields Point. Heterogeneous arrays widen to GeometryArray[Geometry]. Operations narrow too, as far as the element kind guarantees: this point array's arr.buffer(...) is a GeometryArray[Polygon], arr.area is a float64 ndarray while geom.area is float, and so on across the whole surface.

GeometryArray, CellArray, H3VertexArray, and H3EdgeArray share one row-selection rule. An integer or integer-like object implementing __index__ (not a boolean), or a zero-dimensional NumPy integer array, selects one element. A slice, one-dimensional integer array, integer sequence, or boolean mask selects a container of the same family. NumPy selectors with rank greater than one raise TypeError because these row containers have one dimension.

Type checkers preserve the result when an ndarray's rank is known; if its annotation uses a variable-length tuple of integer dimensions, the result is a scalar-or-container union. If NumPy types the shape as Any, annotate the rank when you need a narrower result type.

The public classes and cross-grid protocols can be used directly in your own annotations. Accepted string tokens are already narrowed in every gometry signature, so IDEs and type checkers flag misspellings at the call site:

import gometry as gm
def cell_area(cell: gm.Cell) -> float:
    return cell.area

print(round(cell_area(gm.H3Cell("871fb4662ffffff"))))

NumPy bridge and preview¶

to_numpy explicitly materializes an object ndarray of geometry handles; NumPy's array protocol also supports np.asarray(arr, dtype=object). gm.GeometryArray is the inverse for object ndarrays and __geo_interface__ rows. Coordinate columns live on .coords (.x/.y/.z/.m are read-only float64 ndarrays):

import numpy as np
import gometry as gm

arr = gm.points([0.0, 1.0, 2.0], [0.0, 1.0, 2.0])
obj = arr.to_numpy()
back = gm.GeometryArray(obj)
print("round-trip len:", len(back), "| x column:", repr(np.asarray(back.coords.x)))

round-trip len: 3 | x column: array([0., 1., 2.])

Like scalar geometries, a geographic array previews itself with lonboard when the optional visualization extra is installed. _repr_html_svg renders the capped SVG grid used without it, bounded so a large array never renders thousands of elements:

4 geometries

Vectorized calls & broadcasting¶

Array methods and binary free functions accept arrays and return arrays. A scalar polygon tested against an array of points is evaluated in one native pass, returning a read-only numpy.ndarray that shares Rust-owned storage with no copy:

import gometry as gm

poly = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)
points = gm.points([21.0, 30.0, 21.5], [52.0, 52.0, 52.1], crs=4326)

mask = gm.contains(poly, points)
print(type(mask).__name__, repr(mask))

ndarray array([ True, False,  True])

Result types follow the operation:

Kind	Return shape	Examples
Numbers	`float64` / `int64` / `uint64` ndarray	`area`, `distance`, `line_locate`, spatial keys
Masks	`bool_` ndarray	`contains`, `intersects`, `within`, …
Bounds	`(n, 4)` `float64` ndarray	`bounds` rows (`nan` for empty)
Id pairs	`(left, right)` pair of `int64` ndarrays	`join`, `self_join`
Ragged matches	`Groups` (CSR)	array-form `index.query`
Geometries	`GeometryArray`	`buffer`, `to_crs`, …; overlays via `& / \|`
Coordinates	`Coordinates` view	`geom.coords`, packed point columns
Witness pairs	`(left, right)` pair of point `GeometryArray` columns	`nearest_points`

Text/bytes and ragged diagnostics stay ordinary Python lists because they have no single dense ndarray representation.

gometry adopts NumPy's length-preserving broadcasting cases and refuses everything else — there is no implicit Cartesian product. Exactly three operand shapes are allowed:

Left	Right	Result	Meaning
scalar	array (n)	array (n)	test the scalar against each element
array (n)	scalar	array (n)	test each element against the scalar
array (n)	array (n)	array (n)	pairwise — element i vs element i
array (n)	array (m), n≠m	`GeometryError`	refused — use `index`/`join`

Scalar × array¶

import gometry as gm

poly = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)
points = gm.points([21.0, 30.0], [52.0, 52.0], crs=4326)
print(gm.contains(poly, points).tolist())

[True, False]

Equal-length pairwise¶

When both sides are arrays of the same length, the operation is pairwise, not a cross product: result i compares left i with right i.

import gometry as gm

polys = gm.GeometryArray([
    gm.box(20.0, 51.0, 22.0, 53.0, crs=4326),
    gm.box(29.0, 51.0, 31.0, 53.0, crs=4326),
])
points = gm.points([21.0, 30.0], [52.0, 52.0], crs=4326)

print("pairwise contains:", gm.contains(polys, points).tolist())

pairwise contains: [True, True]

Mismatched lengths are refused¶

Two non-scalar arrays of different lengths raise instead of broadcasting into an n×m grid. The exception identifies the mismatch:

import gometry as gm

polys = gm.GeometryArray([
    gm.box(20.0, 51.0, 22.0, 53.0, crs=4326),
    gm.box(29.0, 51.0, 31.0, 53.0, crs=4326),
])
points = gm.points([21.0, 30.0, 40.0], [52.0, 52.0, 52.0], crs=4326)

try:
    gm.contains(polys, points)
except gm.GeometryError as e:
    print("GeometryError:", e)

GeometryError: GeometryArray operands must have the same length, got 2 and 3

Many-to-many uses an index or join

gm.join(polys, points, predicate="contains") runs a bbox prefilter plus exact refine in Rust and returns matching pairs. See Indexing & joins.

Coming from Shapely? See Migrating. Shapely 2.x already follows NumPy broadcasting and raises on incompatible 1-D shapes (e.g. lengths 2 and 3) rather than a Cartesian product. gometry's differences that matter here are attaching a CRS to arrays and using the same strict equal-length / scalar-vs-array broadcast rules without n×m expansion.

Missing rows¶

None is a missing row; a typed empty geometry is a real geometry value. For example, POINT Z EMPTY retains its point kind and XYZ axes even though it has no coordinates. GeometryArray carries missing rows natively using a validity model. A dense array stores packed columns without a missing-row marker when no row is missing. This distinction is preserved across arrays, dataframe handoff, Arrow, GeoParquet, and most vectorized operations. Batch ingest (from_wkt / from_wkb / from_geojson / GeometryArray([...])) keeps native GeoArrow typing when some rows are missing rather than downgrading the whole column to WKB:

The columns after geometry are is_missing, area, and contains(window):

import gometry as gm
import numpy as np

geoms = gm.GeometryArray([gm.Point(1, 1), None, gm.from_wkt('POINT EMPTY')])
window = gm.box(0, 0, 2, 2)
areas = geoms.area
contains = gm.contains(window, geoms)
coordinate_view = geoms.coords
values = np.asarray(coordinate_view, dtype=np.float64)
indexes = coordinate_view.row_index

for row, geometry in enumerate(geoms.to_wkt()):
    print(row, geometry, geoms.is_missing[row], areas[row], contains[row])
print('coordinate values:', values.tolist())
print('coordinate indexes:', indexes.tolist())

0 POINT (1 1) False 0.0 True
1 None True nan False
2 POINT EMPTY False 0.0 False
coordinate values: [[1.0, 1.0]]
coordinate indexes: [0]

values is the flattened coordinate matrix and indexes repeats each source row for its coordinates. len(geoms) counts logical rows; len(geoms.coords) counts vertices. Missing and typed-empty rows both contribute no vertices, but geoms.is_missing distinguishes a missing row from an empty geometry.

One rule set everywhere: predicates → False, measures → NaN, row-aligned geometry results and ragged groups stay missing, aggregates skip, exports emit None. Fixed-width columns that cannot hold None keep their ndarray shape with sentinels — for example bounds use NaN, topological dimensions use -1, and some curve-key columns use 0. Those sentinels can collide with valid values, so is_missing is authoritative. Flattening and explode-style operations (parts, triangulations, sampled points) keep the missing source row as a None group; its flat payload has no placeholder geometry.

Missing rows enter through any constructor or importer (GeometryArray([g, None]), from_wkt/from_wkb with None entries, RFC 7946 "geometry": null features, Arrow validity bitmaps, geopandas/pandas/polars nulls) and survive every row-aligned operation, slice, fancy indexing (arr[[0, 2]]), pickle, and Arrow round-trip.

Coordinates flatten present rows

arr.coords is a vertex stream, not a row-aligned array. Missing rows contribute no vertices. Its row_index property maps each coordinate back to its source row, or call drop_missing() first when a dense coordinate-only workflow is clearer.

drop_missing() and fill_missing(value) convert back to dense. fill_missing accepts one scalar geometry or a row-aligned GeometryArray fill source; indexes and GeometryCollection construction require dense input. Drop or fill missing rows before passing them to a boundary that cannot represent them.

CellArray¶

CellArray is a homogeneous sequence of typed cells. Hierarchy ops on a cell object return one — H3Cell.children(), .neighbors, and the S2/geohash/tile mirrors materialize CellArray[H3Cell] (or the matching cell type) instead of a bare Python list. Accessors mirror the scalar cell name: cells.token, cells.center, and cells.area stay singular while returning columnar/list results. Cell geometry is explicit through accessors such as .polygon, .center, and .to_polygon(); cells do not carry an OGC geometry_type because they are grid identifiers, not geometry objects.

Construct one directly from a non-empty homogeneous iterable of cell objects and gometry infers the grid type: gm.CellArray([cell_a, cell_b]). Raw integer IDs, string tokens, NumPy arrays, buffers, and empty inputs have no reliable type evidence, so make the grid explicit with type=gm.H3Cell (or the corresponding S2/geohash/tile class). Mixed cell classes and an explicit type that disagrees with a cell object are errors.

Construction preserves input order and duplicates; CellArray is a sequence, not a set.

Bulk cell factories preserve missing input rows as None. Missing rows are included by slicing, reversing, and fancy indexing; None in cells, cells.count(None), and cells.index(None) follow normal Python sequence semantics, while cell/id searches ignore missing rows.

Row-preserving CellArray operations keep one result per logical input row: hierarchy predicates return False for either missing side, geometry accessors retain the geometry-array missing mask, numeric results use NaN, parent retains its mask, ragged hierarchy results contain an empty row, and token contains None. Set-like compact, uncompact, and to_polygon operate on present cells only.

import gometry as gm

cell = gm.H3Cell(gm.Point(21.0, 52.0, crs=4326), resolution=6)
children = cell.children()
print(type(children).__name__, len(children), type(children[0]).__name__)
print("ids:", [c.id for c in children[:3]])

CellArray 7 H3Cell
ids: [608537035086495743, 608537035103272959, 608537035120050175]

H3 topological vertices and directed edges are not cells: H3Cell.vertices returns H3VertexArray and H3Cell.edges returns H3EdgeArray. Those arrays keep the same sequence and uint64 id-column ergonomics but only expose valid topology operations such as vertices.point, edges.origin, edges.destination, edges.line, and edges.length.

Bulk point→cell construction uses the prefixed plural builders — gm.h3_cells, gm.s2_cells, and the geohash/tile twins. Grid cover factories return CellArray for one geometry and Groups of CellArray for geometry arrays; hierarchy transforms on that set may return CellArray when the engine can keep a packed id column.

For exact membership, keep the source geometry separately and use free predicates such as gm.covers or gm.contains.

Groups (ragged results)¶

Some operations return one variable-length result per input row — spatial index queries against an array of probes, or multi-hit predicate scans. Groups stores those results as a flat .values payload plus .offsets that bound each row. The payload keeps its actual type: an int64 or float64 ndarray, a GeometryArray, or a CellArray.

Construct reusable groups with Groups(values, offsets). The flat values array determines the payload kind even when empty, and offsets start at zero and end at len(values). Numeric payloads must be one-dimensional int64 or float64 ndarrays; integer IDs remain exact. Geometry and cell payloads keep their inner missing values. Use missing= for missing groups, whose offset windows must be empty. This distinguishes a missing group from a present group containing missing values.

For example, gm.Groups(np.array([1.0, 3.0, 5.0]), [0, 2, 2, 3]) represents three station groups: [1.0, 3.0], [], and [5.0]. Passing numeric groups to GeometryArray.line_interpolate samples each route at its own explicit stations.

import gometry as gm

pts = gm.points([21.0, 30.0, 40.0], [52.0, 52.0, 52.0], crs=4326)
zones = gm.GeometryArray([
    gm.box(20.0, 51.0, 22.0, 53.0, crs=4326),
    gm.box(29.0, 51.0, 31.0, 53.0, crs=4326),
])
idx = gm.SpatialIndex(pts)
groups = idx.query(zones, predicate="contains")

print(type(groups).__name__, "rows:", len(groups))
print("row 0 ids:", groups[0].tolist())
print("row 1 ids:", groups[1].tolist())
print("offsets:", groups.offsets.tolist())

Groups rows: 2
row 0 ids: [0]
row 1 ids: [1]
offsets: [0, 1, 2]

Groups.values is the logical flat payload window; contiguous source windows can share backing, while gathered values may be materialized. A slice whose first row is source row zero can share its offsets; a slice starting later returns rebased offsets. Groups.counts is materialized from adjacent offsets.

Index hits use read-only int64 value and offset arrays; geometry-valued ragged results wrap per-row GeometryArray slices.

GeometryArray.parts and other ragged array operations return one group per source row. A missing source row is None; a present empty geometry is an empty group. Groups.is_missing identifies the missing rows while offsets continues to map each source row to its flat payload. Groups.to_list() preserves missing rows as None.

Coordinate extraction is the flat numeric exception: geom.coords exposes one vertex stream; convert it with NumPy for a matrix and read view.row_index when source-row indexes are needed.

Scalar vs array index query

A scalar idx.query(geom, predicate=...) returns a dense int64 ndarray. Pass an array of query geometries and you get Groups — one result row per element. See Spatial indexing.

Batching choices¶

Gometry keeps geometry kernels in Rust and exposes packed, vectorized entry points.

Vectorize instead of looping¶

A Python loop calls the scalar API once per object; a vectorized call accepts the whole batch and returns one aligned result.

import gometry as gm

poly = gm.box(20.0, 51.0, 22.0, 53.0, crs=4326)
lons = [20.0, 20.5, 21.0, 30.0]
lats = [51.0, 51.5, 52.0, 10.0]

sample = [gm.Point(lo, la, crs=4326) for lo, la in zip(lons, lats)]
loop_result = [gm.contains(poly, g) for g in sample]

pts = gm.points(lons, lats, crs=4326)
vectorized_result = gm.contains(poly, pts)

print("same answer:", loop_result == vectorized_result.tolist())
print("rows evaluated:", len(vectorized_result))

same answer: True
rows evaluated: 4

Repeated predicates on one geometry¶

When you test one geometry against many others, call the predicate directly. A geometry probed repeatedly builds an edge index after its first few probes and keeps it, so later calls, scalar or vectorized, reuse it. For two arrays, prefer the vectorized predicates or a join.

Index or join for many-to-many work¶

A nested loop over two collections makes O(N·M) predicate calls. The SpatialIndex and gm.join replace it with a bbox prefilter plus exact refine in Rust. Cost depends on index construction, query count, bounding-box selectivity, candidate count, and predicate complexity; adversarial overlapping boxes can still approach a Cartesian product.

import gometry as gm

points = gm.points([21.0, 30.0, 21.5], [52.0, 52.0, 52.1], crs=4326)
areas = gm.GeometryArray([
    gm.box(20.0, 51.0, 22.0, 53.0, crs=4326),
    gm.box(29.0, 51.0, 31.0, 53.0, crs=4326),
])

pairs = gm.join(points, areas, predicate="within")
print("matched pairs:", pairs)

matched pairs: (array([0, 1, 2]), array([0, 1, 0]))

For repeated queries against a fixed set, build the index once and reuse it; see Indexing & joins.

Prefer `equals_identical` for array value identity¶

For array value identity, use equals_identical. It is the vectorized form of scalar ==; frame differences return False, and dimensional empties remain axes-sensitive.

Ingest in bulk, and use the direct parser form¶

Batch readers accept a whole iterable in one call. Direct string or bytes input is scanned in Rust instead of being expanded into a Python object tree first.

import gometry as gm

wkb = [gm.Point(i * 0.01, i * 0.02, crs=4326).to_wkb() for i in range(5)]

loop_result = [gm.from_wkb(b) for b in wkb]

batch_result = gm.from_wkb(wkb)

print("same rows:", [g.to_wkt() for g in loop_result] == batch_result.to_wkt())
print("batch result:", batch_result)

same rows: True
batch result: <GeometryArray[Point] len=5>

The same holds for every reader. from_wkb/from_wkt/from_geojson all take an iterable and return a packed GeometryArray in one pass; the geodesic/CRS batches (destination, point_between, CRS.geodesic*, gm.crs_transform_bounds) are vectorized the same way, so a whole column of ingest and downstream transform stays in Rust:

import numpy as np

import gometry as gm

wkb = gm.points([0.0, 1.0, 2.0], [50.0, 51.0, 52.0], crs=4326).to_wkb(
    include_srid=True
)
geoms = gm.from_wkb(wkb)

destinations = geoms.destination(
    np.array([90.0, 180.0, 270.0]),
    np.array([1000.0, 2000.0, 3000.0]),
)

bounds = gm.crs_transform_bounds(
    4326,
    3857,
    [(0, 50, 1, 51), (1, 51, 2, 52)],
)

print("ingested:", geoms)
print("destinations:", destinations.geometry_type)
print("bounds shape:", bounds.shape)

ingested: <GeometryArray[Point] len=3 EPSG:4326>
destinations: ['Point', 'Point', 'Point']
bounds shape: (2, 4)

For very large files, stream at the file/database layer into bounded batches, then call gometry once per batch — that keeps memory bounded without falling back to per-row parser overhead.

Pass raw GeoJSON to the parser

gm.from_geojson(text) parses the JSON in Rust and skips the discarded properties entirely. Calling json.loads first and passing the resulting dict forces a full Python object tree to be built and re-walked before the native parser sees it. Pass the raw str/bytes (a file's contents, a polars/duckdb binary column) and let the Rust reader do the scan.

Ingest	Per-item or expanded form	Batched/direct form
GeoJSON	`from_geojson(json.loads(s))`	`from_geojson(s)` — `str`/`bytes`
Many points	`GeometryArray([Point(x, y) for ...])`	`gm.points(xs, ys)`
WKB/WKT	`[from_wkb(b) for b in blobs]`	`from_wkb(blobs)` — one batch
Features + attributes	`json.loads` + a Python feature loop	`gm.from_geojson(text, features=True)`

When you need the per-feature properties/id alongside geometry, pass features=True to gm.from_geojson. It parses the whole FeatureCollection in Rust and returns a validated Features(geometries, properties, ids) record — one crossing instead of a Python loop.

Crossing the dataframe / Arrow boundary¶

Per-row interop conversion changes a columnar boundary into Python object work. Bulk converters accept whole columns.

pandas: arr.to_pandas() stores the native array behind a concrete extension dtype; gm.from_pandas() shares that storage zero-copy. It also imports object Series of geometry-like rows; those rows are converted rather than shared.
Polars: arr.to_polars() / gm.from_polars() encode and decode a whole WKB/EWKB Binary Series in native batched calls. The boundary does not require PyArrow.
GeoPandas: gm.from_geopandas() / arr.to_geopandas() convert a GeoSeries in one vectorized WKB step. Avoid rebuilding a GeoSeries element by element ([g.wkb for g in series]); that per-row hop into Shapely is the per-item path the vectorized helpers exist to replace.
Arrow: array.to_arrow() and gm.from_arrow exchange whole columns through GeoArrow. Homogeneous geometry uses native layouts; mixed types can use WKB fallback. Importing an Arrow capsule provider does not require PyArrow.

For dataframe geometry columns, see DataFrames. For portable GeoArrow interchange or GeoParquet storage, see Arrow & storage.