Usage and good practices

Usage and good practices#

GeoUtils supports scalable execution for most raster and point cloud operations, as well as some vector operations.

This summary page introduces scalability concepts, demonstrates usage, and lists good practices. For more details on scalability concepts, support or implementations, refer to the next pages.

GeoUtils relies on two execution backends:

  • Dask, through the rst and pc Xarray accessors and the vct and pc GeoPandas accessors (lazy and chunked execution),

  • Multiprocessing, through its Raster, Vector and PointCloud objects (chunked execution only).

Both backends use the exact same chunked logic, and yield identical results as in-memory operations. Lazy refers to deferred execution using Dask, while chunked refers to processing raster tiles or dataframe partitions separately to limit memory usage. For details on scalability concepts, see the Concept definitions page.

As a rule of thumb:

  • Use Dask to work on Xarray and GeoPandas objects through our accessors rst, vct and pc, and if you want to chain several operations lazily.

  • Use Multiprocessing to work with our Raster, Vector and PointCloud objects, and if you are fine with intermediate writing/reading between steps.

  • Use standard in-memory execution to work efficiently on small datasets, which is possible even if those were loaded from larger files (keep in mind that crop() or subsample() only load requested data, even without Dask/Multiprocessing).

These choices involve both memory and time. Parallel backends have initialization and scheduling costs, so they are often slower for small inputs. See Performance for controlled comparisons and Profiling to measure the behavior of your own data and hardware.

Note

GeoUtils currently targets scalable CPU execution. However, as many of our numerical operations rely on NumPy, SciPy or Numba, those are planned to be linked to their GPU counterparts (CuPy and Numba CUDA).

Using Dask through accessors#

With Dask, raster operations are both chunked and lazy. This behavior is enabled by opening a file with the chunks argument, which returns an Xarray or GeoPandas object backed by Dask arrays.

import geoutils as gu

filename_rast = gu.examples.get_path("exploradores_aster_dem")
ds = gu.open_raster(filename_rast, chunks={"x": 200, "y": 200})
ds

<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.0
    overview_factors:          []
    overview_selected_level:   -1
    overview_selected_factor:  1
    block_shape:               (3, 539)

GeoUtils, through the rst accessor, automatically detects the Dask input and switches to a chunked implementation for the given operation, for example to reproject() to a different resolution:

# Change output resolution
out_res = (ds.rst.res[0] * 2, ds.rst.res[1] / 2)

# Reproject lazily and out-of-memory
ds_reproj = ds.rst.reproject(
    res=out_res,
    resampling="bilinear",
)
ds_reproj

<xarray.DataArray (y: 1236, x: 270)> Size: 1MB
dask.array<concatenate, shape=(1236, 270), dtype=float32, chunksize=(200, 200), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 10kB 4.852e+06 4.852e+06 ... 4.834e+06 4.834e+06
  * x            (x) float64 2kB 6.272e+05 6.273e+05 ... 6.433e+05 6.433e+05
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:             Area
    scale_factor:              1.0
    add_offset:                0.0
    overview_factors:          []
    overview_selected_level:   -1
    overview_selected_factor:  1
    block_shape:               (3, 539)

The resulting raster remains lazy (i.e. a Dask array). Computation (reading input files and applying all chained operations) only happens when later explicitly requested with compute().

For a raster output, one typically wants to write to file while avoid loading it in-memory:

ds_reproj.rst.to_file("reproj_rast.tif", compute=True)

This triggers computation and writes the raster chunk-by-chunk to the output file.

Or, the output can be chained with another operation. For example, we can extract a small random subsample of the reprojected raster:

# Subsample lazily and out-of-memory
sub_ds = ds_reproj.rst.subsample(
    subsample=5000,
)
sub_ds

<xarray.DataArray 'b1' (point: 5000)> Size: 20kB
dask.array<getitem, shape=(5000,), dtype=float32, chunksize=(5000,), chunktype=numpy.ndarray>
Coordinates:
    x        (point) float64 40kB 6.295e+05 6.349e+05 ... 6.341e+05 6.29e+05
    y        (point) float64 40kB 4.838e+06 4.843e+06 ... 4.844e+06 4.845e+06
