Forest cover fraction (CGLS-LC100)

The Copernicus Global Land Service LC100 collection 3 gives, for every 100 m pixel, the fraction of the ground covered by tree canopy, as a percentage. Epoch 2019 is the last the service produced. Where a land-cover class says “tree cover” or not, this says how much, which is what an optical snow product needs: snow under a closed canopy is snow a satellite cannot see.

The default route (source="zenodo") reads one global GeoTIFF on Zenodo with range requests and needs no account. source="gee" is the Earth Engine image collection, which adds the 2015-2018 epochs and the other LC100 cover fractions (esd.auth.login("earthengine")).

The first figure is the map. The second puts tree cover against elevation from the Copernicus DEM, also credential-free, and finds the treeline: the elevation band where the canopy fraction falls to zero.

import matplotlib.pyplot as plt
import numpy as np

import easysnowdata as esd

aoi = (-121.94, 46.72, -121.54, 46.99)  # Mount Rainier, WA

box = esd.parse_aoi(aoi)
grid = box.to_geobox(crs="utm", resolution=100)  # the map below is drawn on this grid
print(box)
print(grid)

AOI(bounds=(-121.9400, 46.7200, -121.5400, 46.9900), clip=True)
GeoBox((306, 310), Affine(100.0, 0.0, 580600.0,
       0.0, -100.0, 5205100.0), CRS('EPSG:32610'))

Where the product comes from: the Zenodo GeoTIFF needs nothing, the Earth Engine collection needs an account.

for src in esd.catalog.get("forest-cover-fraction").sources:
    print(f"{src.id:22} {src.title:45} {', '.join(src.requires) or 'no account'}")

zenodo                 Zenodo GeoTIFF                                no account
gee                    Earth Engine                                  earthengine

The 2019 tree-cover fraction. The 255 sentinel the product uses outside its footprint is masked to NaN and kept as the encoded nodata. The product is geographic; the box is loaded with a 1 km margin and drawn on the AOI’s UTM grid, bilinear because a fraction is continuous. The elevation binning below uses the native pixels.

forest_da = esd.land.forest_cover.load(box.buffer(1000), chunks=None)
print(forest_da)
print(
    f"pixels without a value (255 in the source): {float(forest_da.isnull().mean()):.1%}"
)
esd.plotting.map(
    forest_da.odc.reproject(grid, resampling="bilinear"),
    cmap="Greens",
    vmin=0,
    vmax=100,
    title="CGLS-LC100 tree cover fraction, 2019",
)

<xarray.DataArray 'tree_cover_fraction' (latitude: 291, longitude: 430)> Size: 501kB
array([[ 94.,  94.,  94., ...,  92.,  92., 100.],
       [100., 100.,  94., ...,  93.,  93., 100.],
       [ 96.,  96., 100., ...,  98.,  98., 100.],
       ...,
       [ 93.,  89.,  89., ..., 100.,  95., 100.],
       [ 86.,  95.,  95., ...,  96.,  96., 100.],
       [ 79.,  95.,  95., ...,  90.,  90.,  96.]],
      shape=(291, 430), dtype=float32)
Coordinates:
  * latitude     (latitude) float64 2kB 47.0 47.0 47.0 ... 46.71 46.71 46.71
  * longitude    (longitude) float64 3kB -122.0 -122.0 -122.0 ... -121.5 -121.5
    spatial_ref  int64 8B 0
    time         datetime64[us] 8B 2019-01-01
Attributes: (12/43)
    archive_facility:       VITO NV
    copyright:              Copernicus Service information 2020
    creator:                Dr. Marcel Buchhorn (VITO)
    delivered_product_crs:  WGS84 (EPSG:4326)
    doi:                    10.5281/zenodo.3939050
    file_creation:          Mon Jul 13 19:31:31 2020
    ...                     ...
    product_id:             forest-cover-fraction
    data_citation:          Buchhorn, M., et al. (2020). Copernicus Global La...
    license:                CC BY 4.0
    easysnowdata_version:   0.3.9.dev1+gd97689c2f
    units:                  %
    epoch:                  2019
pixels without a value (255 in the source): 0.0%

<Axes: title={'center': 'CGLS-LC100 tree cover fraction, 2019'}, xlabel='longitude [°]', ylabel='latitude [°]'>

Tree cover against elevation. The Copernicus DEM is 1 arc-second on the same geographic grid family, so it is sampled onto the 100 m forest grid with a bilinear interpolation; the pixels are then binned by 100 m of elevation.

dem_da = esd.terrain.dem.load(aoi, chunks=None).interp_like(forest_da, method="linear")
edges = np.arange(500, 3100, 100)
centres = edges[:-1] + 50
elevation = dem_da.values.ravel()
cover = forest_da.values.ravel()
valid = np.isfinite(elevation) & np.isfinite(cover)
which = np.digitize(elevation[valid], edges) - 1
median, low, high = (np.full(len(centres), np.nan) for _ in range(3))
for i in range(len(centres)):
    values = cover[valid][which == i]
    if len(values) >= 20:
        low[i], median[i], high[i] = np.percentile(values, [25, 50, 75])

bare = centres[np.nan_to_num(median, nan=100) < 5]
treeline = float(bare.min())
print(f"lowest 100 m band with median tree cover below 5 %: {treeline:.0f} m")

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.fill_betweenx(centres, low, high, color="#74c476", alpha=0.4, label="25-75 %")
ax.plot(median, centres, color="#006d2c", linewidth=2, label="median")
ax.axhline(treeline, color="0.3", linestyle="--", linewidth=1)
ax.text(100, treeline, f" treeline ≈ {treeline:.0f} m", va="bottom", ha="right")
ax.set_xlabel(esd.plotting.label(forest_da))
ax.set_ylabel(esd.plotting.label(dem_da, name="elevation"))
ax.set_xlim(0, 100)
ax.set_ylim(edges[0], float(np.nanmax(centres[np.isfinite(median)])) + 150)
ax.set_title("Tree cover by elevation band")
ax.legend(frameon=False, loc="upper right")
for side in ("top", "right"):
    ax.spines[side].set_visible(False)
fig.tight_layout()

lowest 100 m band with median tree cover below 5 %: 1950 m

The Earth Engine route, for the earlier epochs:

epochs_da = esd.land.forest_cover.load(aoi, source="gee", time="2015/2019")
epochs_da.mean(dim=("y", "x")).plot()