1. Background

The measurement of carbon fluxes in forest ecosystems can be a challenging task, but there are alternative methods to assess net ecosystem exchange (NEE) of CO2. One of these methods involves the use of remote sensing techniques, specifically the use of Normalized Difference Vegetation Index (NDVI) and Photochemical Reflectance Index (PRI). By measuring changes in these indices over time, researchers can estimate carbon fluxes at daily, monthly, and yearly time scales.

NDVI is a commonly used index that measures the greenness of vegetation, indicating photosynthetic activity. On the other hand, PRI measures the absorption of solar radiation by the chlorophyll and carotenoid pigments in plants, indicating the efficiency of photosynthesis. By combining these two indices, researchers can estimate the rate of carbon uptake by plants, as well as the amount of carbon released through respiration.

These remote sensing techniques have been used globally to estimate Gross Primary Productivity (GPP) or Net Primary Productivity (NPP) at large spatial scales. However, accurate estimates of these parameters depend on proper parameterization of the models used, including maximum LUE (ε*) and the temperature optima (Topt) of gross photosynthesis. Eddy covariance (EC) technique, which directly measures NEE across the canopy-atmosphere interface, is being used to derive these ecosystem parameters as input to NDVI and PRI-based models.

Remote sensing techniques play a key role in the optimization of ecosystem model parameters, validation of its outputs, and up-scaling of CO2 flux across spatial or temporal scales. Although the use of NDVI and PRI-based models for carbon flux estimation is still in its early stages, it shows great potential for accurate and consistent estimates of carbon fluxes in forest ecosystems.

2. Objectives

How efficient is the Light Use Efficiency (LUE) method based on remote sensing in estimating carbon fluxes in tropical forests using Landsat8 images and the technique developed by Rahman et al. (2000)?

# Load packages
library(rgdal)
## Warning: package 'rgdal' was built under R version 4.2.3
## Carregando pacotes exigidos: sp
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/Casa/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/Casa/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(raster)
## Warning: package 'raster' was built under R version 4.2.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
# Define working directory
setwd("C:/Users/Casa/Desktop/LC08_L1TP_220076_20220901_20220912_02_T1")

# Load image bands as raster objects
blue <- raster("LC08_L1TP_220076_20220901_20220912_02_T1_B2.tif")
green <- raster("LC08_L1TP_220076_20220901_20220912_02_T1_B3.tif")
red <- raster("LC08_L1TP_220076_20220901_20220912_02_T1_B4.tif")
nir <- raster("LC08_L1TP_220076_20220901_20220912_02_T1_B5.tif")
swir1 <- raster("LC08_L1TP_220076_20220901_20220912_02_T1_B6.tif")

# Calculate NDVI
ndvi <- (nir - red) / (nir + red)

# Plot NDVI
plot(ndvi, main = "NDVI")

# Calculate PRI
pri <- (swir1 - red) / (swir1 + red)

# Plot PRI
plot(pri, main = "PRI")

# Calculate CO2Flux using Ryu et al. (2012) model
alpha <- 0.0019
beta <- -0.0696
gamma <- 0.1427
co2flux <- alpha * ndvi + beta * pri + gamma

# Plot CO2Flux
plot(co2flux, main = "CO2Flux")

# Hist distribution co2flux
hist(co2flux, breaks = 1000, col = "green")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

# Calculate minimum and maximum pixel values in co2flux raster
minval <- cellStats(co2flux, stat = "min", na.rm = TRUE)
maxval <- cellStats(co2flux, stat = "max", na.rm = TRUE)

# Define break points for the 5 classes
breaks <- seq(from = minval, to = maxval, length.out = 6)

# Classify pixels into 5 classes
classes <- cut(co2flux, breaks = breaks, include.lowest = TRUE)

# Convert classes to vector
classes <- as.vector(classes)

# Count the number of pixels in each class
table(classes)
## classes
##        1        2        3        4        5 
##       20 10177180 30824550   733082     1358
# Create data frame with NDVI and PRI values
df <- data.frame(NDVI = values(ndvi), PRI = values(pri))

# Histogram of NDVI values
ggplot(df, aes(x = NDVI)) +
  geom_histogram() +
  ggtitle("Histogram of NDVI Values") +
  xlab("NDVI Values") +
  ylab("Frequency")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 19507137 rows containing non-finite values (`stat_bin()`).

# Histogram of PRI values
ggplot(df, aes(x = PRI)) +
  geom_histogram() +
  ggtitle("Histogram of PRI Values") +
  xlab("PRI Values") +
  ylab("Frequency")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 19507895 rows containing non-finite values (`stat_bin()`).

# Histogram of CO2Flux values
library(ggplot2)
ggplot(data.frame(co2flux = as.vector(co2flux)), aes(x = co2flux)) +
  geom_histogram() +
  ggtitle("Histogram of CO2Flux Values") +
  xlab("CO2Flux Values") +
  ylab("Frequency")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 19508061 rows containing non-finite values (`stat_bin()`).

Results

In river regions, the CO2Flux tends to be higher due to the presence of decaying organic matter in the river water and sediment. The decomposition of organic matter releases CO2 into the water, which is then transferred to the atmosphere through the process of diffusion. In addition, biological activity in the water, such as the respiration of aquatic organisms, can also contribute to the increase in CO2Flux. That is, the river courses do not serve as a carbon sink, demonstrated for the Sorocaba region by the high values, which vary between 0.14 and 0.17 μmol CO2 m⁻² s⁻¹.

On the other hand, in forest regions, the CO2Flux tends to be lower due to the presence of vegetation that absorbs CO2 through the process of photosynthesis. During photosynthesis, plants convert carbon dioxide and water into sugars and oxygen, releasing oxygen into the atmosphere and absorbing CO2. In addition, forests often have soils rich in organic matter, which can result in high rates of carbon fixation in the soil. Thus, forest regions 0.125 to 0.135 μmol CO2 m⁻² s⁻¹.

Conclusions

The large water courses in the region do not act as a carbon sink. Therefore, this behavior can be explained by the high rate of purification of water courses by the kinetic energy of water and carbon exchange between water and atmosphere.