library(tidyverse)
library(sf)
library(terra)
library(leaflet)
library(rstac)Mini-project: Regreening App × STAC
From restoration observations to environmental evidence
In this notebook, you will learn how to combine field data from the Regreening App with spatial raster data from the GGW STAC.
If you are new to R, or want extra practice, we encourage checking out Notebook 1 - Basics of R and analysis with climate datasets, or some of the other notebooks available on the Practicals Gallery of the K4GGWA Platform’s Evidence Hub.
You will:
- download your Regreening App project’s polygon data from the Data Reporting System (DRS);
- check and explore the data in R;
- turn the wkt geometry in the csv into spatial polygons;
- access rasters from the GGW STAC;
- extract raster values within your restoration polygons;
- make a simple visualisation and write a short interpretation.
1. Get your Regreening App data
Download your project data
Open the Regreening App Data Reporting System (DRS).
As a Project Manager (PM):
- Open the DRS.
- Select the project you want to analyse.
- Download the project data as a csv file.
- Save the csv somewhere easy to find (a project folder).
The exported file contains information about the restoration plot, trees recorded within it, management practices, uses, location and polygon geometry. The kinds of data available depend on the app module that was used.
A note about the geometry
The Plot Geometry column contains polygons stored as wkt (Well-Known Text).
For example:
POLYGON((41.0162303 11.7928467,
41.0161877 11.7928214,
41.0161457 11.7928436,
41.0160885 11.7928373,
41.0162303 11.7928467))The coordinates are stored as:
longitude latitudeWe will turn these strings into spatial polygons using sf.
2. Set up R
Choose your downloaded csv:
dat <- read_csv("data.csv", show_col_types = FALSE)
glimpse(dat)Rows: 184
Columns: 36
$ `Uploaded Date` <dttm> 2024-05-23 19:21:25, 2024-05-23 1…
$ `Date Collected` <date> 2024-05-23, 2024-05-23, 2026-04-2…
$ `Data Collector` <chr> "Abdu Ali", "Abdu Ali", "Girmay We…
$ `Project Name` <chr> "Regreening Africa", "Regreening A…
$ `Farmer/Group/Institution` <chr> "Abdu Husein", "Abdu Husein", "gir…
$ Country <chr> "Ethiopia", "Ethiopia", "Ethiopia"…
$ `County/Region/Provence` <chr> "Afar", "Afar", "Tigray", "Tigray"…
$ `Subcounty/District/commune/worede` <chr> "Samara", "Samara", "south east En…
$ `Plot Name` <chr> "2bab3532-cea9-4011-b699-a76e66903…
$ `Date FMNR Started` <date> 1987-05-01, 1987-05-01, 2026-04-0…
$ `Plot has Crops` <lgl> FALSE, FALSE, FALSE, TRUE, FALSE, …
$ `Crop Name` <chr> NA, NA, NA, "Tiff", NA, NA, NA, NA…
$ `Plot has fence` <lgl> TRUE, TRUE, TRUE, FALSE, TRUE, TRU…
$ `Fence type` <chr> "LIVE", "LIVE", "LIVE", "SOCIAL_FE…
$ `Other Fence type` <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `Plot Ownership type` <chr> "INDIVIDUAL - CHURCH_OR_MOSQUE - C…
$ `Plot other ownership type` <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `Plot Size` <dbl> 1.35, 1.35, 0.50, 10.60, 0.12, 13.…
$ `Plot Size Unit` <chr> "Hectare", "Hectare", "Hectare", "…
$ `Plot Actual Size` <dbl> 0.00, 0.00, 0.03, 10.06, 0.18, 16.…
$ `Plot Geometry` <chr> "POLYGON((41.0162303 11.7928467,41…
$ `Local Name` <chr> "Burtukan", "Papaye", "seraw", "ac…
$ `Scientific Name` <chr> "Abutilon mauritianum", "Acacia ab…
$ `Management Practices` <chr> "COMPOST_MANURE COPPICING WEEDING …
$ `Other Managements` <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA…
$ `Tree Usage` <chr> "NUTS FODDER SOIL_FERTILITY LEAF_V…
$ `Other Usages` <chr> NA, NA, NA, NA, NA, "fencing", NA,…
$ Height <chr> "LESS_1_5", "BETWEEN_1_5_AND_3", "…
$ `Actual Height` <chr> "50.00", "None", "3.50", "5.00", "…
$ `RCC/CBH` <chr> "1.35", "1.35", "6.00", "70.00", "…
$ Latitude <chr> "11.792856", "11.79289", "13.47751…
$ Longitude <chr> "41.016057", "41.015948", "39.4550…
$ Altitude <chr> "417.318", "417.508", "2137.500", …
$ Accuracy <chr> "4.17", "3.90", "10.48", "1.87", "…
$ Photo <chr> "https://radrs.icraf.org/uploaded/…
$ Comment <chr> "Very important of this trees", "A…# Let's do a quick data check
dim(dat)[1] 184 36dat |>
summarise(
records = n(),
unique_plots = n_distinct(`Plot Name`),
unique_species = n_distinct(`Scientific Name`, na.rm = TRUE),
missing_geometry = sum(is.na(`Plot Geometry`))
)# A tibble: 1 × 4
records unique_plots unique_species missing_geometry
<int> <int> <int> <int>
1 184 52 92 0A single restoration plot can contain multiple tree records, so the number of rows in the csv may be greater than the number of unique polygons.
Check this:
dat |>
count(`Plot Name`) |>
slice_head(n = 10)# A tibble: 10 × 2
`Plot Name` n
<chr> <int>
1 0084a325-00da-41c6-a8b0-8ab2914a594b 2
2 02db3658-0625-4124-90f7-8e6ff972f113 5
3 0679f72b-8b63-400c-a59d-8872fbd0f09a 1
4 074bba73-b9d7-4b7e-82cf-05cdb81dacda 2
5 0fed5cd8-0658-43bf-b4a1-a00eee6bec4d 1
6 2a9db6e0-deb1-4103-ad72-2f9304b8386d 2
7 2af52e21-6ad6-46ae-95fd-9a05d23c4d22 1
8 2b1ae5df-3a87-4299-ac13-08b2ab361436 7
9 2bab3532-cea9-4011-b699-a76e66903f4e 2
10 2c63b49b-7183-45a9-8dc8-29eb7c23e674 3This distinction matters when calculating plot-level statistics.
3. Explore the restoration data
Which tree species have been recorded?
Try a simple bar chart of the most frequently recorded species.
dat |>
filter(!is.na(`Scientific Name`)) |>
count(`Scientific Name`, sort = TRUE) |>
slice_head(n = 20) |>
ggplot(aes(x = reorder(`Scientific Name`, n), y = n)) +
geom_col() +
coord_flip() +
labs(
title = "Most frequently recorded tree species",
x = NULL,
y = "Number of records"
) +
theme_minimal()
How does the data look?
Often, since we are working with citizen-science data, there may be spelling errors, duplications, or other errors we should be aware of, and manage appropriately.
Which management practices are being used?
The Management Practices field can contain several practices in one cell, for example:
PRUNING FENCING WEEDINGWe can separate these into individual practices.
management <- dat |>
filter(!is.na(`Management Practices`)) |>
separate_rows(`Management Practices`, sep = " ") |>
count(`Management Practices`, sort = TRUE)
management# A tibble: 13 × 2
`Management Practices` n
<chr> <int>
1 PRUNING 165
2 FENCING 64
3 THINNING 48
4 COMPOST_MANURE 47
5 WEEDING 45
6 WATERING 44
7 MULCHING 33
8 COPPICING 25
9 POLLARDING 25
10 FIRE 16
11 CHEMICAL_FERTILIZER 13
12 OTHER 2
13 None 1Visualise the most common practices:
management |>
slice_head(n = 10) |>
ggplot(aes(
x = reorder(`Management Practices`, n),
y = n
)) +
geom_col() +
coord_flip() +
labs(
title = "Recorded management practices",
x = NULL,
y = "Number of records"
) +
theme_minimal()
Mini-question
Write down one observation from your exploratory analysis.
For example:
Pruning appears frequently in the monitoring records, while several other practices occur much less often.
Then ask:
Is this pattern a real difference in management, or could it partly reflect how data were collected?
4. Turn the WKT into polygons
Now we move from ordinary tabular data to spatial data.
First, keep one record for each unique restoration plot:
plots <- dat |>
distinct(`Plot Name`, `Plot Geometry`, .keep_all = TRUE) |>
filter(!is.na(`Plot Geometry`))Convert the WKT geometry to an sf object:
plots_sf <- st_as_sf(
plots,
wkt = "Plot Geometry",
crs = 4326
)
plots_sfSimple feature collection with 52 features and 35 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 38.95715 ymin: 11.79282 xmax: 41.01623 ymax: 13.69444
Geodetic CRS: WGS 84
# A tibble: 52 × 36
`Uploaded Date` `Date Collected` `Data Collector` `Project Name`
* <dttm> <date> <chr> <chr>
1 2024-05-23 19:21:25 2024-05-23 Abdu Ali Regreening Africa
2 2026-04-29 12:24:10 2026-04-29 Girmay Welay Regreening Africa
3 2026-05-05 16:51:29 2026-05-05 Abrah Hadgu Regreening Africa
4 2026-05-16 18:16:50 2026-05-14 Haftu Tsegy Regreening Africa
5 2026-05-18 11:38:05 2026-05-18 Gebremedhn Hagos Regreening Africa
6 2026-05-20 18:45:06 2026-05-20 Goitom Haileslasie Regreening Africa
7 2026-05-20 19:13:59 2026-05-20 girmay welay Regreening Africa
8 2026-05-21 18:23:43 2026-05-21 Azab Araya Regreening Africa
9 2026-05-21 19:50:23 2026-05-21 girmay welay Regreening Africa
10 2026-05-20 18:28:51 2026-05-21 Shishay Hagos Regreening Africa
# ℹ 42 more rows
# ℹ 32 more variables: `Farmer/Group/Institution` <chr>, Country <chr>,
# `County/Region/Provence` <chr>, `Subcounty/District/commune/worede` <chr>,
# `Plot Name` <chr>, `Date FMNR Started` <date>, `Plot has Crops` <lgl>,
# `Crop Name` <chr>, `Plot has fence` <lgl>, `Fence type` <chr>,
# `Other Fence type` <chr>, `Plot Ownership type` <chr>,
# `Plot other ownership type` <chr>, `Plot Size` <dbl>, …Check the geometries:
sum(!st_is_valid(plots_sf))[1] 7st_crs(plots_sf)Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]See your restoration plots
leaflet(plots_sf) |>
addProviderTiles("OpenStreetMap.Mapnik") |>
addPolygons(
weight = 2,
fillOpacity = 0.3
)You should now be able to see the polygons you collected in the field.
5. Connect to the GGW STAC
The GGW STAC provides access to spatial raster datasets that can be used alongside the restoration observations.
stac_source <- stac("https://stacapi100.thegrit.earth/")
stac_source###rstac_query
- url: https://stacapi100.thegrit.earth/
- params:
- field(s): version, base_url, endpoint, params, verb, encodeList the available collections:
collections_query <- stac_source |>
rstac::collections()
available_collections <- rstac::get_request(collections_query)
available_collections###Collections
- collections (13 item(s)):
- ls8_250_Africa
- eswatini
- flood_analysis_ggw
- tree_suitability_africa_1km
- africa_climate
- africa_soil_maps_250m
- srilanka_landhealth
- ggw
- somalia
- harmony_in_habitat
- ukpact_landhealth
- SRTM
- twende_landhealth
- field(s): collections, linksThe STAC catalogue contains datasets that can provide environmental context for restoration sites.
6. Find raster data around your project
We can use the extent of our restoration polygons to define a bounding box.
project_bbox <- st_bbox(plots_sf)
project_bbox xmin ymin xmax ymax
38.95715 11.79282 41.01623 13.69444 Convert it to the order expected by STAC:
bbox <- c(
project_bbox["xmin"],
project_bbox["ymin"],
project_bbox["xmax"],
project_bbox["ymax"]
)
bbox xmin ymin xmax ymax
38.95715 11.79282 41.01623 13.69444 Now we can search the STAC.
For example, the reference STAC notebook uses the africa_soil_maps_250m collection:
items_query <- rstac::stac_search(
q = stac_source,
collections = "africa_soil_maps_250m",
bbox = bbox,
limit = 20
)
available_items <- rstac::get_request(items_query)
available_items###Items
- matched feature(s): 8
- features (8 item(s) / 0 not fetched):
- 20b98b663dce4a51bf81127b2c096e6a
- 215917e3cf4d4e32a2206e5fec3cc4c8
- 6b60d0bae25841398ac8c7ea5c04ee77
- 7726efadd9ac48ef8ca804a7a7ad70c5
- 92da9f22bec64b1587de76efb4b81460
- c1fc237fd08e4a249277bae7069211ff
- cbc0ff0331d8465789d323c5b7af9929
- d4568b525a804a74b6be8d7b8b6845d3
- assets: layer, thumbnail
- item's fields:
assets, bbox, collection, geometry, id, links, properties, stac_version, typeYou can inspect the available items/assets:
for (item in available_items$features) {
print(names(item$assets))
}[1] "layer" "thumbnail"
[1] "layer" "thumbnail"
[1] "layer" "thumbnail"
[1] "layer" "thumbnail"
[1] "layer" "thumbnail"
[1] "layer" "thumbnail"
[1] "layer" "thumbnail"
[1] "layer" "thumbnail"7. Read a raster from the STAC
For this exercise we will use the Soil Organic Carbon raster from the K4GGWA STAC.
The COG URL used in the reference notebook is:
soc_url <- "https://stacapi100.thegrit.earth/eodata/cogeo/ggw/Soil_Organic_Carbon_2023_GGW_500m.tif"A Cloud Optimised GeoTIFF (COG) can be accessed remotely without downloading the entire raster.
Create a SpatRaster:
soc <- terra::rast(
paste0("/vsicurl/", soc_url)
)
socclass : SpatRaster
size : 5319, 14584, 1 (nrow, ncol, nlyr)
resolution : 0.004491576, 0.004491576 (x, y)
extent : -17.5441, 47.96105, 3.395632, 27.28633 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source : Soil_Organic_Carbon_2023_GGW_500m.tif
name : Soil_Organic_Carbon_2023_GGW_500m
min value : -9999
max value : 12759Check its basic properties:
list(
crs = crs(soc),
resolution = res(soc),
extent = ext(soc),
layers = nlyr(soc)
)$crs
[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
$resolution
[1] 0.004491576 0.004491576
$extent
SpatExtent : -17.544097498854256, 47.961053019141261, 3.3956317739717932, 27.286326755130467 (xmin, xmax, ymin, ymax)
$layers
[1] 18. Overlay the restoration plots and raster
Before extracting raster values, make sure the polygons and raster use compatible coordinate reference systems.
plots_vect <- terra::vect(plots_sf)
plots_vect <- terra::project(
plots_vect,
soc
)
crs(plots_vect)[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"crs(soc)[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"Now extract the raster values within each restoration polygon.
soc_stats <- terra::extract(
soc,
plots_vect,
fun = mean,
na.rm = TRUE
)
head(soc_stats) ID Soil_Organic_Carbon_2023_GGW_500m
1 1 752
2 2 1046
3 3 2030
4 4 1885
5 5 2372
6 6 2039Each polygon now has a summary of the raster values found inside it.
9. Your first raster + vector analysis
Join the extracted values back to the plot data:
plots_with_soc <- plots_sf |>
mutate(
SOC_mean = soc_stats[, 2]
)
plots_with_soc |>
st_drop_geometry() |>
select(
`Plot Name`,
`Scientific Name`,
`Management Practices`,
SOC_mean
) |>
head()# A tibble: 6 × 4
`Plot Name` `Scientific Name` `Management Practices` SOC_mean
<chr> <chr> <chr> <dbl>
1 2bab3532-cea9-4011-b699-a76… Abutilon mauriti… COMPOST_MANURE COPPIC… 752
2 e920412d-6c4a-491a-b1be-05c… Acacia abyssinica WATERING PRUNING COMP… 1046
3 0679f72b-8b63-400c-a59d-887… Acacia abyssinica COPPICING PRUNING WEE… 2030
4 86688712-7781-48d0-a04d-d4c… acaciaetbaicatum PRUNING FENCING 1885
5 86d496c6-3622-4a62-ac9e-97b… Acacia abyssinica PRUNING 2372
6 3db125ef-5401-47ce-a86b-5bc… delonix regia FENCING PRUNING 2039Challenge 1 — Describe your restoration landscape
Make a histogram of the mean SOC value across your restoration plots.
plots_with_soc |>
st_drop_geometry() |>
ggplot(aes(x = SOC_mean)) +
geom_histogram(bins = 20) +
labs(
title = "Soil organic carbon across restoration plots",
x = "Mean soil organic carbon",
y = "Number of restoration plots"
) +
theme_minimal()
Write a short interpretation
In 2–3 sentences:
What range of environmental conditions do the restoration plots occupy?
Are most plots clustered around similar values, or is there substantial variation?
10. Challenge 2 — Compare environmental conditions between plots
Now explore whether restoration plots with different characteristics occur under different environmental conditions.
For example, compare plots with and without crops.
plots_with_soc |>
st_drop_geometry() |>
ggplot(
aes(
x = `Plot has Crops`,
y = SOC_mean
)
) +
geom_boxplot() +
labs(
title = "Soil organic carbon by cropping status",
x = "Plot has crops",
y = "Mean soil organic carbon"
) +
theme_minimal()
Interpret carefully
You are comparing observed environmental conditions, not demonstrating that crops caused differences in SOC.
Write:
What pattern do you see?
What are two possible explanations?
What additional information would you need before making a causal claim?
11. Challenge 3 — Map an environmental variable
A map can often communicate the relationship between restoration activities and environmental conditions more effectively than a table.
Map the restoration polygons and colour them according to their extracted SOC value.
pal <- colorNumeric(
palette = "YlGn",
domain = plots_with_soc$SOC_mean,
na.color = "transparent"
)
leaflet(plots_with_soc) |>
addProviderTiles("OpenStreetMap.Mapnik") |>
addPolygons(
fillColor = ~pal(SOC_mean),
fillOpacity = 0.7,
color = "white",
weight = 1,
popup = ~paste0(
"<strong>Plot:</strong> ", `Plot Name`,
"<br><strong>Mean SOC:</strong> ",
round(SOC_mean, 2)
)
) |>
addLegend(
pal = pal,
values = ~SOC_mean,
title = "Mean SOC"
)Your task
Modify the map so that it shows a different variable or a different raster from the K4GGWA STAC.
Possible questions include:
- Do restoration plots occur across a wide range of soil conditions?
- Are plots with crops concentrated in particular environmental conditions?
- Do plots with fencing show a different spatial pattern?
- Are particular tree species recorded in particular environmental settings?
12. Challenge 4 — Combine field observations with raster evidence
Choose one field variable and one environmental variable.
For example:
| Field observation | Environmental variable |
|---|---|
| Crops present | Soil organic carbon |
| Fencing | Soil organic carbon |
| Management practice | Soil organic carbon |
| Tree species | Soil organic carbon |
| Plot size | Soil organic carbon |
Create one clear figure.
For example, you could compare SOC between plots with and without fencing:
final_plot <- plots_with_soc |>
st_drop_geometry() |>
ggplot(
aes(
x = `Plot has fence`,
y = SOC_mean
)
) +
geom_boxplot() +
labs(
title = "Environmental conditions across restoration plots",
subtitle = "Soil organic carbon by fencing status",
x = "Plot has fence",
y = "Mean soil organic carbon"
) +
theme_minimal()
final_plot
Save your figure:
ggsave(
"regreening_stac_mini_project.png",
final_plot,
width = 8,
height = 5,
dpi = 300
)13. Write your evidence note
Write 100–150 words answering:
What did you investigate?
What restoration characteristic did you examine?
What environmental conditions did you observe?
What did the STAC raster tell you about the landscape?
What pattern did you find?
Describe the main pattern in your figure.
What might explain it?
Suggest one or two plausible ecological, climatic or land-use explanations.
What can’t you conclude?
Remember:
An association between a restoration characteristic and a raster-derived environmental variable does not by itself demonstrate causation.
Your analysis is an exploratory use of spatial evidence.
14. Looking ahead: analysing change through time
The workflow becomes even more powerful when the STAC contains time-indexed raster data.
Instead of having one raster representing a single period, a STAC collection can contain observations for different dates or time intervals.
Conceptually:
Restoration polygon
│
├── Raster: 2010
│
├── Raster: 2015
│
├── Raster: 2020
│
└── Raster: 2025You can then extract the raster value for the same polygon at each time step.
For example:
Plot 1
2010 12.4
2015 13.1
2020 15.7
2025 17.2This allows questions such as:
- Has vegetation productivity changed inside restoration plots?
- Has soil moisture changed?
- Has land surface temperature changed?
- How has a climate indicator changed before and after restoration?
- Do different restoration approaches occur in areas with different environmental trajectories?
The STAC search can eventually use a time interval, for example:
items_query <- rstac::stac_search(
q = stac_source,
collections = "YOUR_TIME_ENABLED_COLLECTION",
bbox = bbox,
datetime = "2010-01-01/2025-12-31",
limit = 100
)You would then loop over the returned time-indexed items, extract values for the same polygons, and build a plot × time dataset.
Note
Coming soon: As more time-indexed raster products are added to the K4GGWA STAC, this workflow can be extended from describing current environmental conditions to tracking change within restoration areas over time.
15. Your final mini-project
Your final submission should contain:
1. A spatial dataset
Your Regreening App polygons joined with at least one STAC-derived environmental variable.
2. One polished visual
A map, boxplot, histogram, scatterplot or another clear figure.
3. A short evidence note
100–150 words explaining:
What did you investigate → What did you find → Why might it matter → What are the limitations?
4. One idea for using the result
Think beyond the notebook.
Could this analysis help you:
- identify where restoration is taking place under particular environmental conditions?
- compare projects or landscapes?
- support project reporting?
- identify areas for further monitoring?
- communicate restoration progress?
- develop a figure for a report or presentation?
- create a short evidence story for colleagues or stakeholders?
Key takeaway
The important skill is not simply knowing how to make a map or run terra::extract().
It is learning to connect:
What was observed in the field
↓
Where it happened
↓
What environmental conditions occurred there
↓
What pattern the data show
↓
What that evidence can — and cannot — tell us
That is the foundation for using restoration monitoring data together with Earth observation and spatial environmental datasets.