geocatalog

Where geocatalog sits¶

The scoping report placed geocatalog beside pipekit and geopatcher, all three above georeader and below geotoolz and xrtoolz (the picture is drawn in Report 7; the stack-wide rules are in Report 16). That still holds. geocatalog owns the spatiotemporal index end to end, has no pipekit reference at all, loads geopatcher and geocloud only lazily inside the functions that need them, and consumes georeader rather than replacing it.

What is in geocatalog¶

The mental model¶

This is the picture from the repo’s own docs. We build a GeoCatalog once, from a directory of files, a STAC search, a NASA CMR call or a hand-made list. We then hand it GeoSlices to query, and it hands back smaller catalogs (the matching rows). When we want pixels, a loader takes the slice and the matching catalog and returns a GeoTensor or an xr.Dataset. At every step the geometry, time and CRS are explicit, so there is no hidden state (Figure 1).

It works in three layers, in order:

1. Discover. Source adapters yield SourceRows, which are pre-catalog records that carry the query that produced them, whereas the local builders read filenames and metadata directly.

2. Index. A GeoCatalog Protocol implementation. Both backends expose the same surface: query, intersect, union, iter_rows, iter_slices, total_bounds and temporal_extent.

3. Materialise. Loaders take a GeoSlice plus a catalog and return a carrier.

Source layout and public surface¶

The scoping report counted 11 files and about 3,600 lines from geotoolz’s catalog module. The current package is about 14,300 lines, organised since 0.3.0 into ten namespaces by workflow step, each a thin facade over geocatalog._src (Table 1):

There is also a geocatalog command line, with build {raster, xarray, vector}, query, stats, migrate, convert and info.

GeoSlice, the unit of work¶

GeoSlice is a frozen dataclass of four fields (bounds, interval (a closed pd.Interval), resolution and crs), plus an align mode that is excluded from equality and hashing. Because it is frozen, a slice can be cached, hashed, sent across processes, or logged and replayed, and to change one we use dataclasses.replace. The repo’s docs give the reason: “silent mutation of a query that’s been logged is the worst kind of bug”.

A slice carries everything the next stage needs to fetch a chip without consulting the catalog again. The shape, height, width and transform are derived from it, and to_crs() reprojects it.

Grid alignment is the one subtlety. The bounds are corner coordinates in the slice’s CRS, $(x_0, y_0, x_1, y_1) \in \mathbb{R}^{4}$ (map coordinates, not inference variables), and the resolution is $(\Delta x, \Delta y) \in \mathbb{R}_{>0}^{2}$. The pixel counts $n_x, n_y \in \mathbb{N}$ along each axis are then

These expressions silently accept an extent that is not a whole number of pixels ((1) rounds half up), which bites at matchup boundaries, where a chip is stacked against a label raster one pixel short. aligned_shape() is the strict version, and the align= mode chooses what construction does about a residual:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import pandas as pd
from geocatalog import GeoSlice
from geocatalog.grid import is_grid_aligned

iv = pd.Interval(pd.Timestamp("2024-06-01"), pd.Timestamp("2024-06-05"),
                 closed="both")
s = GeoSlice(bounds=(0.0, 0.0, 105.0, 60.0), interval=iv,
             resolution=(10.0, 10.0), crs="EPSG:32630")   # x-extent 105 = 10.5 px

s.shape            # → (6, 11)        rounds half up, silently
s.aligned_shape()  # ValueError: x-extent=105.0 is not an integer multiple of
                   # step=10.0 (nearest is 10 -> residual -5). Use align='snap' ...

snap = GeoSlice(bounds=(0.0, 0.0, 105.0, 60.0), interval=iv,
                resolution=(10.0, 10.0), crs="EPSG:32630", align="snap")
snap.bounds        # → (0.0, 0.0, 110.0, 60.0)   rounds outward, keeps the origin

t = GeoSlice(bounds=(5.0, 0.0, 65.0, 60.0), interval=iv,
             resolution=(10.0, 10.0), crs="EPSG:32630")
is_grid_aligned(s, t, explain=True)
# → {'aligned': False, 'crs_match': True, 'x_res_match': True,
#    'y_res_match': True, 'x_origin_residual': -5.0, 'y_origin_residual': 0.0}

