Statistics

GeoUtils provides implementations to derive statistical estimators relevant to geospatial data on rasters and point clouds, with efficient reduction for large datasets.

Several types of operations are supported:

• Global and grouped statistics derive distributional estimators by bins or categories,

• Spatial statistics (variography) derives spatial correlation estimators as a function of distance.

Zonal statistics are included as grouped statistics where vector features define the categories, with support for fractional polygon-raster intersections.

Any custom estimator can be passed, with certain estimators listed below optimized for larger-than-memory samples (both for grouped and spatial estimators).

Note

Statistical operations accept datasets with different CRS or spatial support (e.g. binning a raster with a point cloud), using efficient co-sampling techniques under-the-hood. For more details on those, see the Sampling feature page.

All methods on this page accept combinations of out-of-memory rasters, point clouds and vectors. They normally follow the chunks or partitions used when opening each input; see Chunked strategies for advanced controls.

Show the code for opening example files

Hide the code for opening example files

import geoutils as gu
import numpy as np
import pandas as pd

# Open an elevation raster and glacier outlines covering the same region
ds = gu.open_raster(gu.examples.get_path("exploradores_aster_dem"))
glaciers = gu.open_vector(gu.examples.get_path("exploradores_rgi_outlines"))

Global and grouped statistics

Global statistics derive estimators on the entire raster or point cloud, while grouped statistics apply the same estimators to continuous bins or discrete categories, with potentially several groupings at once.

Note

Global and grouped statistics use NumPy reductions and optionally Flox’s groupby_reduce(). GeoUtils adds subsampling, weights support, and fractional vector-raster intersection support.

Estimators

GeoUtils accepts any custom estimator, and provides the built-ins below.

“Mergeable” estimators can successively aggregate values per chunk to exactly derive the final value, never requiring to hold the whole sample at once. Others require the group in memory.

Estimator	Meaning	Mergeable
Mean	Arithmetic average	Yes
Standard deviation	Spread around the mean	Yes
Min	Smallest value	Yes
Max	Largest value	Yes
Sum, sum of squares	Sum of values or their squares	Yes
RMSE	Square root of the mean squared value	Yes
Counts and percentages	Valid and total counts, and their relative proportions	Yes
Median	Middle value	No
NMAD	Normalized median absolute deviation (“median-STD”)	No
Percentiles, IQR, LE90	Quantiles and ranges derived from them	No
Custom callable	User-defined calculation	No

Tip

To avoid holding a whole group in memory with a non-mergeable estimator (including a custom function), set a maximum number of samples through subsample, with subsample_per_group=True. This commonly preserves sufficient information because nearby geospatial observations are often highly correlated.

Global statistics

ds.rst.stats() or Raster.stats() with by=None
gdf.pc.stats() or PointCloud.stats() with by=None

Global statistics allow to derive any statistical estimator on any or all raster bands or point cloud columns.

By default, the mean, standard deviation, max, min and count are computed (none of which requires holding large samples in memory; see Estimators). Nodata values are skipped by default, but explicitly reported in value counts.

# Compute the default summary statistics
ds.rst.stats()

{'Min': 317.7101745605469,
 'Max': 3959.903076171875,
 'Mean': 1497.2626953125,
 'Standard deviation': 586.923583984375,
 'Valid count': 324194,
 'Total count': 333102,
 'Percentage valid points': 97.32574406638207}

Requesting one statistic returns a scalar, and requesting several statistics returns a dictionary. Use "all" to request every available estimator.

ds.rst.stats("mean")

1497.2626953125

ds.rst.stats(["mean", "median", "std", "nmad"])

{'mean': 1497.2626953125,
 'median': 1276.063720703125,
 'std': 586.923583984375,
 'nmad': 245.85694885253906}

Custom estimators receive the data array and must handle its NaNs (with Xarray accessor) or NumPy masked array (with a Raster), for example:

def count_high_elevations(data: np.ndarray) -> int:
    """Count observations above 1500 m in the selected elevation data."""

    # Count only finite observations
    values = np.ma.asarray(data).compressed()
    return int(np.count_nonzero(np.isfinite(values) & (values > 1500)))

ds.rst.stats(count_high_elevations)

94060

Statistics also accept Reducer objects, such as Mean() or Median(), including custom definitions from Customize operators and neighborhoods. A single reducer with error_structure= returns its value and an uncertainty summary; see Uncertainty for an example.

Masking

Use mask to restrict the statistics to selected locations, which accepts vector or raster/point-cloud/array mask:

# Summarize elevations inside all glacier outlines together
ds.rst.stats(["mean", "std", "valid count"], mask=glaciers)

{'mean': 1806.7179789431889, 'std': 688.4437277620498, 'valid count': 324194}

Subsampling

Use subsample to cap the sample size, which is especially useful with robust estimators, such as the median and NMAD, that are less sensitive to outliers but need to hold the whole sample in memory (see Estimators):

ds.rst.stats(["mean", "median", "std", "nmad"], subsample=100_000)

{'mean': 1498.96142578125,
 'median': 1276.377685546875,
 'std': 589.5297241210938,
 'nmad': 244.70681762695312}

Multiple values

Multi-band rasters or multi-column point clouds can be summarized at once. Use values to select the data: raster bands (starting at one), or point cloud column names; and name them for the output dictionary by using a mapping:

# Name and summarize two bands of an RGB image
rgb = gu.open_raster(gu.examples.get_path("everest_landsat_rgb"))
# Select band 1/3, naming them "red" and "blue" for the output
rgb.rst.stats(["mean", "std"], values={"red": 1, "blue": 3})

{'red': {'mean': 178.22250366210938, 'std': 75.93692779541016},
 'blue': {'mean': 182.03843688964844, 'std': 72.17345428466797}}

Grouped statistics

ds.rst.stats() or Raster.stats() with by
gdf.pc.stats() or PointCloud.stats() with by

Grouped statistics derive statistical estimators within bins or categories of one or more variables. Bins can be continuous intervals, such as elevation bands, while categories can be land-cover classes or vector zones.

Use by to name the grouping variables. Vector attributes can provide numeric grouping variables; supply bins to treat an attribute as continuous.

Masking, subsampling and selecting multiple values work as for global statistics. Set subsample_per_group=True to apply the subsample independently within each group.

The examples below summarize three RGB bands by elevation and glacier coverage. The RGB rendering, elevation raster, and glacier outlines share the same grid and location:

Show code cell content

Hide code cell content

# Create a three-band terrain rendering on the elevation grid
elevation_values = np.asarray(ds)
valid_elevation = np.isfinite(elevation_values)
lower, upper = np.nanpercentile(elevation_values, [2, 98])
normalized_elevation = np.clip((elevation_values - lower) / (upper - lower), 0, 1)
terrain_rgb_values = plt.get_cmap("terrain")(normalized_elevation)[..., :3]
terrain_rgb_values[~valid_elevation] = np.nan
terrain_rgb = gu.raster.from_array(
    np.moveaxis(terrain_rgb_values, -1, 0),
    transform=ds.rst.transform,
    crs=ds.rst.crs,
)

Show the code for plotting the figure

Hide the code for plotting the figure

fig, axes = plt.subplots(1, 2)
terrain_rgb.rst.plot(ax=axes[0])
glaciers.vct.plot(ref=ds, ax=axes[0], fc="none", ec="k", linewidth=0.6)
axes[0].set_title("Terrain RGB and glaciers")
ds.rst.plot(ax=axes[1], cmap="terrain", cbar_title="Elevation (m)")
glaciers.vct.plot(ref=ds, ax=axes[1], fc="none", ec="k", linewidth=0.6)
axes[1].set_title("Elevation and glaciers")
axes[1].set_yticklabels([])
plt.tight_layout()

Continuous bins

Use bins to define intervals of a continuous variable. For example, summarize the RGB bands within elevation bands:

# Summarize all three color bands within explicit elevation intervals
elevation_edges = [0, 1000, 2000, 3000, 4000]
elevation_stats = terrain_rgb.rst.stats(
    ("mean", "std"),
    by={"elevation": ds},
    values={"red": 1, "green": 2, "blue": 3},
    bins={"elevation": elevation_edges},
)
elevation_stats

value	red			green			blue
statistic	count	mean	std	count	mean	std	count	mean	std
elevation
[0.0, 1000.0)	21312	0.173010	0.023382	21312	0.253980	0.046764	21312	0.653980	0.046764
[1000.0, 2000.0)	249122	0.101944	0.183699	249122	0.616168	0.160946	249122	0.783736	0.174118
[2000.0, 3000.0)	39623	0.779884	0.146292	39623	0.726370	0.196383	39623	0.505689	0.085268
[3000.0, 4000.0)	14137	0.943858	0.070138	14137	0.928138	0.089777	14137	0.924769	0.093986

The output table has one row per group and one column per value–statistic combination. Each value also has a finite count.

An integer, such as bins={"elevation": 10}, creates ten equal-width bins. Numeric edges include each lower edge and the final upper edge. Pass a pandas.IntervalIndex to choose which edges are included:

# Assign boundary elevations to the interval ending at that elevation
right_closed = pd.IntervalIndex.from_breaks(elevation_edges, closed="right")
terrain_rgb.rst.stats(
    "mean",
    by={"elevation": ds},
    values={"red": 1, "green": 2, "blue": 3},
    bins={"elevation": right_closed},
)

value	red		green		blue
statistic	count	mean	count	mean	count	mean
elevation
(0, 1000]	21312	0.173010	21312	0.253980	21312	0.653980
(1000, 2000]	249122	0.101944	249122	0.616168	249122	0.783736
(2000, 3000]	39623	0.779884	39623	0.726370	39623	0.505689
(3000, 4000]	14137	0.943858	14137	0.928138	14137	0.924769

Discrete categories

Use categories to define discrete classes and their order. A boolean glacier mask, for example, separates terrain inside and outside the outlines:

# Compare RGB values inside and outside glaciers
glacier_mask = glaciers.vct.create_mask(ds)
glacier_stats = terrain_rgb.rst.stats(
    ("mean", "std"),
    by={"glacier": glacier_mask},
    values={"red": 1, "green": 2, "blue": 3},
    categories={"glacier": [False, True]},
)
glacier_stats

value	red			green			blue
statistic	count	mean	std	count	mean	std	count	mean	std
glacier
False	163011	0.050738	0.070184	163011	0.536063	0.146492	163011	0.822510	0.133217
True	161183	0.403626	0.367553	161183	0.703745	0.208606	161183	0.671384	0.201477

The same operation applies to land-cover codes or other categorical variables. Values outside the declared categories are excluded. Boolean and Pandas categorical inputs can supply their categories without an explicit declaration.

Zonal statistics

Zonal statistics assign vector feature attributes as categories. Pass (vector, "column") in by; the geometries determine each zone, and the selected column provides its category. A unique feature ID gives one category per feature, while repeated IDs combine several features into one category. Categories are inferred from the column, so no bins or categories argument is needed.

For example, the glacier inventory’s RGIId column identifies each glacier:

# Calculate separate RGB statistics for each glacier outline
glacier_zonal_stats = terrain_rgb.rst.stats(
    ("mean", "std"),
    by={"glacier": (glaciers, "RGIId")},
    values={"red": 1, "green": 2, "blue": 3},
    observed=False,
)
glacier_zonal_stats.sort_values(("red", "count"), ascending=False).head()

