NumPy ABI and Array Interop in Spatial Extensions

Almost every geospatial C-extension hands data to NumPy at some point — a raster band becomes an array, a coordinate set arrives as one, a mask goes back the other way — and that hand-off is governed by a second binary contract that sits alongside the CPython one. This guide sits under the Geospatial C-Extension Fundamentals & ABI Architecture reference and covers the NumPy C-API ABI, the buffer protocol, and the copy-versus-view decisions that determine whether a wheel built today keeps working after the next NumPy major release. It targets NumPy 1.23–2.x, GDAL 3.8+, PROJ 9.3+, Cython 3.0+, and wheels built to the abi3 contract described in C-API vs CPython ABI compatibility.

Prerequisites & Environment #

• NumPy 1.23 or newer available at build time; the build should compile against the oldest NumPy you intend to support, not the newest one installed.
• A build backend that can express a build-time NumPy pin separately from the runtime dependency — scikit-build-core, meson-python or setuptools, wired as described in mastering pyproject.toml for spatial wheels.
• GDAL 3.8+ with its Python bindings, or a binding of your own that exposes raster and coordinate buffers.
• python -c "import numpy; print(numpy.__version__, numpy.__file__)" to confirm which NumPy the build environment actually sees.

# The build-time and runtime NumPy are different questions — check both
python -c "import numpy; print('runtime', numpy.__version__)"
python -c "import numpy; print('headers', numpy.get_include())"

Core Configuration #

The NumPy C-API is versioned by a macro, and the rule is the mirror image of the CPython one: compile against old headers, run against new libraries. Building against NumPy 2 headers with no version target produces a binary that requires NumPy 2 at runtime; building against the same headers with the compatibility macro set to an older API version produces a binary that runs on both.

/* Pin the NumPy C-API version the extension compiles against.
   0x0000000F is the NumPy 1.23 API; anything newer stays unused. */
#define NPY_NO_DEPRECATED_API NPY_1_23_API_VERSION
#define NPY_TARGET_VERSION    NPY_1_23_API_VERSION
#include <numpy/arrayobject.h>

# pyproject.toml — build against the floor, depend on a range at runtime
[build-system]
requires = ["scikit-build-core>=0.9", "numpy>=1.23"]
build-backend = "scikit_build_core.build"

[project]
dependencies = ["numpy>=1.23"]

Two details in that pair matter more than they look. The build requirement must be permissive enough that a modern resolver can install some NumPy, while the compiled artifact targets the floor through the macro rather than through the pin — modern NumPy releases build forward-compatible binaries by default, and the macro makes that explicit rather than incidental. And the runtime dependency must not be capped at the NumPy you built with: an upper bound of <2 on a package that would have worked fine turns every NumPy 2 environment into an unsolvable one.

Step-by-Step Implementation #

1. Set the API floor in the source, before including the NumPy headers, so the choice is visible in the file that depends on it:

   #define NPY_NO_DEPRECATED_API NPY_1_23_API_VERSION
   #define NPY_TARGET_VERSION    NPY_1_23_API_VERSION
   #include <numpy/arrayobject.h>
   

2. Import the array API exactly once per module, in the init function, and check the result — a missing import_array() produces a segfault on the first array call rather than an error:

   PyMODINIT_FUNC PyInit__geospatial_ext(void) {
       if (_import_array() < 0) return NULL;   /* sets an exception on failure */
       return PyModule_Create(&moduledef);
   }
   

3. Point the compiler at the headers through the backend rather than a hard-coded path:

   find_package(Python REQUIRED COMPONENTS Interpreter Development.Module NumPy)
   target_link_libraries(_geospatial_ext PRIVATE Python::NumPy)
   

4. Decide copy or view per entry point and document it, using the decision table in the next section; the choice is part of your API, not an implementation detail.

Copy or View: the Decision That Shapes the API #