The modes are "off" (the default), "warn" (a GridAlignmentWarning), "error" and "snap". Snapping extends xmax right and ymin down, so the affine origin (xmin, ymax) never moves. Two slices are on the same pixel lattice when they share a CRS and resolution and their origins agree modulo the resolution, which is exactly what is_grid_aligned checks.

Why GeoSlice lives here¶

The scoping report locked in this decision, and the code bears it out:

1. The catalog is the primary producer. catalog.iter_slices(...) and CatalogDomain.slices() yield GeoSlices, and every loader consumes one.

2. The patcher does not use it. geopatcher’s unit is a Patch (data plus anchor plus neighbourhood), not a bounded request, and geopatcher does not import geocatalog.

3. A separate package would be over-engineering. GeoSlice is about 330 lines.

Table 2 shows who imports it, which so far is only geocatalog and user code.

The row schema¶

Every backend shares one row schema (Table 3):

Each kind of catalog adds its own columns:

• xarray catalogs add n_timesteps and time_var;

• vector catalogs add layer;

• STAC ingestion adds asset_key, stac_item_id and stac_collection;

• bundles add id, source, collection, assets and provenance.

The IntervalIndex is what makes temporal queries cheap. Combined with an R-tree on geometry (in memory) or bbox-column pushdown (DuckDB), query(slice) is two cheap index lookups, intersected. Times are held as naive UTC throughout, and footprints ingested from STAC or CMR that cross the antimeridian are split into a MultiPolygon at ±180°.

Two backends, one Protocol¶

Both backends implement GeoCatalog, and the Protocol is the contract, so anywhere code accepts a GeoCatalog either one works (isinstance checks work against it). Table 4 compares the two, and Figure 2 shows where each one wins on latency.

open_catalog(source, *, engine="auto", …) picks a backend for a GeoParquet artifact, DuckDB when the extra is installed and in memory otherwise, and engine="duckdb" or "memory" forces one. A remote DuckDB catalog reads only the row groups a query touches. With strict=True it raises CatalogMetadataError on incomplete metadata instead of inferring it.

Set algebra¶

query, intersect and union all return new catalogs and leave their inputs untouched, so they are safe to compose and to cache (Figure 3).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import geopandas as gpd
import pandas as pd
from shapely.geometry import box

from geocatalog import GeoSlice, intersect, open_catalog, query, union
from geocatalog.backends import InMemoryGeoCatalog
from geocatalog.storage import from_geoparquet, to_geoparquet

gdf = gpd.GeoDataFrame(
    {"filepath": ["a.tif", "b.tif"],
     "start_time": pd.to_datetime(["2024-06-01", "2024-06-10"]),
     "end_time": pd.to_datetime(["2024-06-02", "2024-06-11"])},
    geometry=[box(0, 0, 100, 100), box(50, 50, 150, 150)],   # overlap in space,
    crs="EPSG:32630",                                         # not in time
)
cat = InMemoryGeoCatalog(gdf, kind="raster")
cat.total_bounds                  # → (0.0, 0.0, 150.0, 150.0)

s = GeoSlice(bounds=(0.0, 0.0, 60.0, 60.0),
             interval=pd.Interval(pd.Timestamp("2024-06-01"),
                                  pd.Timestamp("2024-06-05"), closed="both"),
             resolution=(10.0, 10.0), crs="EPSG:32630")
list(query(cat, s).gdf["filepath"])        # → ['a.tif']  b.tif misses the window

a = InMemoryGeoCatalog(gdf.iloc[[0]].copy(), kind="raster")
b = InMemoryGeoCatalog(gdf.iloc[[1]].copy(), kind="raster")
len(intersect(a, b))                       # → 0  disjoint in time
len(intersect(a, b, spatial_only=True))    # → 1  the (50, 50, 100, 100) overlap
len(union(a, b))                           # → 2

to_geoparquet(cat, "cat.parquet")          # GeoParquet 1.1, with a bbox column
len(from_geoparquet("cat.parquet"))        # → 2
open_catalog("cat.parquet", engine="memory")   # → an InMemoryGeoCatalog

intersect(…, spatial_only=True) is the tool for pairing imagery with static labels, which have no temporal overlap by construction.

Why the loaders live here¶

The scoping report kept the loaders in geocatalog rather than georeader or geotoolz, and that held:

• They are catalog API. The pattern is “the catalog produces slices, a loader consumes them”, and both sides ship together.

• They keep georeader independent. georeader stays a reader, writer and GeoTensor library that knows nothing about catalogs.

load_raster reprojects each source into the slice’s CRS, range-reads remote files through GDAL, and has an async twin, aload_raster. load_raster_timeseries derives its time steps from the slice. The scoping report expected the loaders to be extras-gated, but at 0.5.1 only the xarray and vector loaders are lazy and the raster loader is part of the base install (Section ??).

What shipped beyond the scope¶

The scoping report described an index with builders and loaders, whereas the current package adds the pieces that turn it into a data-acquisition tool. Most of them come from the repo’s query → matchup → patch design (Table 5):

Dependencies and extras¶

The scoping report wanted an index-only install with no GDAL, made of pandas, geopandas, pyproj, shapely and rasterio (for GeoSlice’s affine transforms and windows). Loaders, DuckDB and streaming would be extras. Table 6 and Table 7 set that against what shipped.

No pipekit dependency. This held exactly, and it is the cleanest decision in the report. The catalog is not an Operator but a queryable index: pipekit composes operators that transform a carrier, whereas geocatalog finds which files match a query. We can open a catalog without pipekit, and a pipekit pipeline can read from a source that is not a catalog. The two meet in the consumer libraries, which is where the missing operators belong (Section ??).

The patcher bridge¶

CatalogDomain lets a geopatcher patcher walk a catalog without either package importing the other, and it shipped as designed in that respect. geopatcher’s Domain is a runtime_checkable Protocol with two members, crs and bounds. CatalogDomain has both, so it satisfies Domain structurally: isinstance(CatalogDomain(...), geopatcher.Domain) returns True at runtime, and neither package imports the other at module level.

1
2
3
4
5
from geocatalog.patch import CatalogDomain

dom = CatalogDomain(cat, resolution=(10.0, 10.0))
dom.bounds, dom.crs.to_epsg()    # → ((0.0, 0.0, 150.0, 150.0), 32630)
len(dom.slices())                # → 2  one GeoSlice per catalog row

field_for is the second half of the bridge, which the scoping report did not have. With the [patch] extra it hands a catalog to geopatcher as a RasterField (or a list of them), importing geopatcher lazily inside the call:

1
2
3
4
from geocatalog.patch import field_for

field = field_for(catalog, slice_, asset="visual")   # → geopatcher.RasterField
# geopatcher's SpatialPatcher then splits `field` into patches.

Migration from geotoolz¶

The scoping report planned a transitional period with three steps:

• move geotoolz.catalog and geotoolz.types.GeoSlice into geocatalog;

• have geotoolz re-export everything (including import geocatalog as catalog) for two minor versions, then deprecate the re-exports;

• keep only the pipekit-shaped operators in a new geotoolz.catalog_ops.

1
2
3
4
5
6
class CatalogPipeline(Operator):
    """ETL driver: iterate a catalog, apply an op per scene, maybe in parallel."""
    def __init__(self, catalog: GeoCatalog, op: Operator, *, n_workers=1): ...

class CatalogReader(Operator):
    """Read one GeoSlice into a GeoTensor, inline in a Sequential."""

The first step happened and the rest did not (Table 8):

The effect is that geocatalog is used on its own: we open a catalog and load slices ourselves, then hand the results to geotoolz or xrtoolz operators. The scoping report’s cross-library examples still work in that form, with today’s names:

1
2
3
4
5
6
7
8
9
10
11
from geocatalog import open_catalog
from geocatalog.load import load_raster, load_xarray

cat = open_catalog("s3://bucket/scenes.parquet")      # DuckDB if installed
for sl in cat.iter_slices(resolution=(60.0, 60.0)):   # one GeoSlice per row
    gt = load_raster(cat, sl)                         # → GeoTensor, for geotoolz
    ...

for sl in reanalysis.iter_slices(resolution=(0.25, 0.25)):
    ds = load_xarray(reanalysis, sl)                  # → xr.Dataset, for xrtoolz
    ...

CatalogPipeline is now assumed by three reports (2, 3 and this one) and by use cases 2 and 12, and it exists nowhere, which makes it the most important open item here.

The monorepo development model¶