value	red			green			blue
statistic	count	mean	std	count	mean	std	count	mean	std
glacier
RGI60-17.15831	91912	0.410381	0.356404	91912	0.695759	0.210875	91912	0.654508	0.197729
RGI60-17.15825	22870	0.303591	0.376234	22870	0.687101	0.194452	22870	0.748345	0.205063
RGI60-17.15833	14502	0.281847	0.297474	14502	0.599426	0.254005	14502	0.680087	0.164838
RGI60-17.15808	10337	0.869822	0.133905	10337	0.833372	0.171398	10337	0.825561	0.179432
RGI60-17.15836	6492	0.754298	0.157189	6492	0.738341	0.192489	6492	0.573033	0.180187

Locations outside all zones are excluded, and each statistic uses the finite values in its category. observed=False retains categories without observations, with zero counts and undefined statistics.

The same grouping applies to point measurements:

# Summarize a sample of point elevations within each glacier
points = ds.rst.to_pointcloud(data_name="elevation", subsample=2000, random_state=42)
points.pc.stats(("mean", "std"), by={"glacier": (glaciers, "RGIId")}).head()

value	elevation
statistic	count	mean	std
glacier
RGI60-17.08503	4	1399.216583	83.531431
RGI60-17.08517	8	1382.775574	42.089448
RGI60-17.08519	5	1307.857349	62.000728
RGI60-17.08631	2	1389.033325	57.048950
RGI60-17.08642	7	1292.840297	34.394481

Raster cells are assigned by their centres, and points by their coordinates. Grouping assigns each location to one category; use non-overlapping zones for independent feature statistics. Vector and GeoDataFrame inputs are both accepted, including when the raster values are backed by Dask.

Note

by={"zone": (zones, "id")} groups by feature IDs. A bare vector, by={"inside": zones}, groups the union of all features into inside/outside categories. mask=zones instead restricts all groups to locations inside the features.

Fractional zonal statistics

For polygons cutting through raster cells, fractional=True weights each value by the exact covered fraction of its cell. The default center-based assignment described above remains useful for categorical membership; fractional coverage better describes an average over polygon area.

# Weight edge cells by the area each glacier polygon covers
fractional_glacier_stats = ds.rst.stats(
    ("mean", "sum"), by={"glacier": (glaciers, "RGIId")}, fractional=True
)
fractional_glacier_stats.head()

value	band_1
statistic	count	mean	sum
glacier
RGI60-17.08440	173.068273	1336.546053	2.313137e+05
RGI60-17.08503	723.939777	1442.353758	1.044177e+06
RGI60-17.08517	589.405006	1389.539434	8.190015e+05
RGI60-17.08519	649.770004	1354.089846	8.798470e+05
RGI60-17.08613	39.999975	1396.363853	5.585452e+04

The mean normalizes coverage weights, while the sum adds coverage-weighted values. Choose the statistic according to whether values represent means, totals or densities. In fractional summaries, count is the sum of valid covered-cell fractions, which can differ from the number of intersected cells.

This mode currently requires one vector grouper on raster values. Polygons with the same label are combined; different labels are evaluated independently, allowing a cell to contribute to several overlapping zones. overlap_backend= selects ExactExtract or Shapely, as for Fractional rasterization.

Fractional zonal statistics currently require eager raster values. Crop the raster to the area of interest before loading it; ordinary center-based zonal statistics support Dask inputs.

Several grouping variables

Several grouping variables define their joint bins or categories. For example, summarize the RGB bands jointly by elevation interval and glacier membership:

# Retain empty combinations so the complete comparison grid remains visible
joint_stats = terrain_rgb.rst.stats(
    ("mean", "std"),
    by={"elevation": ds, "glacier": glacier_mask},
    values={"red": 1, "green": 2, "blue": 3},
    bins={"elevation": elevation_edges},
    categories={"glacier": [False, True]},
    observed=False,
)
joint_stats.head()

	value	red			green			blue
	statistic	count	mean	std	count	mean	std	count	mean	std
elevation	glacier
[0.0, 1000.0)	False	16370	0.176947	0.023314	16370	0.246106	0.046628	16370	0.646106	0.046628
	True	4942	0.159968	0.018324	4942	0.280065	0.036647	4942	0.680065	0.036647
