Concept definitions#
This section describes scalability concepts that are important to grasp to efficiently manipulate geospatial objects.
Scalable execution relies on three complementary mechanisms:
Deferred I/O with implicit loading: Operations that update data-related metadata without loading the underlying array or geometries,
Chunked execution: Operations that process data tile-by-tile to limit memory usage,
Lazy execution: Operations whose computation is deferred until explicitly requested.
These mechanisms are often combined (e.g., Dask operations are always chunked and lazy) but are conceptually independent.
Note
The above concepts only apply to operations that interact with the underlying data arrays or geometries of objects.
Naturally, all metadata operations (e.g., accessing crs, bbox, or info()) have no effect on the array, and therefore do not trigger any loading or scalable execution.
Deferred I/O and implicit loading#
Deferred input/output refers to operations that modify only internal I/O metadata, avoiding reading the data entirely and postponing loading.
Typical examples include crop(), copy(), and translate(),
which behave similarly as Xarray’s sel(), copy(), or assign_coords().
When using the Xarray rst accessor, this behavior follows the native Xarray deferred I/O model. The Raster class implements the
same behavior so that both APIs have consistent semantics.
An important aspect of deferred I/O is that it works with both in-memory (NumPy) and scalable backends (Dask), allowing to extract parts of large files without any chunked or lazy considerations.
import geoutils as gu
# We open the dataset without Dask backend (same behaviour with Raster class)
filename_rast = gu.examples.get_path("exploradores_aster_dem")
ds = gu.open_raster(filename_rast)
# We crop the data
ds_cropped = ds.rst.icrop((0, 0, 100, 100))
# Neither input nor output raster are loaded yet
print(f"Is input raster loaded after cropping (deferred I/O)? {ds.rst.is_loaded}")
print(f"Is output raster loaded after cropping (deferred I/O)? {ds_cropped.rst.is_loaded}")
Is input raster loaded after cropping (deferred I/O)? False
Is output raster loaded after cropping (deferred I/O)? False
This behaviour pairs intrinsically with implicit loading: When an object is opened, only metadata is loaded.
Accessing data, or calling operations that require the underlying array or geometries will implicitly load the data into memory.
# Is the above raster loaded?
print(f"Is raster loaded before data operation? {ds.rst.is_loaded}")
# We compute statistics, which loads the array
ds.rst.stats()
# The raster is now loaded
print(f"Is raster loaded after data operation? {ds.rst.is_loaded}")
Is raster loaded before data operation? False
Is raster loaded after data operation? True
Chunked execution#
Chunked execution refers to processing raster data tile-by-tile instead of loading the full array into memory.
This enables out-of-core execution, allowing datasets larger than available RAM to be processed safely.
In GeoUtils, chunked execution is implemented through two backends:
Both backends read and process raster chunks sequentially, keeping peak memory usage proportional to the chunk size rather than the full dataset size. Chunked execution therefore allows GeoUtils to scale to large datasets while maintaining a predictable memory footprint. For a list of expected memory usage per operation, see the Supported operations page.
# Open raster (data is not loaded)
rast = gu.Raster(filename_rast)
# Create Multiprocessing config, output filepath optional (temporary file by default)
mp_config = gu.multiproc.MultiprocConfig(chunks=200)
# Filter raster with a gaussian in a chunked manner through Multiprocessing
rast_filt = rast.filter(method="gaussian", sigma=4, mp_config=mp_config)
# The operation happened out-of-memory in chunk-by-chunk
print(f"Temporary raster file created during operation: {rast_filt.name}")
print(f"Is input raster loaded after filtering? {rast.is_loaded}")
print(f"Is output raster loaded after filtering? {rast_filt.is_loaded}")
Temporary raster file created during operation: /tmp/tmpvzdn4kq0
Is input raster loaded after filtering? False
Is output raster loaded after filtering? False
Lazy execution#
Lazy execution refers to deferring computation until results are explicitly requested.
In GeoUtils, lazy execution is available through the Xarray rst accessor with Dask arrays.
Operations build a Dask computation graph instead of executing immediately. The computation is triggered only when required, for example when calling
compute() or when writing results to disk. It is particularly useful when chaining multiple raster operations, because intermediate results do not need to be materialized or
written/read from disk (which costs extra I/O time, often much longer than compute time).
Lazy execution always relies on chunked execution, but the reverse is not true: chunked processing can also run eagerly, as in the Multiprocessing backend.
# Open raster lazily with chunks (automatically enables Dask)
ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200})
print("Input is lazy (Dask arrays):\n")
ds
Input is lazy (Dask arrays):
<xarray.DataArray (y: 618, x: 539)> Size: 1MB
dask.array<getitem, shape=(618, 539), dtype=float32, chunksize=(200, 200), chunktype=numpy.ndarray>
Coordinates:
* y (y) float64 5kB 4.852e+06 4.852e+06 ... 4.834e+06 4.834e+06
* x (x) float64 4kB 6.272e+05 6.272e+05 ... 6.433e+05 6.433e+05
band int64 8B 1
spatial_ref int64 8B 0
Attributes:
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0If the output returned is still a Dask array, the operation was lazy.
# Interpolate 30 points from array in chunk-by-chunk
import numpy as np
rng = np.random.default_rng(seed=42)
x = rng.uniform(ds.rst.bbox.left, ds.rst.bbox.right, size=30)
y = rng.uniform(ds.rst.bbox.bottom, ds.rst.bbox.top, size=30)
ds_interp = ds.rst.interp_at_points((x, y), as_array=True)
# Result is still lazy
print("Result is still lazy after raster interpolation:\n")
ds_interp
Result is still lazy after raster interpolation:
|
||||||||||||||||
We can materialize it with compute(), which now returns a NumPy array:
ds_interp = ds_interp.compute()
ds_interp
array([1235.2104 , 1413.3643 , 1015.6273 , 1179.5364 , 1426.9363 ,
1244.6381 , 1087.3882 , 1007.41235, 2838.135 , 1281.0834 ,
1633.231 , 1150.2612 , 1218.9752 , 1043.8933 , 1304.0006 ,
1203.2421 , 1137.6985 , 3652.7976 , 1174.9348 , 1176.458 ,
939.13885, 1645.7578 , 1315.4735 , 1046.2162 , 1109.4808 ,
1414.9052 , 1317.2467 , 1293.689 , 2714.8916 , 1073.3386 ],
dtype=float32)