Extracting SIGRID 1 km plots for a whole country, tagged with province (ADM1) and district (ADM2) names

I put together a small Python script that automates a common request: take the SIGRID 1000 m grid and pull out every plot that falls inside a given country, with the province and district name attached to each plot. Optionally it also splits the result into one file per SIGRID density (1x1 km, 2x2 km, ... up to 100x100 km).

The script is hosted next to the grids so you can grab it directly:

https://www.openforis.org/fileadmin/SIGRID_1000m_grids/sigrid_country.py

You only need to give it a 2-letter ISO country code, for example bj for Benin or ye for Yemen. Everything else - downloading the right grid tiles, downloading the country boundaries, and doing the spatial join - is automatic.

What you get

For a country you get a single CSV, <iso2>_plots_with_districts.csv, with the original SIGRID columns plus two new ones, ADM1NM (province/department) and ADM2NM (district), and containing only the plots that fall inside the country:

CE_ID, yCoordinate, xCoordinate, ADM1NM, ADM2NM, grid_1, grid_2, grid_3, ..., grid_100

With the --split-by-density option you additionally get one CSV per density grid, named by country and resolution - Benin_1x1km.csv, Benin_2x2km.csv, Benin_3x3km.csv, ... each holding only the plots flagged for that grid.

Prerequisites

• Python 3.8+

• Two libraries:

pip install pandas geopandas

geopandas pulls in shapely, pyproj and fiona, which handle the geometry and the point-in-polygon test.

Usage

Download sigrid_country.py, then run:

Benin - combined country file only

python sigrid_country.py bj

Benin - also split into one file per density grid

python sigrid_country.py bj --split-by-density

Choose an output directory

python sigrid_country.py ke --out ./results --split-by-density

The first run downloads the grid tiles it needs and caches them under ./sigrid_cache/, so re-runs, or a second country in the same tiles, are fast and do not re-download.

How it works

The script chains together the two data sources that already exist on the Open Foris / EarthMap side:

1. Country boundaries, UN 2023, ADM2. It downloads the country's administrative boundaries as GeoJSON from EarthMap's public store, using the ISO2 code:

https://firebasestorage.googleapis.com/v0/b/production-earthmap.appspot.com/o/boundaries%2FUN2023%2F{iso2}_adm2.geojson?alt=media

These features carry the ADM1NM and ADM2NM attributes, which is where the province and district names come from.

1. Which SIGRID tiles to fetch. The SIGRID 1000 m grids are published as 10 x 10 degree tiles, for example SIGRID_x_0_10_y_0_10_1000m_1_subgrid.csv.zip. Each ZIP holds three subgrid_0/1/2 CSVs. The script reads the bounding box of the country boundaries and works out exactly which tiles overlap it, then downloads only those. Benin needs 2 tiles; Yemen needs 2; larger countries need more.

2. Filter + spatial join. For each subgrid CSV it first drops everything outside the country's bounding box, which is fast, then runs a point-in-polygon join using geopandas.sjoin, predicate="within", CRS EPSG:4326, against the ADM2 polygons. Plots that do not fall inside any district are dropped; the rest get their ADM1NM / ADM2NM values. Points sitting exactly on a shared border are assigned to a single district to avoid duplicates.

3. Output. Matched plots from all tiles are streamed into one combined CSV. If --split-by-density is set, that CSV is then split by the grid_N flag columns.

About the density grids

Each SIGRID plot carries boolean flags grid_1 ... grid_100. grid_N = true means the plot is a sampling point of the N x N km grid. grid_1 is true for every plot, the full 1 km grid; coarser grids are progressively sparser. The densities actually present are:

1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 20, 25, 30, 50, 100 km

As a sanity check, the number of plots roughly scales with 1/N^2. For example, the 2x2 km grid has about a quarter of the 1x1 km plots, the 10x10 km grid about 1/100, and so on.

Worked examples

Benin (bj): 115,259 plots across all 12 departments and 77 districts. Two tiles downloaded, around 30 s of processing.

Yemen (ye): 453,565 plots across 22 governorates. With --split-by-density this produced 17 files, from Yemen_1x1km.csv, with 453,565 plots, down to Yemen_100x100km.csv, with 46 plots.

Notes & troubleshooting

• ISO2 code must be lowercase and match the boundary file name, for example bj, ye, ke. If the script reports it cannot download the GeoJSON, double-check the code.

• Ocean / empty tiles. Some 10-degree tiles simply do not exist on the server, for example open ocean. The script skips a missing tile rather than failing.

• Coordinate system. SIGRID coordinates are plain lon/lat, WGS84 / EPSG:4326. The boundaries are reprojected to the same CRS before the join, so no manual projection is needed.

• Accented names. ADM names come straight from the UN 2023 boundaries. In the output CSV files, these names keep their original diacritics, for example Oueme or Hadramawt in simplified ASCII form here. Files are written UTF-8.

• Large countries. Each tile is around 11 MB zipped but around 85 MB per subgrid CSV uncompressed, and a big country can span many tiles. Make sure you have disk space for the cache, and expect longer run times. The script reads one subgrid at a time and streams output to keep memory use modest.

• A small CSV quirk in the source files, two extra trailing commas per row, is handled internally. You do not need to pre-clean anything.