[1000.0, 2000.0)	False	146472	0.035876	0.053845	146472	0.568366	0.114671	146472	0.842661	0.124306
	True	102650	0.196218	0.250288	102650	0.684378	0.190234	102650	0.699655	0.198757
[2000.0, 3000.0)	False	169	0.706486	0.146503	169	0.625404	0.189296	169	0.445042	0.079588

The rows form a pandas.MultiIndex. Vector zones can be combined with continuous bins in the same way, by including both variables in by. observed=False retains empty combinations; the default, observed=True, omits combinations without eligible locations.

Use plot_grouped_stats() to plot a statistic and its sample counts for one or two grouping variables. min_count hides estimates based on too few observations.

Show the code for plotting the figure

Hide the code for plotting the figure

# Display the mean red value and the number of observations supporting it
fig, ax = plt.subplots(figsize=(7, 4))
fig.subplots_adjust(right=0.78)
axes = gu.stats.plot_grouped_stats(joint_stats, value="red", statistic="mean", min_count=100, ax=ax)
axes["statistic"].set(xlabel="Elevation (m)", ylabel="Inside glacier")

# Place the color scale beyond the count panels so every elevation interval stays visible
axes["colorbar"].set_position((0.84, 0.11, 0.025, 0.5))
_ = axes["colorbar"].set_ylabel("Mean red value")

To summarize only one variable, make a separate grouping call. Statistics such as medians cannot generally be combined from subgroup results.

Returning group masks

Set return_masks=True when the geographic membership of each group is also needed. The same row key retrieves its mask:

_, elevation_masks = terrain_rgb.rst.stats(
    "mean",
    by={"elevation": ds},
    values={"red": 1},
    bins={"elevation": elevation_edges},
    return_masks=True,
)
first_interval = next(iter(elevation_masks))
first_mask = elevation_masks[first_interval]
first_interval

Interval(0.0, 1000.0, closed='left')

Spatial statistics (variography)

Spatial statistics allows to describe how differences between values change with distance (i.e. spatial autocorrelation).

GeoUtils performs the pair sampling and lag-bin reduction efficiently, then relies on SciKit-GStat’s semivariance estimators to derive estimators and fit models. Models can then be passed to other GeoUtils functions or exported to covariance formats such as GSTools or GPyTorch for kriging, random field generation and other analyses.

Estimators

For mergeable estimators listed below, we also implement a specific aggregation routine to derive them without holding more than a chunk in memory.

Estimator	Meaning	Mergeable
Matheron	Half the mean squared pair difference	Yes
Cressie–Hawkins	Robust transform of the mean square-root difference	Yes
MinMax	Difference between the maximum and minimum, normalized by the mean	Yes
Dowd	Robust median-based semivariance	No
Genton	Robust order statistic of within-bin differences	No
Shannon entropy	Entropy of the pair differences	No
Percentile	Selected percentile of the pair differences	No
Custom callable	User-defined calculation	No

Variography

ds.rst.variogram() or Raster.variogram()
gdf.pc.variogram() or PointCloud.variogram()

An empirical variogram groups sampled pairs by distance and estimates their semivariance. A fitted model describes this spatial variability with a correlation range, a structured variance (partial sill) and an optional nugget.

Install geoutils[opt] for the optional geostatistical backends. GeoUtils uses SciKit-GStat estimators and models and returns a Variogram with the empirical bins and fitted parameters.

# Estimate elevation variability over distances up to 5 km using one pair sample
variogram = ds.rst.variogram(
    n_pairs=20_000,
    n_lags=12,
    max_lag=5000,
    model="spherical",
    random_state=42,
)
variogram.to_dataframe().head()

	lag	semivariance	count	semivariance_error	bin_lower_edge	bin_edge	fitted_semivariance
0	36.932452	93.599361	2246	NaN	30.000000	45.948921	221.807359
1	64.400595	310.161025	1965	NaN	45.948921	70.376778	386.752414
2	90.819636	604.938206	1226	NaN	70.376778	107.791234	545.364144
3	132.226880	1116.537030	2248	NaN	107.791234	165.096362	793.861520
4	206.291245	1798.160500	1822	NaN	165.096362	252.866657	1237.894383