Every function that returns raster or coordinate data makes one of two promises. A copy hands back memory NumPy owns, independent of any dataset; a view hands back an array whose buffer is owned by native code, which is faster and imposes a lifetime rule on the caller. Getting this wrong in either direction is a real defect: a needless copy of a large raster doubles peak memory, and an undocumented view segfaults when the underlying dataset is closed.

The mechanism for a safe view is small but unforgiving. Wrap the native pointer with PyArray_SimpleNewFromData, then immediately attach an owner with PyArray_SetBaseObject so the array holds a reference to whatever keeps the memory alive — usually a capsule wrapping the dataset handle. Skipping the second call produces an array that looks correct and dereferences freed memory as soon as the dataset is closed, which is the single most common cause of intermittent crashes in hand-written spatial bindings.

npy_intp dims[2] = {rows, cols};
PyObject *arr = PyArray_SimpleNewFromData(2, dims, NPY_FLOAT32, native_buffer);
if (!arr) return NULL;
PyObject *owner = PyCapsule_New(dataset, "gdal.Dataset", release_dataset);
if (!owner) { Py_DECREF(arr); return NULL; }
if (PyArray_SetBaseObject((PyArrayObject *)arr, owner) < 0) {   /* steals owner */
    Py_DECREF(arr);
    return NULL;
}
PyArray_CLEARFLAGS((PyArrayObject *)arr, NPY_ARRAY_WRITEABLE);  /* if read-only */
return arr;

GDAL’s own bindings take the conservative route: ReadAsArray copies, so the array is independent of the dataset and safe to use after Close(). That is the right default for a general-purpose API, and it is why the memory advice in memory management in geospatial extensions treats arrays as ordinary Python objects rather than as native handles.

The Buffer Protocol as an Alternative #

Not every interop path needs the NumPy C-API. Python’s buffer protocol is part of the CPython Stable ABI, understands strides and formats, and is supported by NumPy, memoryview, array.array and most scientific libraries. An extension that consumes arrays through PyObject_GetBuffer and produces them through a small wrapper never links against NumPy at all — which removes an entire ABI contract from the build.

For the specific case that dominates geospatial extensions — passing a contiguous block of doubles to PROJ and getting it back — the buffer protocol is enough, and choosing it removes the NumPy build dependency from that code path entirely. The pattern is the same one described in releasing the GIL during coordinate transforms: acquire the buffer, copy or use in place, release it before doing anything else.

Where the NumPy C-API earns its cost is anywhere dtypes matter. A raster with a nodata value, a band stored as unsigned 16-bit with a scale factor, a mask that has to become a boolean array — all of these are dtype operations, and reimplementing them over raw buffers is both slower and more error-prone than using the machinery that exists.

Verification #

# 1. The wheel does not hard-require the NumPy it was built with
python -c "import numpy, mypkg; print(numpy.__version__, 'ok')"
# expected: prints the installed NumPy and 'ok' on both NumPy 1.x and 2.x

# 2. No NumPy symbol is dynamically required by the extension
nm -D --undefined-only mypkg/_geospatial_ext*.so | grep -ci numpy
# expected: 0 — the array API is loaded through a function table, not linked

# 3. A returned view keeps its owner alive
python - <<'PY'
import gc, mypkg
arr = mypkg.band_view("scene.tif", 1)
assert arr.base is not None, "view has no base object — it will dangle"
del mypkg; gc.collect()
print(float(arr[0, 0]))   # must not crash
PY

An extension that imports on both NumPy majors, links no NumPy symbols, and returns views with a base object has satisfied the whole contract. The second check is the one that surprises people: NumPy’s C-API is dispatched through a function pointer table populated by import_array(), so a correctly built extension has no NumPy entries in its dynamic symbol table at all.

Optimization & Edge Cases #