The scoping report chose to develop everything in one “pipekit-ecosystem” monorepo with uv workspaces and publish separate packages, for four reasons:

• cross-package changes are common early on;

• GeoSlice changes ripple through every consumer;

• fixtures can be shared;

• one CI can run the integration matrix.

What exists is three monorepos, each a uv workspace with per-package releases (Table 9; Report 16 maps every distribution in them):

The operating pattern held inside each workspace. Packages resolve to local paths during development, and release-please cuts per-package tags (geotoolz-catalog-v0.5.1). Cross-package pins appear in published metadata, for example patch = [geotoolz-patcher>=0.6]. The catalog landed in the geotoolz workspace, next to the patcher and cloud I/O it works with most, rather than beside pipekit. geostack itself is where all three workspaces meet (as submodules in references/).

Effort and timing¶

The scoping report estimated a week of focused work, done in one coordinated v1.0 release alongside the pipekit and geopatcher extractions, so that users would absorb one breaking change rather than three. The record is in Table 10.

There was no single v1.0. The extraction was gradual, with breaking changes absorbed per package (geocatalog 0.3.0, geopatcher’s matching reorganisation), and geotoolz is at 0.10.2. The split unblocked what the report hoped, which is that catalog backends, schema migrations and new sources are released on their own cycle. There have been twelve geotoolz-catalog releases since July (0.1.0 to 0.5.1), each tagged on its own.

Trade-offs¶

The scoping report’s list is below, with how each item turned out.

What got better:

1. Discoverability, partly. The repo and the import name are geocatalog, but the distribution from the monorepo is geotoolz-catalog, which is not what someone searching for a catalog would find (Section ??).

2. A minimal install for narrow use: this did not happen (Section ??).

3. GeoSlice as a cross-library type: it is geocatalog.GeoSlice, but no other library imports it yet.

4. Independent backend work: this worked, as DuckDB parity, GeoParquet 1.1 pushdown and schema migrations all shipped in catalog releases.

What got worse, and the mitigations:

1. One more package to release: release-please makes it one tag.

2. Version pinning between packages: this is largely moot, because geotoolz does not depend on geocatalog yet.

3. Documentation has to live somewhere: it lives in the standalone repo (mkdocs, with the figures and the Lake Tahoe notebook), separate from the code since July.

4. The Protocol bridge is implicit: it is documented in the repo’s concepts page, and we can check it with isinstance.

Open questions¶

1. Which repo is canonical? The code lives in the geotoolz monorepo, whereas the docs, figures and notebook live in the standalone repo, frozen at 0.0.3. Either the docs move into the monorepo (and the standalone repo is archived with a pointer), or the standalone repo becomes the docs site for geotoolz-catalog.

2. geocatalog or geotoolz-catalog on PyPI? Resolved: geotoolz-catalog (Report 16 §F).

3. Should there be an index-only install? The scoped base had no georeader, and no rasterio beyond what GeoSlice needs. If the narrow use case matters, we would move georeader and the raster loader behind an extra.

4. Where do CatalogPipeline and CatalogReader live? They could live in geotoolz, as scoped, which needs geotoolz to depend on geocatalog. Alternatively they could live in an optional geocatalog[pipekit] module, which breaks the no-pipekit rule only behind an extra. Reports 2 and 3 both depend on the answer.

5. Should geotoolz and xrtoolz adopt geocatalog at all? Today the catalog is used standalone and the domain libraries never see a GeoSlice. The scoped re-exports may not be worth a dependency, whereas catalog-aware operators might be.

6. Is slices() the right bridge? CatalogDomain.slices() returns every row’s slice eagerly, whereas the scoped rows_for(anchor) was lazy and query-driven, which matters for large catalogs walked by a patcher.

7. A third backend? The scoping report floated PostGIS, but nothing in the design notes picks it up, and DuckDB over GeoParquet covers the scale cases so far.

Upcoming work¶

The follow-ups are listed in Table 11.

References¶

Sources¶

• The geocatalog repo: concepts, quickstart, recipes and the design notes on exact grid alignment and query → matchup → patch; we took the set-algebra and backend figures on this page from there.

• geotoolz-catalog 0.5.1 in the geotoolz monorepo, and its changelog.

Citations to be added: GeoParquet 1.1, DuckDB, STAC, and NASA CMR, as entries in references.bib.