Show the code for plotting the figure

Hide the code for plotting the figure

fig, ax = plt.subplots(figsize=(6, 3.5))
variogram.plot(ax=ax)
ax.set(xlabel="Distance (m)", ylabel="Elevation semivariance (m²)")
plt.tight_layout()

Distance uses coordinate units; semivariance uses squared value units. This example describes variability in terrain elevation. To estimate an error variogram, use error measurements or elevation differences on stable terrain.

Sampling and fitting

Logarithmic distance sampling represents short and long separations efficiently. n_pairs sets the sample size and n_lags sets the number of distance bins. The default uses one sample, with chunk scheduling handled by the sampling backend. See Pair-sampling for pair selection.

For an advanced estimate of sampling variability, set n_runs=3 or more to repeat sampling with shared distance bins. The result includes the standard error across those independent estimates. This describes sampling variability; it does not give confidence intervals on fitted parameters. Repetitions run sequentially, while each sample retains its backend’s chunk scheduling.

Fitting is optional. Use .fit() to fit or refit the empirical bins without repeating the sampling:

# Fit another model to the retained empirical variogram
refitted = variogram.fit("gaussian")
refitted.model.to_dict()

{'model_name': 'gaussian',
 'effective_range': 2162.60354793847,
 'partial_sill': 13109.31789197035,
 'nugget': 0.0,
 'smoothness': None,
 'shape': None,
 'active_dims': None,
 'components': []}

Use model=["spherical", "gaussian"] to fit a sum of models at different spatial scales. from_pairs() estimates an empirical variogram from an existing pair sample. The result stores bin statistics and model parameters, so the pair arrays can be discarded.

Reusing models

A Variogram stores model parameters and evaluates semivariance, covariance and correlation. Export it to a dataframe, Xarray dataset or dictionary to retain its bins, counts and fitted model:

# Keep a portable result and evaluate correlation at selected distances
dataset = variogram.to_xarray()
restored = gu.Variogram.from_dict(variogram.to_dict())
restored.correlation(np.array([0, 100, 1000]))

array([1.        , 0.96290467, 0.63653792])

Use from_model() to define a model from known parameters. Convert it to GSTools or GPyTorch to use the same spatial model in kriging, random field simulation or Gaussian process calculations:

known = gu.Variogram.from_model("gaussian", effective_range=500, partial_sill=4, nugget=0.2)
gstools_model = known.to_gstools(dim=2).model

# GPyTorch receives the nugget separately as observation noise
conversion = known.to_gpytorch()
kernel, noise = conversion.kernel, conversion.noise

Kriging

points.pc.krige() or PointCloud.krige()
ds.rst.krige() or Raster.krige()

Kriging predicts from observations and a fitted model of their spatial covariance. This differs from distance-only interpolation: the variogram determines how measurements jointly contribute to the prediction. GeoUtils supports ordinary kriging through GSTools and GPyTorch.

# Grid a small elevation sample using the fitted terrain variogram
kriging_points = ds.rst.subsample(400, random_state=42)
kriged = kriging_points.pc.krige(
    variogram, shape=(24, 24), bounds=tuple(ds.rst.bbox), max_overlap=1000, backend="gstools"
)
kriged.rst.shape

(24, 24)

The model’s longest effective range defines the default search distance. max_overlap caps that distance, in CRS units, to bound the observations and memory needed per output. The same search applies at chunk boundaries. exact=True preserves a source value when the prediction location matches it.

Supply error_structure= to account for observation errors in the fitted prediction. A variogram of the measured field and an error model for those measurements describe different quantities; a terrain variogram is not automatically an elevation error model. Use propagate() to request output uncertainty, as described in Uncertainty.

For simulations from an error model instead of predictions from measurements, see Uncertainty.

Chunked strategies

In construction