New Caledonia with loss and gain

Get forest cover change from TMF

The function .get_fcc_loss_gain() can be used to download forest cover change from the Tropical Moist Forest product.

This function considers both forest loss and gain (or regrowth) to derive the forest cover change map.

We will use New Caledonia as a case study.

import os
from pathlib import Path
import time

import ee
import geefcc
import numpy as np
import pandas as pd
from tabulate import tabulate
from osgeo import gdal
import folium
from dotenv import load_dotenv

gdal.UseExceptions()

# Initialize GEE
ee.Initialize(project="deforisk", opt_url="https://earthengine-highvolume.googleapis.com")

We want to estimate and map the forest cover change for the period 2015–2025, considering regrowth of at least 10 years as being forest, which is debatable [Bourgoin et al., 2024, Poorter et al., 2016].

# Download data from GEE
min_years = 10
out_dir = Path(f"out_tmf_{min_years}yr")
ofile = out_dir / "fcc_tmf.tif"
if not ofile.is_file():
    start_time = time.time()
    geefcc.get_fcc_loss_gain(
        aoi=(163.5, -23, 168.15, -19.51),
        year1=2015,
        year2=2025,
        min_years=min_years,
        tile_size=1,
        crop_to_aoi=True,
        parallel=True,
        output_file=ofile,
    )
    end_time = time.time()

We estimate the computation time to download 20 1-degree tiles using several cores.

elapsed_time = (end_time - start_time) / 60
print('Execution time:', round(elapsed_time, 2), 'minutes')

Execution time: 1.25 minutes

Plot the forest cover change map

We plot the forest cover change map. The raster is automatically resampled to a coarser resolution before plotting as it exceeds the default pixel threshold. Country borders and the tile grid are overlaid directly from their vector files.

geefcc.plot_fcc_loss_gain(
    input_file=ofile,
    output_file=f"fcc_tmf_{min_years}yr.png",
    title="Forest cover change 2015 \u2013 2025, TMF",
    dpi=200,
    borders=Path("data") / "borders_NCL.gpkg",
    grid=out_dir / "grid.gpkg",
    xlim=(163, 169),
    ylim=(-23.25, -18.75),
)

Lines in black represent country borders. One degree tiles in grey cover the whole buffer and were used to download the data in parallel.

Interactive map

Reproject in Web Mercator

ifile = out_dir / "fcc_tmf.tif"
ofile = out_dir / "fcc_tmf_epsg3857.tif"
ds = gdal.Warp(
    str(ofile),
    str(ifile),
    xRes=30, yRes=30,
    targetAlignedPixels=True,   # aligne la grille sur des multiples de 30 m
    dstSRS="EPSG:3857",
    resampleAlg="near",
    format="GTiff",
    creationOptions=["COMPRESS=DEFLATE", "TILED=YES"],
)
ds = None                  # ferme le dataset et force l'écriture sur disque

Create tiles

import subprocess
subprocess.run([
    "gdal2tiles.py", "--xyz", "-p", "mercator",
    "-z", "5-12", "-r", "near", "--processes=4", "-w", "none",
    str(ofile), "map/tiles_fcc"),
], check=True)

Interactive map with folium

# Mapbox token
# Add MAPBOX_PUBLIC_TOKEN=<your token> in a .env file
load_dotenv()
mapbox_token = os.getenv("MAPBOX_PUBLIC_TOKEN")
mapbox_token = "MAPBOX_PUBLIC_TOKEN"  # Uncomment for GitHub

# Mapbox
base = folium.TileLayer(
    tiles=("https://api.mapbox.com/styles/v1/mapbox/light-v11/"
            f"tiles/{{z}}/{{x}}/{{y}}?access_token={mapbox_token}"),
    name="Mapbox",
    attr="© Mapbox, © <a href='https://www.openstreetmap.org/copyright'>OpenStreetMap</a> contributors",
    max_zoom=14
)

# Add tiles
tile_layer = folium.TileLayer(
    tiles=("tiles_fcc/{z}/{x}/{y}.png"),
    overlay=True,
    name="FCC TMF 2015 \u2013 2025",
    min_zoom=5,
    max_zoom=14,
    max_native_zoom=12,
    attr="EC-JRC"
)

# Create empty map
m = folium.Map(
    tiles=None,
    location=[-21.3, 166.2],
    zoom_start=8,
    max_zoom=14
)
tile_layer.add_to(m)
base.add_to(m)

# Ajouter layer control
folium.LayerControl().add_to(m)

# Save to html
m.save("map/fcc_tmf_interactive.html")

Area per class of forest cover change

We use function stat_fcc_loss_gain() to reproject the raster and compute the number of pixels per class and the corresponding area (in ha). We use projection UTM zone 58S (EPSG code 32758) for New Caledonia.

res_df = geefcc.stat_fcc_loss_gain(
    input_file=ofile,
    epsg=32758,
    output_file=f"fcc_statistics_{min_years}yr.csv",
)
tabulate(res_df, headers=res_df.columns, tablefmt="orgtbl", showindex=False)

category	label	count	area_ha
0	stable non-forest	201155951	18104036
1	stable forest	9349449	841450
2	forest –> deforested	177846	16006
3	non-forest –> old regrowth	69727	6275
4	forest –> old regrowth (via deforestation)	0	0
5	stable old-regrowth	12317	1109
6	old regrowth –> deforested	1522	137
7	stable old afforestation	187208	16849
8	old afforestation –> deforested	51	5
9	non-forest –> old afforestation	412777	37150

Deforestation and regrowth estimates

We can then estimate gross loss, gross gain and net loss in forest cover change for the period 2015–2025. Category 4 (forest –> old regrowth via deforestation) is accounted for both gross loss and gain.