• Strides are not optional. A raster read with a line stride, or a band interleaved with others, is not C-contiguous. Requesting PyBUF_C_CONTIGUOUS and rejecting anything else is simplest; accepting strided input and honouring the strides is faster for the caller. Pick one and say which.
• Endianness travels with the dtype. GDAL will happily hand you big-endian data from a source file; NumPy represents that as a byte-order flag rather than converting. Code that reads the buffer directly must respect the flag or produce silently wrong numbers.
• float32 versus float64 is a real decision. Coordinates need double precision; raster values usually do not. Converting a large raster to float64 because the transform code expects doubles can double memory for no accuracy gain.
• Structured dtypes are a trap in C. A record array’s field offsets depend on alignment rules that differ between compilers and platforms. Prefer separate arrays for separate fields when the data crosses the boundary.
• NumPy 2 removed long-deprecated aliases. np.float_, np.unicode_ and several C macros are gone. Building against the floor with NPY_NO_DEPRECATED_API set catches these at compile time rather than in a user’s environment.

Troubleshooting #

ValueError: numpy.dtype size changed, may indicate binary incompatibility. The extension was compiled against a newer NumPy than the one running. Rebuild with NPY_TARGET_VERSION set to your floor, and check that the build environment is not silently installing a newer NumPy than you intended.

SystemError: <built-in function …> returned NULL without setting an error. Nearly always a missing or failed import_array(). The array API table is null, the first call dereferences it, and the failure surfaces far from the cause.

Segmentation fault when a returned array is used after the dataset is closed. A view without a base object. Attach an owner with PyArray_SetBaseObject, or return a copy — those are the only two correct options.

ImportError: numpy.core.multiarray failed to import. The runtime NumPy is older than the API version the extension was compiled for. This is the opposite of the first case and has the same fix: build against the floor.

Frequently Asked Questions #

Should a spatial package cap its NumPy dependency? #

Only when a specific incompatibility is known and documented. A speculative numpy<2 cap is the single most disruptive thing a scientific package can do to its users’ environments, because it propagates transitively and makes otherwise valid solutions unsolvable. Build forward-compatible and remove the cap.

Does using NumPy break the abi3 promise? #

No. The two contracts are independent: Py_LIMITED_API governs which CPython symbols you may use, and the NumPy API version governs which array-API entries you may call. A wheel can be cp39-abi3 and still use the NumPy C-API, because the array API is resolved through a table at import rather than through dynamic linking.

Is Cython easier than writing this by hand? #

For most projects, yes. Cython’s typed memoryviews handle the buffer acquisition, strides and release correctly, and its NumPy integration sets the API version macros for you. The rules described here still apply — they are properties of the contracts, not of the language you write against them in — but the number of places you can get them wrong drops sharply.

How do I test both NumPy majors without two environments? #

You need two environments, and they are cheap: a matrix cell that installs numpy==1.26 and one that installs the latest 2.x, both importing the built wheel and running one array round-trip. That pair catches the compile-floor mistakes immediately, and it is the only way to be confident before a user tells you.

Dtypes, Nodata and the Places Precision Leaks #

Array interop is not only a memory question; it is where a raster’s semantics either survive the boundary or quietly do not. Three properties routinely get lost, and each produces results that look plausible.

The first is nodata. GDAL records a nodata value per band as a number, not as a mask, so a band read into a plain array arrives with the sentinel embedded in the data. Code that computes a mean over that array averages the sentinel along with the measurements — and because a typical nodata value is something like -9999, the answer is wrong by an amount that depends on how much of the scene is empty. The fix is to convert at the boundary: read the nodata value alongside the data and return a masked array, or return the mask separately and document that callers must apply it. What you must not do is return a bare array and assume the caller knows.

The second is scale and offset. Many raster products store physical values as integers with a scale factor and an offset in the metadata; the raw array is not in the units the user thinks it is. Reading the band without applying them yields numbers that are internally consistent and off by a multiplicative factor, which is the hardest kind of error to notice because every downstream computation still works. Either apply the transform and return floats, or return the raw integers with the scale and offset attached — again, the choice is part of the API and belongs in the documentation rather than in the reader’s head.

