Raster machine learning with xarray and rioxarray

This tutorial shows how to build a raster classification workflow with xarray, rioxarray, and scikit-learn. It follows the same broad steps as the Pyspatialml landcover tutorial: open raster predictors, align them, sample labelled pixels, fit a classifier, and write a classified raster.

The example uses the Landsat bands and labelled-pixel raster bundled with Pyspatialml. You can use the example data without using the Raster class.

Install the additional packages

Xarray and rioxarray are not required dependencies of Pyspatialml. Install them along with Dask to run the chunked prediction example:

python -m pip install xarray rioxarray dask

Open and align the predictor rasters

Each input GeoTIFF becomes a two-dimensional DataArray. rio.reproject_match puts every layer on the same grid as the first layer. This step is useful even when the inputs are expected to be aligned: it makes the target grid explicit.

import numpy as np
import os
import tempfile
import xarray as xr
import rioxarray
import matplotlib.pyplot as plt

from pyspatialml.datasets import nc

paths = {
    "band_10": nc.band1,
    "band_20": nc.band2,
    "band_30": nc.band3,
    "band_40": nc.band4,
    "band_50": nc.band5,
    "band_70": nc.band7,
}

layers = {
    name: rioxarray.open_rasterio(path, masked=True).squeeze("band", drop=True)
    for name, path in paths.items()
}

reference = next(iter(layers.values()))
layers = {
    name: layer.rio.reproject_match(reference)
    for name, layer in layers.items()
}

# A Dataset gives each predictor a name. Stack the variables on a feature axis
# for scikit-learn, which expects rows=samples and columns=features.
predictor_ds = xr.Dataset(layers)
predictors = predictor_ds.to_array(dim="feature")
predictors = predictors.rio.write_crs(reference.rio.crs)

predictors

Plot one predictor to inspect its spatial pattern. Xarray uses the raster’s coordinate labels for the axes.

fig, ax = plt.subplots(figsize=(7, 5))
predictor_ds["band_30"].plot(ax=ax, cmap="Greens")
ax.set_title("Landsat band 3")
plt.show()

If your rasters already share the same CRS, transform, width, and height, you can skip reprojection and open them directly. Reprojecting large rasters has a cost, so check the grids before applying it to a large collection.

Build training data from labelled pixels

Open the labelled-pixel raster and match it to the predictor grid. The example assumes zero represents unlabelled pixels, as in the bundled dataset. Convert the labelled cells to a table, remove missing predictor values, then fit a scikit-learn model.

labels = rioxarray.open_rasterio(nc.labelled_pixels, masked=True)
labels = labels.squeeze("band", drop=True).rio.reproject_match(reference)

# Extract the labelled pixels and their predictor values by 
# stacking both the labelled pixel raster and the predictor dataset into a single DataFrame.
training = xr.Dataset({"class": labels, **layers})
training = training.where(training["class"] > 0)
training = training.to_dataframe().dropna()

feature_names = list(paths)
X = training[feature_names].to_numpy()
y = training["class"].to_numpy()

from sklearn.ensemble import RandomForestClassifier

model = RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)
model.fit(X, y)

RandomForestClassifier(n_estimators=200, n_jobs=-1, random_state=42)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

RandomForestClassifier

?Documentation for RandomForestClassifieriFitted