Dimensions without coordinates: point
Attributes:
    crs:         PROJCRS["WGS 84 / UTM zone 18S",BASEGEOGCRS["WGS 84",ENSEMBL...
    geometry_z:  False

The output point cloud is again lazy, although selecting a fixed sample size can run preliminary counting tasks. Calling compute() returns an in-memory DataArray:

sub_ds.compute()

<xarray.DataArray 'b1' (point: 5000)> Size: 20kB
array([2562.606 , 1260.0782, 2305.1023, ..., 1169.9612, 1366.7542,
       1431.2871], shape=(5000,), dtype=float32)
Coordinates:
    x        (point) float64 40kB 6.295e+05 6.349e+05 ... 6.341e+05 6.29e+05
    y        (point) float64 40kB 4.838e+06 4.843e+06 ... 4.844e+06 4.845e+06
Dimensions without coordinates: point
Attributes:
    crs:         PROJCRS["WGS 84 / UTM zone 18S",BASEGEOGCRS["WGS 84",ENSEMBL...
    geometry_z:  False

Lazy execution is particularly useful to chain several operations without ever loading the full raster into memory or writing any file to disk.

Using Multiprocessing through GeoUtils objects#

Chunked execution can also be enabled on the Raster object using Multiprocessing. Multiprocessing performs chunked execution (only loads chunks in-memory), but is not lazy (it runs immediately).

rast = gu.Raster(filename_rast)
rast

Raster(
  data=not_loaded; shape on disk (1, 618, 539); will load (618, 539)
  transform=| 30.00, 0.00, 627175.00|
            | 0.00,-30.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32718
  nodata=-9999.0)

By passing a MultiprocConfig configuration to an operation, the out-of-memory behaviour is triggered. The configuration requires chunks, and can optionally define the output file and cluster to use (defaults to temporary instances for both, see Multiprocessing configuration).

from geoutils.multiproc import MultiprocConfig

# Optional: specify output file (otherwise a temporary file is created)
mp_config = MultiprocConfig(chunks=200, outfile="reproj_rast.tif")

# Reproject out-of-memory, reading and writing chunk-by-chunk
rast_reproj_mp = rast.reproject(
    res=out_res,
    resampling="bilinear",
    mp_config=mp_config,
)
rast_reproj_mp

Raster(
  data=not_loaded; shape on disk (1, 1236, 270); will load (1236, 270)
  transform=| 60.00, 0.00, 627175.00|
            | 0.00,-15.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32718
  nodata=-9999.0)

If the output is a Raster, it is written to disk out-of-memory, and the returned object is a Raster of that file without data loaded. This keeps syntax consistent with in-memory code, and allows to easily chain operations.

Point outputs are written to a GeoPackage and returned as an unloaded PointCloud.

# Subsample out-of-memory and return an unloaded point cloud
point_config = MultiprocConfig(chunks=200, outfile="subsample.gpkg")
samp_rast_mp = rast_reproj_mp.subsample(
    subsample=5000,
    mp_config=point_config,
)
samp_rast_mp

PointCloud(
  ds=               b1                geometry
       0     1217.582031  POINT (642595 4840445)
       1     1113.102539  POINT (627715 4851995)
       2     1201.984741  POINT (643195 4836995)
       3     1543.900146  POINT (632575 4847795)
       4     1283.852051  POINT (634795 4851680)
       ...           ...                     ...
       4995  1301.445557  POINT (628855 4844330)
       4996  3275.488525  POINT (628495 4837610)
       4997  1194.393066  POINT (627175 4846265)
       4998  1452.622314  POINT (630595 4846910)
       4999  2392.880127  POINT (632395 4835015)
       
       [5000 rows x 2 columns]
  crs=EPSG:32718
  bbox=BoundingBox(left=np.float64(627175.0), bottom=np.float64(4833560.0), right=np.float64(643255.0), top=np.float64(4852085.0)))

This backend is convenient when working directly with Raster objects and performing step-by-step processing.

Good practices with chunked and lazy operations#

  • Downsample or use file overviews when opening/plotting/analyzing data (see Downsampling and overviews on how to do this in GeoUtils),

  • Choose chunk sizes large enough to reduce scheduling overhead, but small enough to fit comfortably in memory,

  • Check that your data files have disk chunksizes (otherwise loads all) and use a multiple of them as in-memory chunking (see Formats and cloud storage to handle disk chunks in GeoUtils),

  • If memory is the limiting factor for you, use a single-threaded scheduler through Dask (dask.config.set(scheduler='single-threaded')) or sequential chunk processing through multiprocessing (default cluster) to avoid parallel chunk being loaded,

  • If speed is the limiting factor for you, use parallelized processes through Dask (see Dask scheduler configuration) or through a multiprocessing cluster,

  • Insert breakpoints (for example by writing intermediate results to disk) to prevent building overly large Dask graphs.

For more guidance on chunk sizing and performance, see the Dask array best practices.

Finally, ensure an operation is supported in chunks in GeoUtils or elsewhere. The full description of supported methods is available on the Supported operations page.