GHSL — settlement classes (GHS-SMOD)¶

GHS-SMOD is the Degree of Urbanisation settlement model — a categorical grid (urban centre / clusters / rural / water). It has no tiled resolution, so this downloads the whole-globe 1 km file (~100 MB) once and crops it to a Portugal-sized AOI. The output carries a class colour table and a .legend.json sidecar, and reprojection uses nearest-neighbour so class codes are never blended.

The code cells are marked NBVAL_SKIP so the docs test suite does not re-download the global file; the saved outputs below are real.

Setup¶

Consolidate the imports up front. earthlens provides the unified EarthLens entry point and the GHSL Catalog (class codes, labels, and colours); pyramids reads the downloaded GeoTIFF; matplotlib renders the categorical map.

In [ ]:

# NBVAL_SKIP
import json
import tempfile

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import BoundaryNorm, ListedColormap
from pyramids.dataset import Dataset

from earthlens import EarthLens
from earthlens.ghsl import Catalog

# NBVAL_SKIP import json import tempfile import matplotlib.pyplot as plt import numpy as np from matplotlib.colors import BoundaryNorm, ListedColormap from pyramids.dataset import Dataset from earthlens import EarthLens from earthlens.ghsl import Catalog

1 · Download and crop GHS-SMOD¶

Build the request first — source, the GHS_SMOD product, the 2020 epoch, and a Portugal-sized AOI — into a temporary output directory. Keeping construction on its own line makes the request easy to read and re-run.

In [ ]:

# NBVAL_SKIP
out = tempfile.mkdtemp()
app = EarthLens(
    data_source='ghsl',
    variables=['GHS_SMOD'],
    start='2020-01-01',
    end='2020-12-31',
    aoi=[-9.6, 36.9, -6.0, 42.2],
    path=out,
)

# NBVAL_SKIP out = tempfile.mkdtemp() app = EarthLens( data_source='ghsl', variables=['GHS_SMOD'], start='2020-01-01', end='2020-12-31', aoi=[-9.6, 36.9, -6.0, 42.2], path=out, )

download() fetches the global 1 km file, crops it to the AOI with nearest-neighbour reprojection, and returns the list of written paths.

In [ ]:

# NBVAL_SKIP
paths = app.download(progress_bar=False)
paths

# NBVAL_SKIP paths = app.download(progress_bar=False) paths

2 · The class legend¶

The categorical raster ships a class-code → label legend sidecar next to the GeoTIFF. Reading it shows the settlement classes encoded in the pixel values.

In [ ]:

# NBVAL_SKIP
legend = json.loads(paths[0].with_suffix('.legend.json').read_text())
legend

# NBVAL_SKIP legend = json.loads(paths[0].with_suffix('.legend.json').read_text()) legend

3 · Plot the settlement classes¶

Look up the product metadata from the GHSL Catalog and read the cropped raster into a 2-D array, dropping the leading band axis if present.

In [ ]:

# NBVAL_SKIP
prod = Catalog().get('GHS_SMOD')
arr = np.asarray(Dataset.read_file(str(paths[0])).read_array(), dtype='float32')
arr = arr[0] if arr.ndim == 3 else arr

# NBVAL_SKIP prod = Catalog().get('GHS_SMOD') arr = np.asarray(Dataset.read_file(str(paths[0])).read_array(), dtype='float32') arr = arr[0] if arr.ndim == 3 else arr

Build a discrete colormap and a BoundaryNorm from the catalog's class codes so each settlement class maps to its own GHSL legend colour, then render the map with a labelled colour bar.

In [ ]:

# NBVAL_SKIP
codes = sorted(prod.legend)
cmap = ListedColormap([prod.colors[c] for c in codes])
norm = BoundaryNorm([c - 0.5 for c in codes] + [codes[-1] + 0.5], cmap.N)

plt.figure(figsize=(6, 6))
im = plt.imshow(arr, cmap=cmap, norm=norm)
cbar = plt.colorbar(im, ticks=codes, shrink=0.8)
cbar.ax.set_yticklabels([prod.legend[c] for c in codes])
plt.title('GHS-SMOD 2020 (1 km) — Portugal')
plt.axis('off')
plt.show()

# NBVAL_SKIP codes = sorted(prod.legend) cmap = ListedColormap([prod.colors[c] for c in codes]) norm = BoundaryNorm([c - 0.5 for c in codes] + [codes[-1] + 0.5], cmap.N) plt.figure(figsize=(6, 6)) im = plt.imshow(arr, cmap=cmap, norm=norm) cbar = plt.colorbar(im, ticks=codes, shrink=0.8) cbar.ax.set_yticklabels([prod.legend[c] for c in codes]) plt.title('GHS-SMOD 2020 (1 km) — Portugal') plt.axis('off') plt.show()