forest_t1 = res_df.loc[[1, 2, 4, 5, 6, 7, 8], "area_ha"].sum()
forest_t2 = res_df.loc[[1, 3, 4, 5, 7, 9], "area_ha"].sum()
gross_loss = - res_df.loc[[2, 4, 6, 8], "area_ha"].sum()
gross_gain = res_df.loc[[3, 4, 9], "area_ha"].sum()
lossgain_df = pd.DataFrame({
    "label": ["forest_t1", "forest_t2", "gross loss", "gross gain", "net change"],
    "area_ha": [forest_t1, forest_t2, gross_loss, gross_gain, gross_gain + gross_loss],
})
time = 2025 - 2015
lossgain_df["annual_change_ha"] = (lossgain_df["area_ha"] / time).round().astype(int)
ratio = lossgain_df["area_ha"] / forest_t1
lossgain_df["annual_change_perc"] = round(100 * (1 - pow((1 - ratio), 1 / time)), 2)
lossgain_df.iloc[:2, 2:4] = np.nan

# Export
lossgain_df.to_csv(f"loss_gain_statistics_{min_years}yr.csv", index=False)
tabulate(lossgain_df, headers=lossgain_df.columns, tablefmt="orgtbl", showindex=False)

label	area_ha	annual_change_ha	annual_change_perc
forest_t1	875556	nan	nan
forest_t2	902833	nan	nan
gross loss	-16148	-1615	-0.18
gross gain	43425	4342	0.51
net change	27277	2728	0.32

When considering regrowth of at least 10 years, which is short for forest recovery [Bourgoin et al., 2024], the gain (4342 ha/yr) compensates the forest cover loss (-1615 ha/yr), and the net change is positive (2728 ha/yr, corresponding to 0.32 %/yr).

When disregarding afforestation (categories 8 and 9, which might be several things including missed forest in the past), figures are different:

forest_t1 = res_df.loc[[1, 2, 4, 5, 6, 7], "area_ha"].sum()
forest_t2 = res_df.loc[[1, 3, 4, 5, 7], "area_ha"].sum()
gross_loss = - res_df.loc[[2, 4, 6], "area_ha"].sum()
gross_gain = res_df.loc[[3, 4], "area_ha"].sum()
lossgain_df = pd.DataFrame({
    "label": ["forest_t1", "forest_t2", "gross loss", "gross gain", "net change"],
    "area_ha": [forest_t1, forest_t2, gross_loss, gross_gain, gross_gain + gross_loss],
})
time = 2025 - 2015
lossgain_df["annual_change_ha"] = (lossgain_df["area_ha"] / time).round().astype(int)
ratio = lossgain_df["area_ha"] / forest_t1
lossgain_df["annual_change_perc"] = round(100 * (1 - pow((1 - ratio), 1 / time)), 2)
lossgain_df.iloc[:2, 2:4] = np.nan

# Export
lossgain_df.to_csv(f"loss_gain_statistics_{min_years}yr_noaff.csv", index=False)
tabulate(lossgain_df, headers=lossgain_df.columns, tablefmt="orgtbl", showindex=False)

label	area_ha	annual_change_ha	annual_change_perc
forest_t1	875551	nan	nan
forest_t2	865683	nan	nan
gross loss	-16143	-1614	-0.18
gross gain	6275	628	0.07
net change	-9868	-987	-0.11

In this case, the gain (628 ha/yr) does not compensate for the loss (-1614 ha/yr), and the net change is negative (-987 ha/yr, corresponding to -0.11 %/yr).

References

[1] (1,2)

C. Bourgoin, G. Ceccherini, M. Girardello, C. Vancutsem, V. Avitabile, P. S. A. Beck, R. Beuchle, L. Blanc, G. Duveiller, M. Migliavacca, G. Vieilledent, A. Cescatti, and F. Achard. Human degradation of tropical moist forests is greater than previously estimated. Nature, 631(8021):570–576, July 2024. doi:10.1038/s41586-024-07629-0.

[2]

L. Poorter, F. Bongers, T. M. Aide, A. M. Almeyda Zambrano, P. Balvanera, J. M. Becknell, V. Boukili, P. H. S. Brancalion, E. N. Broadbent, R. L. Chazdon, D. Craven, J. S. de Almeida-Cortez, G. A. L. Cabral, B. H. J. de Jong, J. S. Denslow, D. H. Dent, S. J. DeWalt, J. M. Dupuy, S. M. Durán, M. M. Espírito-Santo, M. C. Fandino, R. G. César, J. S. Hall, J. L. Hernandez-Stefanoni, C. C. Jakovac, A. B. Junqueira, D. Kennard, S. G. Letcher, J.-C. Licona, M. Lohbeck, E. Marín-Spiotta, M. Martínez-Ramos, P. Massoca, J. A. Meave, R. Mesquita, F. Mora, R. Muñoz, R. Muscarella, Y. R. F. Nunes, S. Ochoa-Gaona, A. A. de Oliveira, E. Orihuela-Belmonte, M. Peña-Claros, E. A. Pérez-García, D. Piotto, J. S. Powers, J. Rodríguez-Velázquez, I. E. Romero-Pérez, J. Ruíz, J. G. Saldarriaga, A. Sanchez-Azofeifa, N. B. Schwartz, M. K. Steininger, N. G. Swenson, M. Toledo, M. Uriarte, M. van Breugel, H. van der Wal, M. D. M. Veloso, H. F. M. Vester, A. Vicentini, I. C. G. Vieira, T. V. Bentos, G. B. Williamson, and D. M. A. Rozendaal. Biomass resilience of neotropical secondary forests. Nature, 530(7589):211–214, February 2016. doi:10.1038/nature16512.