Generate area-weighted statistics from raster coverage datasets for APEX boundary regions for use in the APEX Geospatial Explorer.
apex-build-statistics intersects polygon boundaries with raster datasets (such as land cover or soil products) and computes aggregated statistics for each region. Results are exported as FlatGeobuf (FGB) datasets with metadata for efficient downstream analysis and visualisation.
The APEX Geospatial Explorer includes functionality to display dataset statistics based on polygon areas, some of which represent hierarchical administrative regions.
This tool prepares statistics files (such as the area of each land-cover class or soil property range per region) by analysing raster datasets using vector boundaries prepared by the apex-convert-boundaries repository.
The APEX workflow consists of two complementary repositories:
-
apex-convert-boundaries
Prepares hierarchical boundary datasets. -
apex-build-statistics
Generates raster coverage statistics using those prepared boundaries.
Together they form the complete geospatial processing pipeline.
flowchart LR
A[Source Boundary Datasets<br>GAUL / NUTS / Natura2000]
--> B[apex-convert-boundaries<br>Standardise and enrich boundaries]
B --> C[APEX Boundary Dataset<br>FlatGeobuf + hierarchical metadata]
C --> D[apex-build-statistics<br>Raster summarisation engine]
E[Raster Coverage Datasets<br>WorldCover / WorldSoils / others]
--> D
D --> F[Coverage Statistics Outputs<br>Area-weighted summaries<br>per region and coverage]
apex-build-statistics consumes the prepared boundary datasets and intersects them with raster coverage datasets.
Examples of supported coverages include:
- ESA WorldCover
- global soil property datasets
- other environmental raster layers
The tool:
- identifies raster tiles intersecting each region
- processes raster blocks efficiently
- classifies pixel values into bins
- computes area-weighted statistics
- Raster–vector intersection using Rasterio, GeoPandas, and Shapely
- Area-weighted raster statistics per polygon
- Support for multiple raster coverages
- Efficient block-based raster processing
- Automatic tiling of large raster datasets
- CLI interface with progress reporting
- Hierarchical outputs by APEX level
- FlatGeobuf outputs optimised for large datasets
- Python 3.10+
- GDAL / OGR
- Rasterio
- GeoPandas
- NumPy
git clone https://github.com/sparkgeo/apex-build-statistics.git
cd apex-build-statistics
pip install -e .For development dependencies:
pip install -e ".[dev]"build-statistics coveragesbuild-statistics process boundaries.fgb worldcoverSpecify an output directory:
build-statistics process boundaries.fgb worldcover -o output/Boundary datasets must contain APEX attributes used by the processing pipeline. These can be prepared using the apex-convert-boundaries repository.
| Field | Description |
|---|---|
apex:id |
Unique region identifier |
apex:name |
Region name |
apex:level |
Hierarchy level |
apex:children |
Child region references |
apex:total_area |
Polygon area |
Supported formats include any GDAL-compatible vector format:
- FlatGeobuf
- GeoPackage
- Shapefile
- GeoJSON
Results are written as FlatGeobuf datasets grouped by hierarchy level.
Example outputs:
stats.worldcover.boundaries.epsg4326.level00.fgb
stats.worldcover.boundaries.epsg4326.level01.fgb
stats.worldcover.boundaries.epsg4326.level02.fgb
Each output dataset contains:
- Input geometry
- Original APEX attributes
- Statistic columns for each raster bin
- Metadata describing the coverage
| Coverage | Description |
|---|---|
worldcover |
ESA WorldCover land cover dataset |
worldsoils |
Global soil property datasets |
Additional coverages can be added easily.
src/build_statistics/
│
├── cli.py
├── process.py
├── models.py
├── utils.py
├── spinner.py
├── exceptions.py
└── coverages/
├── worldcover.py
└── worldsoils.py
Create a coverage class in:
src/build_statistics/coverages/
Example:
from build_statistics.models import BaseCoverageClass
class ExampleCoverage(BaseCoverageClass):
registry_key = "example"
urls = [...]
units = "example_units"
bin_edges = [0, 10, 20, 30]
bin_labels = ["0-10", "10-20", "20-30"]
bounds = (-180, -90, 180, 90)
def get_tile_keys(self, df):
...Run formatting and linting:
pre-commit run --all-filesInstall development dependencies:
pip install -e ".[dev]"Contributions are welcome.
Typical workflow:
- Fork the repository
- Create a feature branch
- Add tests or documentation
- Submit a pull request
This project is licensed under the terms of the license included in the repository.
- APEX geospatial boundary datasets
- ESA WorldCover
- Open-source geospatial Python ecosystem (Rasterio, GeoPandas, GDAL)