FLEX MAAP
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.
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 "BiomassSimulated"
- Acquisition date between 2017-01-01 and 2017-12-01 (YYYY-MM-dd)
- With bounding box [112.291824, -3.33663, 112.94324, -2.66676]
from pystac_client import Client
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 = 'BiomassSimulated',
bbox = [112.291824, -3.33663, 112.94324, -2.66676],
datetime = ['2017-01-01T00:00:00Z', '2017-12-01T00:00:00Z'] # Search for granules in date range
)
print(f'{len(results.item_collection_as_dict()['features'])} granules found')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 10 hours. Then, you will be asked to log in with your credentials to generate a new token.
Below you can find a sample of code that shows you how to use your token.
access_token = ' ' # Copy here your token
data = results.item_collection_as_dict()
n_products = 5 # 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.