Parameters

	n_estimators n_estimators: int, default=100 The number of trees in the forest. .. versionchanged:: 0.22 The default value of ``n_estimators`` changed from 10 to 100 in 0.22.	200
	n_jobs n_jobs: int, default=None The number of jobs to run in parallel. :meth:`fit`, :meth:`predict`, :meth:`decision_path` and :meth:`apply` are all parallelized over the trees. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary <n_jobs>` for more details.	-1
	random_state random_state: int, RandomState instance or None, default=None Controls both the randomness of the bootstrapping of the samples used when building trees (if ``bootstrap=True``) and the sampling of the features to consider when looking for the best split at each node (if ``max_features < n_features``). See :term:`Glossary <random_state>` for details.	42
	criterion criterion: {"gini", "entropy", "log_loss"}, default="gini" The function to measure the quality of a split. Supported criteria are "gini" for the Gini impurity and "log_loss" and "entropy" both for the Shannon information gain, see :ref:`tree_mathematical_formulation`. Note: This parameter is tree-specific.	'gini'
	max_depth max_depth: int, default=None The maximum depth of the tree. If None, then nodes are expanded until all leaves are pure or until all leaves contain less than min_samples_split samples.	None
	min_samples_split min_samples_split: int or float, default=2 The minimum number of samples required to split an internal node: - If int, then consider `min_samples_split` as the minimum number. - If float, then `min_samples_split` is a fraction and `ceil(min_samples_split * n_samples)` are the minimum number of samples for each split. .. versionchanged:: 0.18 Added float values for fractions.	2
	min_samples_leaf min_samples_leaf: int or float, default=1 The minimum number of samples required to be at a leaf node. A split point at any depth will only be considered if it leaves at least ``min_samples_leaf`` training samples in each of the left and right branches. This may have the effect of smoothing the model, especially in regression. - If int, then consider `min_samples_leaf` as the minimum number. - If float, then `min_samples_leaf` is a fraction and `ceil(min_samples_leaf * n_samples)` are the minimum number of samples for each node. .. versionchanged:: 0.18 Added float values for fractions.	1
	min_weight_fraction_leaf min_weight_fraction_leaf: float, default=0.0 The minimum weighted fraction of the sum total of weights (of all the input samples) required to be at a leaf node. Samples have equal weight when sample_weight is not provided.	0.0
	max_features max_features: {"sqrt", "log2", None}, int or float, default="sqrt" The number of features to consider when looking for the best split: - If int, then consider `max_features` features at each split. - If float, then `max_features` is a fraction and `max(1, int(max_features * n_features_in_))` features are considered at each split. - If "sqrt", then `max_features=sqrt(n_features)`. - If "log2", then `max_features=log2(n_features)`. - If None, then `max_features=n_features`. .. versionchanged:: 1.1 The default of `max_features` changed from `"auto"` to `"sqrt"`. Note: the search for a split does not stop until at least one valid partition of the node samples is found, even if it requires to effectively inspect more than ``max_features`` features.	'sqrt'
	max_leaf_nodes max_leaf_nodes: int, default=None Grow trees with ``max_leaf_nodes`` in best-first fashion. Best nodes are defined as relative reduction in impurity. If None then unlimited number of leaf nodes.	None
	min_impurity_decrease min_impurity_decrease: float, default=0.0 A node will be split if this split induces a decrease of the impurity greater than or equal to this value. The weighted impurity decrease equation is the following:: N_t / N * (impurity - N_t_R / N_t * right_impurity - N_t_L / N_t * left_impurity) where ``N`` is the total number of samples, ``N_t`` is the number of samples at the current node, ``N_t_L`` is the number of samples in the left child, and ``N_t_R`` is the number of samples in the right child. ``N``, ``N_t``, ``N_t_R`` and ``N_t_L`` all refer to the weighted sum, if ``sample_weight`` is passed. .. versionadded:: 0.19	0.0
	bootstrap bootstrap: bool, default=True Whether bootstrap samples are used when building trees. If False, the whole dataset is used to build each tree.	True
	oob_score oob_score: bool or callable, default=False Whether to use out-of-bag samples to estimate the generalization score. By default, :func:`~sklearn.metrics.accuracy_score` is used. Provide a callable with signature `metric(y_true, y_pred)` to use a custom metric. Only available if `bootstrap=True`. For an illustration of out-of-bag (OOB) error estimation, see the example :ref:`sphx_glr_auto_examples_ensemble_plot_ensemble_oob.py`.	False
	verbose verbose: int, default=0 Controls the verbosity when fitting and predicting.	0
	warm_start warm_start: bool, default=False When set to ``True``, reuse the solution of the previous call to fit and add more estimators to the ensemble, otherwise, just fit a whole new forest. See :term:`Glossary <warm_start>` and :ref:`tree_ensemble_warm_start` for details.	False
	class_weight class_weight: {"balanced", "balanced_subsample"}, dict or list of dicts, default=None Weights associated with classes in the form ``{class_label: weight}``. If not given, all classes are supposed to have weight one. For multi-output problems, a list of dicts can be provided in the same order as the columns of y. Note that for multioutput (including multilabel) weights should be defined for each class of every column in its own dict. For example, for four-class multilabel classification weights should be [{0: 1, 1: 1}, {0: 1, 1: 5}, {0: 1, 1: 1}, {0: 1, 1: 1}] instead of [{1:1}, {2:5}, {3:1}, {4:1}]. The "balanced" mode uses the values of y to automatically adjust weights inversely proportional to class frequencies in the input data as ``n_samples / (n_classes * np.bincount(y))`` The "balanced_subsample" mode is the same as "balanced" except that weights are computed based on the bootstrap sample for every tree grown. For multi-output, the weights of each column of y will be multiplied. Note that these weights will be multiplied with sample_weight (passed through the fit method) if sample_weight is specified.	None
	ccp_alpha ccp_alpha: non-negative float, default=0.0 Complexity parameter used for Minimal Cost-Complexity Pruning. The subtree with the largest cost complexity that is smaller than ``ccp_alpha`` will be chosen. By default, no pruning is performed. See :ref:`minimal_cost_complexity_pruning` for details. See :ref:`sphx_glr_auto_examples_tree_plot_cost_complexity_pruning.py` for an example of such pruning. .. versionadded:: 0.22	0.0
	max_samples max_samples: int or float, default=None If bootstrap is True, the number of samples to draw from X to train each base estimator. - If None (default), then draw `X.shape[0]` samples irrespective of `sample_weight`. - If int, then draw `max_samples` samples. - If float, then draw `max_samples * X.shape[0]` unweighted samples or `max_samples * sample_weight.sum()` weighted samples. .. versionadded:: 0.22 .. versionchanged:: 1.9 Float `max_samples` is relative to `sample_weight.sum()` instead of `X.shape[0]` for weighted samples.	None
	monotonic_cst monotonic_cst: array-like of int of shape (n_features,), default=None Indicates the monotonicity constraint to enforce on each feature. - 1: monotonically increasing - 0: no constraint - -1: monotonically decreasing If monotonic_cst is None, no constraints are applied. Monotonicity constraints are not supported for: - multiclass classifications (i.e. when `n_classes > 2`), - multioutput classifications (i.e. when `n_outputs_ > 1`). The constraints hold over the probability of the positive class. Read more in the :ref:`User Guide <monotonic_cst_gbdt>`. .. versionadded:: 1.4	None

Fitted attributes

Name	Type	Value
classes_ classes_: ndarray of shape (n_classes,) or a list of such arrays The classes labels (single output problem), or a list of arrays of class labels (multi-output problem).	ndarray[float32](6,)	[1.,3.,4.,5.,6.,7.]
estimator_ estimator_: :class:`~sklearn.tree.DecisionTreeClassifier` The child estimator template used to create the collection of fitted sub-estimators. .. versionadded:: 1.2 `base_estimator_` was renamed to `estimator_`.	DecisionTreeClassifier	DecisionTreeClassifier()
estimators_ estimators_: list of DecisionTreeClassifier The collection of fitted sub-estimators.	list	[DecisionTreeC...te=1608637542), DecisionTreeC...te=1273642419), DecisionTreeC...te=1935803228), DecisionTreeC...ate=787846414), ...]
estimators_samples_ estimators_samples_: list of arrays The subset of drawn samples (i.e., the in-bag samples) for each base estimator. Each subset is defined by an array of the indices selected. .. versionadded:: 1.4	list	[array([2380, ..., dtype=int32), array([ 686, ..., dtype=int32), array([1366, ..., dtype=int32), array([1591, ..., dtype=int32), ...]
feature_importances_ feature_importances_: ndarray of shape (n_features,) The impurity-based feature importances. The higher, the more important the feature. The importance of a feature is computed as the (normalized) total reduction of the criterion brought by that feature. It is also known as the Gini importance. Warning: impurity-based feature importances can be misleading for high cardinality features (many unique values). See :func:`sklearn.inspection.permutation_importance` as an alternative.	ndarray[float64](6,)	[0.17,0.19,0.13,0.25,0.13,0.12]
n_classes_ n_classes_: int or list The number of classes (single output problem), or a list containing the number of classes for each output (multi-output problem).	int	6
n_features_in_ n_features_in_: int Number of features seen during :term:`fit`. .. versionadded:: 0.24	int	6
n_outputs_ n_outputs_: int The number of outputs when ``fit`` is performed.	int	1

Plot the labelled pixels to inspect their coverage and class values before looking at the model output.

fig, ax = plt.subplots(figsize=(7, 5))
labels.plot(ax=ax, cmap="tab10", add_colorbar=True)
ax.set_title("Training labels")
plt.show()

This example uses the labelled raster directly. For vector training points, reproject the GeoDataFrame to the raster CRS, then sample each layer at the point coordinates (for example with DataArray.sel(x=..., y=..., method="nearest")). For polygons, choose an explicit pixel-inclusion rule and aggregate or sample the pixels inside each geometry.

Predict in Dask chunks

Scikit-learn estimators accept two-dimensional arrays, while the predictors are organized as (feature, y, x). The function below moves the feature axis to the end for each spatial chunk, predicts only pixels with valid values, and returns a two-dimensional chunk. Dask delays the work until the output is written or computed.

def predict_chunk(block):
    # With the feature dimension declared as a core dimension, xarray passes
    # each block as (y, x, feature).
    height, width, n_features = block.shape
    pixels = block.reshape(-1, n_features)
    valid = np.isfinite(pixels).all(axis=1)

    result = np.full(pixels.shape[0], np.nan, dtype="float32")
    if valid.any():
        result[valid] = model.predict(pixels[valid]).astype("float32")
    return result.reshape(height, width)

chunked = predictors.chunk({"feature": -1, "y": 512, "x": 512})
prediction = xr.apply_ufunc(
    predict_chunk,
    chunked,
    input_core_dims=[["feature"]],
    output_core_dims=[[]],
    dask="parallelized",
    vectorize=False,
    output_dtypes=[np.float32],
)
prediction = prediction.rio.write_crs(reference.rio.crs)
prediction = prediction.rio.write_transform(reference.rio.transform())

output = prediction.fillna(0).astype("int16")
output_path = os.path.join(tempfile.gettempdir(), "landcover_prediction.tif")
output.rio.to_raster(
    output_path,
    dtype="int16",
    nodata=0,
    compute=True,
)

The example writes to the system temporary directory so rendering the tutorial does not leave generated rasters in the documentation source folder. Set output_path to a project path when you want to keep the result.

Plot the classified result. Pixels without valid predictor values are shown as the nodata background.

fig, ax = plt.subplots(figsize=(7, 5))
prediction.plot(ax=ax, cmap="tab10", robust=True, add_colorbar=True)
ax.set_title("Predicted landcover")
plt.show()

Chunk sizes are workload-dependent. Each task needs enough memory for its feature block and the estimator’s temporary arrays. If the prediction does not fit comfortably in memory, use smaller spatial chunks.

How this maps to Pyspatialml

Workflow step	xarray and rioxarray approach
Open and name raster layers	`open_rasterio` and `xr.Dataset`
Align grids	`rio.reproject_match`
Extract labelled pixels	Mask and convert a Dataset to a DataFrame
Apply a scikit-learn model	Reshape feature and spatial dimensions; use `apply_ufunc` with Dask chunks
Write a georeferenced output	`rio.to_raster`

This approach gives you xarray’s labeled dimensions and its Dask integration, while making the model-to-raster steps explicit in your code. Pyspatialml provides a more direct raster-stack API for sampling, extraction, and windowed estimator prediction.

Which approach should you use?

Use xarray and rioxarray when you want labeled dimensions, time-series or multidimensional analysis, or Dask-based computation across chunks. This is also a good fit when your project already uses xarray. Use Pyspatialml when you prefer its raster-stack interface and want sampling, extraction, and scikit-learn prediction packaged as raster methods. Both approaches work with GeoTIFFs and NumPy arrays, so you can also use Pyspatialml for a focused step inside a larger xarray workflow.