The third is precision on the coordinate side. Coordinates are doubles for a reason: a float32 easting in a projected coordinate system with values in the hundreds of thousands of metres has roughly centimetre resolution, and after a couple of transformations the error is visible. Raster values are frequently fine in 32-bit, but coordinate arrays should stay in 64-bit from end to end, and an interop layer that silently downcasts is introducing error the user cannot see or correct.

# Return data with its semantics intact rather than a bare array
import numpy as np
from osgeo import gdal

def read_band(path, index=1):
    ds = gdal.Open(path)
    try:
        band = ds.GetRasterBand(index)
        arr = band.ReadAsArray()                     # copy, safe after close
        nodata = band.GetNoDataValue()
        scale = band.GetScale() or 1.0
        offset = band.GetOffset() or 0.0
        band = None
        if nodata is not None:
            arr = np.ma.masked_equal(arr, nodata)
        if scale != 1.0 or offset != 0.0:
            arr = arr * scale + offset               # promotes to float
        return arr
    finally:
        ds.Close()

A fourth issue is worth naming even though it is not a precision problem: shape conventions. GDAL indexes rasters as (row, column) with the origin at the top-left, and a multi-band read returns (band, row, column). Several other libraries in the ecosystem use (x, y) ordering or place the band axis last. Neither is wrong, and an interop layer that transposes silently — or that fails to transpose when the surrounding code expects the other convention — produces a georeferenced image that is subtly mirrored or rotated. State the convention in the docstring of every function that returns an array, and assert it in a test with a deliberately asymmetric fixture, because a square test raster will pass either way.

What Changes With NumPy 2 #

The NumPy 2 transition is the practical reason most spatial packages have had to think about this contract recently, and its lessons generalise to whatever the next major release brings.

The removals that bite hardest in geospatial code are the ones that were aliases rather than features: np.float_, np.unicode_, np.NaN and several np.lib paths that packaging code used for version checks. None of these were doing anything a modern spelling does not, and all of them appeared in code old enough that nobody was reading it. Building against the floor with NPY_NO_DEPRECATED_API turns their C equivalents into compile errors; the Python-level equivalents need a linter or a test run on the new major.

On the binary side, the important change is that extensions built against NumPy 2 headers are, by default, forward and backward compatible with a declared floor — which is the behaviour this page recommends making explicit. The failure mode people actually hit is the opposite of the one they expect: not “my NumPy 2 wheel breaks on NumPy 1” but “my wheel pinned numpy<2 and now nothing in the user’s environment can resolve”. A speculative upper bound propagates through the dependency graph and takes packages with it, which is why the guidance from the scientific ecosystem is so consistent about removing them.

There is a second-order effect specific to this domain. Because GDAL’s own Python bindings, rasterio, pyproj and shapely all interoperate through arrays, a single package with a hard NumPy cap can block an entire spatial environment from upgrading. The practical obligation for a spatial package is therefore slightly stronger than for an isolated library: build forward-compatible, test on both majors in CI, and treat a cap as an incident to be resolved rather than a safe default.

Testing both majors costs one extra matrix cell. Install the built wheel into an environment with the oldest NumPy you claim to support, and into one with the newest, and in each run a round-trip that reads a band, applies a mask and transforms a coordinate array. That single test catches the compile-floor mistakes, the dtype-promotion changes and the alias removals in one pass, which is a better return than any amount of reading release notes.

Related #

• C-API vs CPython ABI compatibility — the companion contract, and why the two are enforced in different places.
• Memory management in geospatial extensions — what a view pins, and why an array can hold hundreds of megabytes of native memory alive.
• Releasing the GIL during coordinate transforms — the buffer-copy discipline that array interop has to respect.
• Mastering pyproject.toml for spatial wheels — where the build-time and runtime NumPy requirements are declared.

Further Reading #

• The NumPy C-API documentation on ABI compatibility and NPY_TARGET_VERSION.