BIOMASS MAAP
It is recommended to save all of your files and folders within my-private-bucket, since everything that is outside of this folder is permanently deleted after seven days from the last usage.
Search and Download Data Through Catalog API
The following notebook provides some usage example of data visualization and download from the STAC catalog. If you would like to have a deeper understanding about the STAC catalog, collections and products, you can refer to the official ESA STAC API guide. The following examples are written in Python.
Catalog data visualization
Here you will understand how to perform search on the STAC catalog through the STAC API interface. The catalog can also be visualized here via browser.
import requests
import pandas as pd
import json
from typing import Any, Dict
from pystac_client import Client
from IPython.display import Markdown as md
URL_LANDING_PAGE = "https://catalog.maap.eo.esa.int/catalogue/"
api = Client.open(URL_LANDING_PAGE) # Connection to the catalogCollection Search
There are different ways you can search for a specific collection within the catalog. Some usage examples are reported below.
Search by free text
from pystac_client import Client, ConformanceClasses
import urllib
value = 'Biomass' # Search for the word "biomass" within the products
params = { 'q': value }
URL = f'{URL_LANDING_PAGE}collections?{urllib.parse.urlencode(params)}'
response = requests.get(URL)
data = json.loads(response.text)
df = pd.json_normalize(data, record_path=['collections'])
df[['id']] # Display the id of the collection
Search by title
value = 'Biomass Simulated data'
params = { 'filter': "title='" + value + "'"}
URL = f'{URL_LANDING_PAGE}collections?{urllib.parse.urlencode(params)}'
response = requests.get(URL)
data = json.loads(response.text)
df = pd.json_normalize(data, record_path=['collections'])
df[['id', 'title']]
Search by platform
URL = URL_LANDING_PAGE + "collections"+ "?filter=platform='BIOMASS'"
response = requests.get(URL)
data = json.loads(response.text)
df = pd.json_normalize(data, record_path=['collections'])
df[['title', 'summaries.platform']]
Search by organization
URL = URL_LANDING_PAGE + "collections"+ "?filter=organisationName='ESA/ESRIN'"
df = pd.json_normalize(data, record_path=['collections'])
df[['title', 'providers']]
Search by bounding box
URL = URL_LANDING_PAGE + "collections"+ "?bbox=14.90,37.700,14.99,37.780" # Longitude, Latitude
response = requests.get(URL)
data = json.loads(response.text)
df = pd.json_normalize(data, record_path=['collections'])
df[['id', 'extent.spatial.bbox']]
Search by temporal extent
URL = URL_LANDING_PAGE + "collections"+ "?datetime=" + '2020-01-01T00:00:00.000Z/2024-12-31T23:59:59.999Z'
response = requests.get(URL)
data = json.loads(response.text)
df = pd.json_normalize(data, record_path=['collections'])
df[['id', 'extent.temporal.interval']]
Granule Search
Search by geometry
from pystac_client import Client
URL_LANDING_PAGE = "https://catalog.maap.eo.esa.int/catalogue/"
api = Client.open(URL_LANDING_PAGE)
aoi_as_dict: Dict[str, Any] = {
"type": "Polygon",
"coordinates": [
[
[
112.82476,
-2.66676
],
[
112.291824,
-2.778783
],
[
112.409676,
-3.33663
],
[
112.94324,
-3.224744
],
[
112.82476,
-2.66676
]
]
]
}
results = api.search(
method = 'GET',
max_items = 5, # Maximum number of granules to take
collections = 'BiomassSimulated', # Search for granules belonging to collection ID (e.g. BiomassLevel1cIOC, BiomassAux...)
intersects = aoi_as_dict,
datetime = ['2015-01-01T00:00:00Z', '2020-01-02T00:00:00Z'] # Search for granules in date range
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')Search by bounding box
results = api.search(
method = 'GET',
max_items = 10, # Maximum number of granules to take
collections = 'BiomassSimulated', # Search for granules belonging to collection ID (e.g. BiomassLevel1cIOC, BiomassAux...)
bbox = [112.291824, -3.33663, 112.94324, -2.66676],
# datetime = ['2015-01-01T00:00:00Z', '2020-01-02T00:00:00Z'] # Search for granules in date range
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')Search by temporal extent
results = api.search(
method = 'GET',
max_items = 50, # Maximum number of granules to take
collections = 'BiomassSimulated', # Search for granules belonging to collection ID (e.g. BiomassLevel1cIOC, BiomassAux...)
datetime = ['2017-01-01T00:00:00Z', '2017-12-02T00:00:00Z'] # Search for granules in date range
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')Search by identifier
product_id = ['BIO_S1_DGM__1S_20170101T222105_20170101T222114_I_G01_M01_C01_T010_F308_01_D6K4O3',
'BIO_S1_DGM__1S_20170101T222150_20170101T222211_I_G01_M01_C01_T011_F001_01_D5U6F7'] # Insert one or more product id to search
results = api.search(
method = 'GET',
collections = 'BiomassSimulated',
ids = product_id
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')Search with filter
results = api.search(
method = 'GET',
max_items = 10, # Maximum number of granules to take
collections = 'BiomassSimulated',
# bbox = [112.291824, -3.33663, 112.94324, -2.66676],
# datetime = ['2015-01-01T00:00:00Z', '2020-01-02T00:00:00Z'] # Search for granules in date range
filter="productType='S1_DGM__1S' and instrument='P-SAR'" # Many other filters can be applied
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')Catalog Data Download
Now that you have seen how to search for collections and products, let's see how you can download data from catalog directly via code.
Token to access the catalog
To download data from the STAC catalog directly via code, you need to generate an access token to pass into the client. This token can be generated by following this page and it will last for 90 days. Then, you will be asked to log in with your credentials to generate a new token.
It is highly recommended to create a .txt file (e.g. credentials.txt) containing all the authentication information that you need. In Coding, you can create a new .txt file by clicking on "Text File" from the Launcher:
The credentials.txt file should be located in the same directory as the script.
Now you can populate your file with your keys and tokens in this format:
CLIENT_ID=offline-token
CLIENT_SECRET=p1eL7uonXs6MDxtGbgKdPVRAmnGxHpVE
OFFLINE_TOKEN=your_esamaap_longlasting_token_hereYou only have to change OFFLINE_TOKEN, pasting your long lasting token generated before. The other two entries should be kept as they are in the above example.
Be careful to not use quotes around the values and ensure there are no trailing spaces or extra characters.
Example Code
Below you can find a sample of code that shows you how to use your token.
Let's consider the case in which (for example) you want to download all the products that satisfy your search requirements:
- Belong to the collection "BiomassLevel1aIOC"
- Acquisition date between 2025-01-01 and 2025-12-31 (YYYY-MM-dd)
- With bounding box [-57.099045, -20.350586, -56.448875, -19.722197]
import os
import requests
from pathlib import Path
from pystac_client import Client
# --- Path to credentials.txt ---
CREDENTIALS_FILE = Path('path/to/credentials.txt').resolve().parent / "credentials.txt" # Insert the .txt path
def load_credentials(file_path=CREDENTIALS_FILE):
"""Read key-value pairs from a credentials file into a dictionary."""
creds = {}
if not file_path.exists():
raise FileNotFoundError(f"Credentials file not found: {file_path}")
with open(file_path, "r") as f:
for line in f:
line = line.strip()
if not line or line.startswith("#"):
continue
if "=" not in line:
continue
key, value = line.split("=", 1)
creds[key.strip()] = value.strip()
return creds
# --- ESA MAAP API ---
def get_token():
"""Use OFFLINE_TOKEN to fetch a short-lived access token."""
creds = load_credentials()
OFFLINE_TOKEN = creds.get("OFFLINE_TOKEN")
CLIENT_ID = creds.get("CLIENT_ID")
CLIENT_SECRET = creds.get("CLIENT_SECRET")
print(CLIENT_SECRET)
if not all([OFFLINE_TOKEN, CLIENT_ID, CLIENT_SECRET]):
raise ValueError("Missing OFFLINE_TOKEN, CLIENT_ID, or CLIENT_SECRET in credentials file")
url = "https://iam.maap.eo.esa.int/realms/esa-maap/protocol/openid-connect/token"
data = {
"client_id": CLIENT_ID,
"client_secret": CLIENT_SECRET,
"grant_type": "refresh_token",
"refresh_token": OFFLINE_TOKEN,
"scope": "offline_access openid"
}
response = requests.post(url, data=data)
response.raise_for_status()
response_json = response.json()
access_token = response_json.get('access_token')
if not access_token:
raise RuntimeError("Failed to retrieve access token from IAM response")
return access_token
if __name__ == "__main__":
URL_LANDING_PAGE = "https://catalog.maap.eo.esa.int/catalogue/"
api = Client.open(URL_LANDING_PAGE) # Connection to the catalog
results = api.search(
method = 'GET',
max_items = 20, # Maximum number of granules to take
collections = 'BiomassLevel1aIOC',
bbox = [-57.099045, -20.350586, -56.448875, -19.722197],
datetime = ['2025-01-01T00:00:00Z', '2025-12-31T00:00:00Z'] # Search for granules in date range
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')
access_token = get_token()
data = results.item_collection_as_dict()
n_products = 2 # Number of products to download
for n in range(0,n_products):
file_url = data['features'][n]['assets']['product']['href']
file_path = "./my_products/" + data['features'][n]['assets']['product']['file:local_path'] # Change ./my_products/ with your desired path to download data
if access_token:
print("Access token verified!")
else:
print("Failed to retrieve access token.")
exit(2)
try:
headers = {"Authorization": f"Bearer {access_token}"}
response = requests.get(file_url, headers=headers, stream=True)
response.raise_for_status() # Raise an exception for bad status codes
with open(file_path, "wb") as f:
for chunk in response.iter_content(chunk_size=8192):
f.write(chunk)
print(f"File downloaded successfully to {file_path}")
print('')
except requests.exceptions.RequestException as e:
print(f"Error downloading file: {e}")
print('')Science Examples
In this chapter you will find a variety of real-world examples that may help you create a custom Jupyter Notebook for your research. These notebooks have dataset specific analysis and visualization examples, and will also provide examples of analysis, data access, tabulation and map-visualization.
AGB reference maps in R using airborne laser scanning (ALS)
This workflow is a demonstration of one possible approach for creating aboveground biomass reference maps in R using airborne laser scanning (ALS) and field plot inventory data collected over plots located across the Canal Zone in central Panama. This version tests the approach over two flightlines of ALS collected at Panama Pacifico.
It is highly recommended to run this tutorial within the "Rtools add-on" environment, which already includes useful R packages, such as lidR. Running this tutorial outside of this environment may lead to errors.
Using local relationships between forest inventory data and ALS, this workflow demonstrates ALS processing and generation of some illustratory plots for ALS, canopy products, and plot shapefiles. Users of this workflow will:
- Prepare the workspace.
- Access the ALS data and read it into a workspace.
- Process the ALS data to produce terrain and canopy products and ALS metrics.
- Access the field plot data and associated geospatial data.
- Generate figures illustrating the datasets.
Prepare the workspace
- Load the necessary packages in R:
install.packages('httr2')
install.packages('randomForest')
install.packages('ggpubr')
install.packages('RStoolbox')
install.packages('RCSF')
install.packages('lasR', repos = c('https://r-lidar.r-universe.dev',
'https://cloud.r-project.org'),
dependencies = TRUE)
install.packages('remotes')
remotes::install_github('njtierney/broomstick')
install.packages('tidyterra')
# load in the required R packages
suppressMessages({
library(broom)
library(broomstick)
library(ggpubr)
library(sf)
library(future)
library(terra)
library(lidR)
library(RCSF)
library(RStoolbox)
library(httr2)
library(randomForest)
library(tidyterra)
library(tidyverse)
library(lasR)
})Access the ALS data and read it into the workspace as a LAScatalog
# make the subdirectories if they don't exist
root_dir = "./data/ALS/BCI_Gigante"
dir_list = c(
paste0(root_dir, "/output/_1_classified_points"),
paste0(root_dir, "/output/_2_DTM_raster"),
paste0(root_dir, "/output/_3_normalized_points"),
paste0(root_dir, "/output/_4_CHM_raster"),
paste0(root_dir, "/output/_5_metrics_raster")
)
for (i in seq_along(dir_list)) {
dir_path = dir_list[i]
ifelse(!dir.exists(file.path(dir_path)),
dir.create(file.path(dir_path), recursive = TRUE),
FALSE)
}Process the ALS data to produce terrain and canopy products
This workflow produces five levels of products, including two transformations of the point cloud and three categories of raster data:
Classified_points
In the first step, read in the ALS flightlines from their source; performing some filtering on read. Only read one point per 2x2 cm voxel, drop points below the lowest elevation expected at the site, only read points with scan angles +/- 15 degrees from nadir, and drop likely noise points with very high numbers of returns or return numbers. Next, classify isolated points occurring above or below the canopy as noise. These points will affect canopy metrics downstream, and should be removed. They can represent real returns from birds or aerosols or can be artifacts, especially near water features. Next, classify ground returns using a cloth simulation filter algorithm. Ground classification is necessary to model terrain because point cloud Z values represent absolute elevation rather than canopy height. To quantify canopy height, it is necessary to subtract the elevation of terrain.
DTM_raster
From classified ground returns, produce a raster of terrain elevations by producing a triangualted irregular network (TIN) from the lowest value within each pixel. Also perform some smoothing to produce a realistic model for canopy height normalization.
Normalized_points
Subtract the terrain value from the Z values in the point cloud to produce a new point cloud where Z represents normalized canopy height.
CHM_raster
Subtract the terrain value from the Z values in the point cloud to produce a new point cloud where Z represents normalized canopy height.
Metrics_raster
Generate canopy metrics as predictors for modeling of aboveground biomass at a pixel resolution of 50 m:
- z_p99: the height below which 99% of returns occur,
- z_mean: the mean height of returns,
- z_cv: the coefficient of variation of return heights,
- z_above5: canopy cover above 5 m,
- a_mean: mean absolute scan angle,
- z_count: the number of returns.
Then generate derived metrics:
- slope: the mean slop derived from the DTM,
- elevation: the mean elevation from the DTM,
- volume: canopy volume as the sum of pixel heights from the smoothed CHM,
- LAI: leaf area index above 5 m; derived from canopy cover adjusted for scan angle
# _1_classified_points
path_als = paste0(root_dir, "/input")
if (length(list.files(dir_list[1])) == 0) {
# only read one point per 2x2 cm voxel,
# drop points below the lowest elevation expected at the site
# only read points with scan angles +/- 15 degrees from nadir,
# drop likely noise points with very high numbers of returns or return numbers
pipeline = reader(filter = "-thin_with_voxel .02
-keep_z_above 20
-drop_abs_scan_angle_above 15
-keep_return 1 2 3 4 5 6 7
-drop_number_of_returns 15
-drop_number_of_returns 14
-drop_number_of_returns 13
-drop_number_of_returns 12
-drop_number_of_returns 11
-drop_number_of_returns 10
-drop_number_of_returns 9
-drop_number_of_returns 8
-drop_withheld"
) +
# classify points as noise if they have fewer than 200 points in their surrounding 10x10 m voxels,
# for big clouds of above- or belowground noise
classify_with_ivf(res = 10, n = 200L, class = 18L) +
# classify points as noise if they have fewer than 30 points in their surrounding 5x5 m voxels
classify_with_ivf(res = 5, n = 30L, class = 18L) +
# classify ground points with a cloth simulation filter algorithm
classify_with_csf(slope_smooth = FALSE,
class_threshold = 0.05,
cloth_resolution = 0.3,
rigidness = 3L,
iterations = 500L,
time_step = 0.6,
class = 2L,
# only let the algorithm use last returns that aren't noise
filter = "Classification != 18") +
# write the point clouds with ground and noise points classified
write_las(paste0(root_dir, "/output/_1_classified_points/*_ground.las"))
exec(pipeline,
on = path_als,
ncores = 16,
with = list(chunk = 600),
progress = FALSE)
}
# _2_DTM_raster
if (length(list.files(dir_list[2])) == 0) {
pipeline = reader_las(filter = "-keep_class 2") +
# create a point cloud of only the lowest ground point in each 50 cm grid cell
filter_with_grid(1, operator = "min", filter = keep_ground()) +
# generate a triangulation of these points
(tri = triangulate(max_edge = 120, filter = "",
#use_low = TRUE,
use_attribute = "Z")) +
# rasterize the triangulation at a 1-m resolution
(dtm = rasterize(1, tri, ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM.tif"))) +
# apply a laplachian filter to replace spikes and pits with median values
(spikefree1 = pit_fill(dtm,
lap_size = 4L,
thr_lap = 0.05,
thr_spk = -0.05,
med_size = 4L,
dil_radius = 1L,
ofile = "")) +
# perform this filter again at a coarser scale to capture larger spikes
pit_fill(spikefree1,
lap_size = 9L,
thr_lap = 0.2,
thr_spk = -0.08,
med_size = 9L,
dil_radius = 1L,
ofile = paste0(root_dir, "/output/_2_DTM_raster/DTM_spikefree.tif"))
lasR::exec(pipeline,
on = paste0(root_dir, "/output/_1_classified_points"),
ncores = 16,
progress = FALSE)
DTM = rast(paste0(root_dir, "/output/_2_DTM_raster/DTM_spikefree.tif"))
DTM[!is.finite(DTM)] = NA
#smooth the DTM a little bit
DTM_smooth = DTM %>%
terra::crop(ext(readLAScatalog(path_als))) %>%
terra::focal(w=5, fun="median", na.policy="only", na.rm=TRUE) |>
terra::focal(w = terra::focalMat(DTM, .5, "Gauss"), na.policy="all", na.rm=FALSE,
filename = paste0(root_dir, "/output/_2_DTM_raster/DTM_smooth.tif"), overwrite = TRUE)
}
# _3_normalize_points
if (length(list.files(dir_list[3])) == 0) {
pipeline =
lasR::reader_las(filter = "-drop_class 18") +
callback(fun = function(data) {
data |> dplyr::mutate(ScanAngle = abs(ScanAngle / 166.666667))
},
expose = "a") +
(dtm_smooth = load_raster(paste0(root_dir, "/output/_2_DTM_raster/DTM_smooth.tif"))) +
transform_with(dtm_smooth, "-") +
write_las(paste0(root_dir, "/output/_3_normalized_points/*_norm.las"))
exec(pipeline,
on = paste0(root_dir, "/output/_1_classified_points"),
ncores = 20,
with = list(buffer = 10),
progress = FALSE)
}
# _4_CHM_raster
if (length(list.files(dir_list[4])) == 0) {
pipeline =
lasR::reader_las(filter = "-drop_class 18") +
chm(res = 1, tin = FALSE, ofile = paste0(root_dir, "/output/_4_CHM_raster/CHM.tif")) +
delete_points(keep_z_above(0)) +
rasterize(50, operators = c("z_p99", "z_mean", "z_cv", "z_above5", "a_mean", "z_count"),
ofile = paste0(root_dir, "/output/_5_metrics_raster/metrics_50m.tif"))
exec(pipeline,
on = paste0(root_dir, "/output/_3_normalized_points"),
ncores = 12,
progress = FALSE)
# make a smooth, clean CHM for visualization and volume extraction
CHM = rast(paste0(root_dir, "/output/_4_CHM_raster/CHM.tif"))
DTM_smooth = rast(paste0(root_dir, "/output/_2_DTM_raster/DTM_smooth.tif"))
CHM_smooth = CHM %>%
# any values below 0, set to 0
terra::clamp(lower = 0, values = TRUE) %>%
# fill in NA values with their neighbors' min value
terra::focal(w=3, fun="min", na.policy="only", na.rm=TRUE, expand = FALSE) %>%
ifel(!is.na(crop(DTM, .)) & is.na(.), 0, .) %>%
# apply a median filter only to cells where there are big changes in Z value
ifel(terra::terrain(., "TRI") > 5,
terra::focal(., w = 3, fun = "median", na.policy = "omit", na.rm = TRUE, expand = FALSE),
., filename = paste0(root_dir, "/output/_4_CHM_raster/CHM_smooth.tif"), overwrite = TRUE)
}
# _5_metrics_raster
# the metrics produced in the previous step, CHM, and DTM
metrics_50m = rast(paste0(root_dir, "/output/_5_metrics_raster/metrics_50m.tif")) |>
# number of points / (50 * 50 m) to get point density
mutate(density = z_count / 2500) |>
# at least 12 pts / m^2 and mean height above 2 m (forest)
filter(density > 10 & z_mean > 2)
CHM_smooth = rast(paste0(root_dir, "/output/_4_CHM_raster/CHM_smooth.tif")) |> crop(metrics_50m)
DTM_smooth = rast(paste0(root_dir, "/output/_2_DTM_raster/DTM_smooth.tif")) |> crop(metrics_50m)
#---------# slope
# produce slopes from 1 m DTM, then average over 50 m pixels
slope = terra::terrain(DTM_smooth, v = "slope") %>%
terra::resample(metrics_50m, method = "average")
#---------# roughness
# produce terrain roughness from 1 m DTM, then average over 50 m pixels
elevation = DTM_smooth %>%
terra::resample(metrics_50m, method = "average")
#---------# volume
# sum canopy height volume over each 50 m pixel
volume = terra::resample(CHM_smooth, metrics_50m, method = "sum")
#---------# LAI
# convert fractional cover to LAI, according to scan angle
# LAI cannot be calculated with a fractional cover of 100%
# for pixels with no points below the threshold, add 1 artificial point below it
cover_new = metrics_50m$z_count / (metrics_50m$z_count + 1)
cover = metrics_50m$z_above5 %>%
ifel(. == 1, cover_new, .)
LAI = (-1/.5) * log(1 - cover) * cos(metrics_50m$a_mean /180*pi)
# combine the metrics of interest as predictors
predictors = c(metrics_50m$z_p99, volume, LAI, metrics_50m$z_cv, slope, elevation)
names(predictors) = c("Z_p99", "Volume", "LAI_5m", "Z_CV", "Slope", "Elevation")
names(predictors)
writeRaster(predictors, filename = paste0(root_dir, "/output/_5_metrics_raster/predictors_50m.tif"),
overwrite = TRUE)
Access the field plot data and associated geospatial data
Read in plot shapefiles as sf (simple feature) object in the sf package.
small_plots = read_sf('data/GIS/CTFS_Plots_Polygons.shp') |>
mutate(NAME = DESC_,
ORIGIN_X = Origin_X,
ORIGIN_Y = Origin_Y) |>
dplyr::select(NAME, AREA_HA, ORIGIN_X, ORIGIN_Y, geometry) |>
st_set_crs(32617)
plots = small_plots |>
filter(NAME %in% c("PanamaPacifico"))
# this buffered version will be used for cropping raster data for plotting
plots_buf_200m = plots |>
st_centroid() |>
st_buffer(200) |>
suppressWarnings()Generate figures illustrating the datasets
# define some variables for plotting point cloud transects in ggplot
theme_transect = theme(
panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_text(margin = ggplot2::margin(16, 3, 3, 3)),
axis.title.y = element_text(margin = ggplot2::margin(3, 16, 3, 3)),
legend.position = "inside",
legend.position.inside = c(.987, .95),
legend.justification = c(1, 1),
panel.background = element_rect(fill = "white"),
legend.key = element_rect(fill = "white", color = "white"),
legend.title = element_blank(),
legend.direction = "vertical",
legend.key.width = unit(.3, "in"),
legend.key.height = unit(.3, "in"),
legend.background = element_rect(color = "black", fill = "white", linewidth = .5),
legend.margin = ggplot2::margin(6, 6, 6, 6))
# illustrate a cross section of the point cloud that passes through the center of a plot
centroid = st_centroid(plots[1])
# centroid of the plot, plus and minus 100 m longitude
point1 = c(centroid$geometry[[1]][1] - 150, centroid$geometry[[1]][2])
point2 = c(centroid$geometry[[1]][1] + 150, centroid$geometry[[1]][2])
Warning message:
“st_centroid assumes attributes are constant over geometries”
transect = data.frame(lat = c(point1[2], point2[2]),
lon = c(point1[1], point2[1]),
type = "transect") |>
st_as_sf(coords = c("lon", "lat"), crs = 32617) |>
dplyr::group_by(type) %>%
dplyr::summarise() %>%
st_cast("LINESTRING")
# identify tile that overlaps with transect
path_transect = (readLAScatalog(paste0(root_dir, "/output/_1_classified_points/")) |>
catalog_intersect(transect))$filename
las = readLAS(path_transect)
las_tr = clip_transect(las, point1, point2, width = 2, xz = TRUE)
range_las_tr = max(las_tr$Z) - min(las_tr$Z)
fig = ggplot(payload(las_tr), aes(X, Z, color = factor(Classification))) +
geom_point(size = 1.8, alpha = .7, shape = 20) +
geom_point(data = filter(payload(las_tr), Classification == 18),
size = 5, alpha = .2) +
coord_equal() +
theme_bw(base_size = 15) +
scale_x_continuous(expand = expansion(mult = c(.05, .05)),
breaks = seq(0, 250, 25)) +
scale_y_continuous(breaks = seq(0, 250, 25),
limits = c(min(las_tr$Z) + -15, max(las_tr$Z) + 20)) +
theme_transect +
labs(x = "Transect distance (m)", y = "Elevation (m)") +
scale_color_manual(labels = c("0" = "Unclassified",
"2" = "Ground",
"18" = "Noise"),
values = c("darkseagreen", "red4", "blueviolet"))
ggsave(plot = fig,
filename = "./figures/cross_section_classified.jpeg",
device = jpeg,
width = 14,
height = 5,
units = "in",
dpi = 300,
bg = "white")
# plot classified points
IRdisplay::display_jpeg(file = "./figures/cross_section_classified.jpeg")
A cross section of a point cloud illustrating points classified as ground with the cloth simulation filter.
# plot transect canopy height
path_transect = (readLAScatalog(paste0(root_dir, "/output/_3_normalized_points/")) |>
catalog_intersect(transect))$filename
las = readLAS(path_transect)
las_tr_norm = clip_transect(las, point1, point2, width = 2, xz = TRUE)
fig = ggplot(payload(las_tr_norm), aes(X,Z)) +
geom_point(size = 1.8, alpha = .7, shape = 20, color = "darkseagreen") +
coord_equal() +
theme_bw(base_size = 15) +
scale_x_continuous(expand = expansion(mult = c(.05, .05)),
breaks = seq(0, 250, 50)) +
scale_y_continuous(expand = expansion(mult = c(.15, .15)),
breaks = seq(0, 250, 25),
limits = c(-15,
range_las_tr + 20)) +
theme_transect +
labs(x = "Transect distance (m)", y = "Canopy height (m)") +
guides(legend.position = "none")
ggsave(plot = fig,
filename = "./figures/cross_section_normalized.jpeg",
device = jpeg,
width = 14,
height = 5,
units = "in",
dpi = 300,
bg = "white")
# plot normalized points
IRdisplay::display_jpeg(file = "./figures/cross_section_normalized.jpeg") 
A cross section of a point cloud illustrating points with Z values representing canopy height, after normalization to a digital terrain model.
# define some variables for mapping in ggplot
theme_map = theme(panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
axis.title = element_blank(),
axis.text.x = element_text(angle = 90),
axis.ticks = element_blank(),
legend.position = "inside",
legend.position.inside = c(.98, .02),
legend.justification = c(1, 0),
panel.background = element_rect(fill = "white"),
legend.key = element_rect(fill = "white", color = "white"),
legend.title = element_text(hjust = 0.5),
legend.direction = "horizontal",
legend.key.width = unit(.35, "in"),
legend.key.height = unit(.1, "in"),
legend.background = element_rect(color = "black", fill = "white", linewidth = .5),
legend.margin = ggplot2::margin(4, 10, 4, 10))
cols_gradient <- c(low = "#002070", mid = "#DD6677", high = "#F1F1A1")
cols_gradient_biomass <- c(low = "#002070", mid = "#309040", high = "#D0E358")
cols_gradient_metrics <- c(low = "#300010", mid = "#505060", high = "beige")
preds_50m = rast(paste0(root_dir, "/output/_5_metrics_raster/predictors_50m.tif"))
CHM_smooth = rast(paste0(root_dir, "/output/_4_CHM_raster/CHM_smooth.tif")) |> crop(preds_50m)
DTM_smooth = rast(paste0(root_dir, "/output/_2_DTM_raster/DTM_smooth.tif")) |> crop(preds_50m)
# plot the DTM and CHM for the 50-ha plot
suppressMessages({
p1 = ggR(DTM_smooth |> terra::crop(st_buffer(plots[1,], 140)),
geom_raster = TRUE) +
scale_fill_gradientn(colors = cols_gradient,
na.value = "grey70",
name = "Elevation (m)") +
geom_sf(data = plots[1,], color = "white", fill = NA, linewidth = 1) +
geom_sf(dat = transect, color = "white", linewidth = .8, linetype = 2) +
scale_x_continuous(expand = expansion(mult = c(0, 0))) +
scale_y_continuous(expand = expansion(mult = c(0, 0))) +
theme_bw(base_size = 8) +
theme_map +
guides(fill = guide_colourbar(title.position="top"))
p2 = ggR(CHM_smooth |> terra::crop(st_buffer(plots[1,], 140)),
geom_raster = TRUE) +
scale_fill_gradientn(colors = cols_gradient,
na.value = "grey20",
name = "Canopy height (m)") +
geom_sf(data = plots[1,], color = "white", fill = NA, linewidth = 1) +
geom_sf(dat = transect, color = "white", linewidth = .8, linetype = 2) +
scale_x_continuous(expand = expansion(mult = c(0, 0))) +
scale_y_continuous(expand = expansion(mult = c(0, 0))) +
theme_bw(base_size = 8) +
theme_map +
guides(fill = guide_colourbar(title.position="top"))
fig = ggarrange(p1, p2,
ncol = 2, nrow = 1)
})
ggsave(plot = fig,
filename = "./figures/DTM_CHM_50ha.jpeg",
device = jpeg,
width = 12,
height = 6,
units = "in",
dpi = 300,
bg = "white")
IRdisplay::display_jpeg(file = "./figures/DTM_CHM_50ha.jpeg")
A subset of the digital terrain model and the canopy height model; plot outline shown in white; transect shown as white dotted line.
Biomass and NISAR Datasets Overlap
In this example it is showed how you can query ESA Biomass and NASA NISAR granules to plot the footprint polygons into an interactive Folium map, as two toggleable layers.
It is recommended to use the MAAP environment.
Pre-requisites
from folium import GeoJson
import geopandas as gpd
import pystac_client
import pandas as pd
import foliumInputs Configuration
In the next cell you need to define the parameters to search both catalogs (ESA catalog and NASA catalog) and compute overlaps.
-
BBOX defines the area of interest as (min_lon, min_lat, max_lon, max_lon) and can be used to spatially filter both datasets.
-
NISAR_DT sets the datetime range for the NISAR STAC search (tighten this first to avoid timeouts).
-
BIOMASS_DT sets the datetime range for the BIOMASS STAC search.
-
NISAR_STAC_URL is the STAC endpoint used to query NISAR items (CMR-STAC / ASF).
-
NISAR_COLLECTION is the NISAR collection ID used in the STAC search (e.g.,
NISAR_L2_GSLC_BETA_V1_1). -
BIOMASS_STAC_URL is the STAC endpoint used to query BIOMASS items (ESA MAAP STAC).
-
BIOMASS_COLLECTION is the BIOMASS collection name used in the STAC search (e.g.,
BiomassLevel1b).
# Common query parameters (edit to your AOI/time window)
BBOX = [-180, -90, 180, 90] # [minx, miny, maxx, maxy]
NISAR_DT = "2025-09-01/2025-12-31" # tighten first to avoid timeouts
BIOMASS_DT = "2025-01-01/.." # adjust if needed
NISAR_STAC_URL = "https://cmr.earthdata.nasa.gov/stac/ASF"
NISAR_COLLECTION = "NISAR_L2_GSLC_BETA_V1_1"
BIOMASS_STAC_URL = "https://catalog.maap.eo.esa.int/catalogue/"
BIOMASS_COLLECTION = "BiomassLevel1a"Perform the query to the catalogs and convert the returned items into a GeoDataFrame
Here the returned STAC items are converted into gdf_nisar and gdf_biomass, a GeoDataFrame whose geometry columns contain the true NISAR and Biomass footprint polygons and whose ID/title field is kept for labeling and later joins.
No data files are downloaded: only metadata footprints are used.
nisar_catalog = pystac_client.Client.open(NISAR_STAC_URL)
nisar_search = nisar_catalog.search(
collections=[NISAR_COLLECTION],
bbox=BBOX,
datetime=NISAR_DT,
max_items=500,
)
nisar_items = list(nisar_search.items())
print("NISAR items returned:", len(nisar_items))
# Convert STAC Items → GeoDataFrame
nisar_features = []
for it in nisar_items:
nisar_features.append({
"type": "Feature",
"geometry": it.geometry,
"properties": {
"nisar_id": it.id,
**(it.properties or {}),
},
})
gdf_nisar = gpd.GeoDataFrame.from_features(nisar_features, crs="EPSG:4326")
gdf_nisar = gdf_nisar[["nisar_id", "geometry"]]
gdf_nisar.head()
biomass_catalog = pystac_client.Client.open(BIOMASS_STAC_URL)
biomass_search = biomass_catalog.search(
collections=[BIOMASS_COLLECTION],
bbox=BBOX,
datetime=BIOMASS_DT,
max_items=2000,
method="GET",
)
biomass_items = list(biomass_search.items())
print("BIOMASS items returned:", len(biomass_items))
biomass_features = []
for it in biomass_items:
if it.geometry is None:
raise ValueError(f"Item {it.id} has no geometry (footprint).")
props = it.properties or {}
biomass_features.append({
"type": "Feature",
"geometry": it.geometry,
"properties": {
"biomass_id": it.id,
"start_datetime": props.get("start_datetime"),
"end_datetime": props.get("end_datetime"),
"datetime": props.get("datetime"),
},
})
gdf_biomass = gpd.GeoDataFrame.from_features(biomass_features, crs="EPSG:4326")
gdf_biomass = gdf_biomass[["biomass_id", "start_datetime", "end_datetime", "datetime", "geometry"]]
gdf_biomass.head()
Biomass and NISAR Footprint Layers
Here it is created an interactive Folium map and an overlay of the two GeoDataFrames.
## Ensure both are WGS84 for folium
gdf_nisar = gdf_nisar.to_crs("EPSG:4326")
gdf_biomass = gdf_biomass.to_crs("EPSG:4326")
# Center/zoom using combined bounds
combined = gpd.GeoSeries(list(gdf_nisar.geometry) + list(gdf_biomass.geometry), crs="EPSG:4326")
minx, miny, maxx, maxy = combined.total_bounds
center = [(miny + maxy) / 2, (minx + maxx) / 2]
m = folium.Map(location=center, tiles="OpenStreetMap", zoom_start=3)
def style_nisar(_):
return {"color": "#1f77b4", "weight": 2, "fillOpacity": 0.15} # blue
def style_biomass(_):
return {"color": "#ff7f0e", "weight": 1, "fillOpacity": 0.10} # orange
GeoJson(
gdf_nisar.__geo_interface__,
name=f"NISAR ({len(gdf_nisar)})",
tooltip=folium.GeoJsonTooltip(fields=["nisar_id"]),
style_function=style_nisar,
).add_to(m)
GeoJson(
gdf_biomass.__geo_interface__,
name=f"BIOMASS ({len(gdf_biomass)})",
tooltip=folium.GeoJsonTooltip(fields=["biomass_id"]),
style_function=style_biomass,
).add_to(m)
folium.LayerControl(collapsed=False).add_to(m)
m.fit_bounds([[miny, minx], [maxy, maxx]])
m
Find the overlap
To identify which Biomass footprint polygons intersect which NISAR footprint polygons, you can use GeoPandas spatial join. Then, the total number of intersection pairs found is printed out.
# Spatial join to find overlapping pairs
pairs = gpd.sjoin(
gdf_nisar.reset_index(drop=True),
gdf_biomass.reset_index(drop=True),
how="inner",
predicate="intersects",
)
print("Overlap pairs:", len(pairs))
# Compact table of matches (deduped)
matches = pairs[["biomass_id", "nisar_id"]].drop_duplicates().reset_index(drop=True)
matches.head(25)
Finally, you need to create the actual overlap polygons and then visualise only those overlaps on a map.
right_geom = gdf_biomass.loc[pairs["index_right"], "geometry"].values
right_geom = gpd.GeoSeries(right_geom, index=pairs.index, crs="EPSG:4326")
overlap_geom = pairs.geometry.intersection(right_geom)
overlap = gpd.GeoDataFrame(
pairs[["nisar_id", "biomass_id"]].copy(),
geometry=overlap_geom,
crs="EPSG:4326",
)
overlap = overlap[~overlap.geometry.is_empty]
# Map: overlap only
m_overlap = folium.Map(tiles="OpenStreetMap")
def style_overlap(_):
return {"color": "#2ca02c", "weight": 2, "fillColor": "#2ca02c", "fillOpacity": 0.35}
GeoJson(
overlap.__geo_interface__,
name=f"Overlap ({len(overlap)})",
tooltip=folium.GeoJsonTooltip(fields=["nisar_id", "biomass_id"]),
style_function=style_overlap,
).add_to(m_overlap)
if len(overlap) > 0:
minx, miny, maxx, maxy = overlap.total_bounds
pad = 0.05
m_overlap.fit_bounds([[miny - pad, minx - pad], [maxy + pad, maxx + pad]])
m_overlap
No comments to display
No comments to display