R & Chemometrics

Ejercicio 1: Utilidad de resta de matrices en espectroscopia

2026-05-27T21:14:33.809+02:00

Para seguir este espectro se dispone de la matriz de espectros "fish_1_19" que es una matriz de espectros que preparé en el post anterior y vamos a calcular el espectro medio, simplemente haciendo la media de cada columna (absorbancias a las 19 longitudes de onda) con Excel.

Extendemos la fórmula a las demás longitudes de onda

Y a ese espectro medio lo vamos a llamar el vector "Raw Mean Spectrum" para usarlo a continuación.

Y a continuación hacemos la resta de la matriz "Fish_1_19" de "Raw_Mean_Spectrum".

Simulando espectros de un NIR de 19 filtros en Excel

2026-05-27T16:56:09.184+02:00

A partir de un equipo NIR de rango 400 a 2498 cada dos nanómetros hemos seleccionado las 19 longitudes de onda correspondientes a los 19 filtros de un equipo Bran + Lubbe InfraAlyzer (ver anterior post). Por supuesto que es una gran perdida de datos, pero la idea de los siguientes posts es hacer análisis matriciales más simples para entender conceptos quimiométricos en Excel y de este modo será mucho mas pedagógico.

En este caso usamos un gráfico de scatter por líneas, pero si visualmente lo deseamos podemos quitarlos y poner solo los puntos.

Planificando hacer una base de datos espectral para análisis matricial en Excel.

2026-05-26T22:42:02.156+02:00

En los próximos posts vamos a preparar una matriz de datos espectral con un fichero de demo en el que seleccionaremos 19 longitudes de onda determinadas, simulando ser un equipo NIR de filtros Bran+Luebbe .

Las longitudes de onda serán: 1445, 1680, 1722, 1734, 1759, 1778, 1818, 1940, 1982, 2100, 2139, 2180, 2190, 2208, 2230, 2270, 2310, 2336 y 2348 nm.

1445 (O-H / Agua))

1680-1780 (C-H) 1720 - 1760 (Grasas / Lípidos)

1940 (O-H / Agua)

1980-2180 (N-H / Proteína)

2100 (Almidón / Proteína)

2170 - 2180 (Proteína)

2270-2348 (C-H Combinaciones)

Con estas selecciones trabajaremos con matrices mas pequeñas y menos intercorrelacionadas para poder entender mejor el análisis quimiométrico en el caso de trabajar con Excel.

Conocimientos básicos: Resta de matrices en Excel

2026-05-26T14:57:11.743+02:00

De gran utilidad es la resta de espectros en quimiometría NIR, por ejemplo, para conocer como ha variado una muestra en el tiempo (restando el espectro actual del original, nos indicará si la muestra a variado en humedad, color,..etc.), o por ejemplo para hacer un centrado espectral que es de lo más común en cálculos de algoritmos como por ejemplo Componentes Principales.

La resta matricial es muy simple en Excel, y usamos a modo de ejemplo las matrices A y B del post anterior:

siendo el resultado:

Conocimientos básicos: Suma de matrices en Excel

2026-05-25T20:28:49.596+02:00

Siguiendo con los temas de matrices en Excel, vamos a realizar en este caso la suma de matrices. Para ello usamos la matriz A del post anterior y vamos a sumarla con otra matriz que debe de cumplir la condición de que tiene el mismo numero de columnas y el mismo número de filas. Si quisiésemos llevar esto al análisis espectral sería como hacer una suma de unos espectros con otros.

Ahora lo que hacemos es definir la nueva matriz como B, seleccionar el rango de columnas y filas adecuado y aplicar la formula A + B, siendo A y B las matrices a sumar:

Al presionar "enter" nos dará el resultado de la suma:

Conocimientos básicos: Definiendo matrices en Excel

2026-05-25T17:55:30.586+02:00

Simplemente rellena un rango de celdas en Excel, si imaginamos que se trata de datos de espectros rellenamos datos un rango k.n, siendo k el numero de puntos de datos y n el número de muestras. En Excel podemos ir a fórmulas y en el administrados de nombre (Name Manager), lo nombramos con una letra mayúscula, por ejemplo A.

Usemos por ejemplo los datos que aparecen en el libro "Chemometrics in Excel":

Para representarlos gráficamente simplemente seleccionamos la matriz y en Insertar, seleccionamos los gráficos en 2D.

Como veis hemos representado una matriz A y con ella hemos visualizado los espectros. Imaginaros que se trata de un equipo con tres filtros de tres determinadas longitudes de onda y por los que hemos analizado tres muestras. Mas adelante lo haremos con espectros reales de equipos con un mayor numero de longitudes de onda.

Podemos hacer ahora una operación sencilla con esa matriz (multiplicarla por 3). Para ello seleccionamos un rango de celdas de 3.4 y en la primera celda escribimos =3*A y al presionar "enter" nos sale la nueva matriz.

Entender bien como funciona el análisis matricial es importante a la hora de comprender los fundamentos de la quimiometría.

Descarga y trabaja con datos de suelos en NIR

2026-05-25T15:12:18.075+02:00

Una de las páginas donde te puedes descargar datos de espectros de suelos con sus métodos de referencia es en

ZENODO: Agroscope_SoilSpectralLibrary_2024.

You can download soil spectra analysis to work with from:

ZENODO: Agroscope_SoilSpectralLibrary_2024

Git and Github for R Users

2025-10-06T19:24:00.003+02:00

Introduction to machine learning with {tidymodels}

2025-10-02T21:04:00.004+02:00

How to create the project with all the necessary data to follow the workshop:

1 - In RStudio select: File - New Project - Version Control - Git

2 - In the repository URL paste: https://github.com/nrennie/nhs-r-tidymodels.git

3 - Choose the folder and Create the project

Get started with tidymodels and classification of penguin data

2025-10-02T20:38:00.001+02:00

Guía rápida para el uso de tidymodels

2025-09-27T20:15:00.001+02:00

Interesante presentación en el grupo de usuarios de R de Madrid para el uso de modelos "tidy" a cargo de Jesús Herranz Valera.

Podéis descargaros las presentación así como el código que se muestra en la misma desde github.

Iremos comentando el uso de este tipo de modelos en próximos posts.

COMO INTERPRETAR ANÁLISIS DE SUELOS PARTE 1

2025-06-27T17:22:00.005+02:00

Digital Soil Spectroscopy

2025-06-19T10:52:00.005+02:00

Determinación de fósforo en suelos (Método Olsen)

2025-06-19T10:44:00.000+02:00

Standard Normal Variate (SNV)

2025-06-17T17:32:00.006+02:00

As in previous posts, we are going to follow the Soil spectroscopy training material and we are going to use the reference material provided.

Let´s start loading the libraries and the data:

library(tidyverse)
library(prospectr)
load("C:/BLOG/Workspaces/NIR Soil Tutorial/post5.RData")

In the training material we can read a clear explanation of what SNV does:

“Standard normal variate is a simple and widely used method that works by centering each individual spectrum to zero and then dividing each spectral band value by the standard deviation of the entire spectrum. The SNV method processes each observation independently. The disadvantage of SNV is that it can be sensitive to noise”.

So for each spectrum we calculate the mean of all data points reflectance values and also their standard deviation. Then we subtract for every data point the mean and divide the result by the standard deviation.

Using the prospectr package, we can use the standardNormalVariate function to apply the SNV transformation to the spectra matrix.

snv_spectra <- standardNormalVariate(dat$spc)

Let´s use first classical R function matplot to plot the spectra before and after the SNV transformation:

matplot(colnames(dat$spc),
  t(snv_spectra),
  type = "l",
  col = "grey",
  lwd = 0.5,
  xlab = "Wavelength (nm)",
  ylab = "Reflectance",
  main = "Spectra after SNV"
)

Now with ggplot2:

# Create the long-format data frame for the different derivatives
snv_spectra_long <- data.frame(
sample = rep(1:nrow(dat), each = ncol(snv_spectra)),
oc = rep(dat$Organic_Carbon, each = ncol(snv_spectra)),
clay = rep(dat$Clay, each = ncol(snv_spectra)),
silt = rep(dat$Silt, each = ncol(snv_spectra)),
sand = rep(dat$Sand, each = ncol(snv_spectra)),
wavelength = rep(my_wavelengths, nrow(snv_spectra)),
absorbance = as.vector(t(snv_spectra))
)


ggplot(snv_spectra_long, aes(x = wavelength, y = absorbance, 
                             group = sample, color = oc)) +
  geom_line(alpha = 0.5) + # Set alpha to 0.5 for transparency
  scale_color_gradient(low = 'blue', high = 'red') +
  theme_minimal() +
  labs(x = 'Wavelength (nm)', y = 'Reflectance', color = 'Organic Carbon (%)')

Bibliography:

Soil spectroscopy training material Wadoux, A., Ramirez-Lopez, L., Ge, Y., Barra, I. & Peng, Y. 2025. A course on applied data analytics for soil analysis with infrared spectroscopy – Soil spectroscopy training manual 2. Rome, FAO.

NIR Chemometrics for soil analysis

2025-06-13T17:04:00.007+02:00

Adding noise to the spectra

2025-06-11T16:15:00.005+02:00

We have seen in the previous post how to remove artifacts in the spectra using the spliceCorrection function from the prospectr package. But there are other sources of noise in the spectra which depends of several factors, as the quality of the spectrophotometer (can use a signal to noise ratio meassurement to test it), the scanning time, resolution, sample presentation, …. In this post we continue with the “course on applied data analytics for soil analysis with infrared spectroscopy” using the reference material provided. And this time we are going to see the part of the course where we add noise to the spectra and then we apply two methods to remove it. The methods are the Moving Window Average and the Savitzky-Golay Filter.

First we load the libraries we are going to use:

library(tidyverse)
load("C:/BLOG/Workspaces/NIR Soil Tutorial/post4.RData")

Adding noise to the first spectrum

We create some Gaussian noise to the first spectrum in the dataset (mean = 0, standard deviation = 0.0025), and we are gong to add it to the first spectrum in the dataset. Take in account that we are going to select the dataset dat$spc, which is the one we created in the previous post, so we are going to use the spectra without the artifacts.

#Set a seed for reproducibility of random noise
set.seed(801124)
# Add random noise to the first sample’s spectrum
my_noisy_spc <- dat$spc[1, ] + rnorm(ncol(dat$spc), sd = 0.0025)

If we just want to se the noise that we are adding:

ggplot(data = data.frame(wavelength = my_wavelengths, absorbance = rnorm(ncol(dat$spc), sd = 0.0025)),
       aes(x = wavelength, y = absorbance)) +
  geom_line(color = "blue") +
  labs(title = "Noisy Spectrum",
       x = "Wavelength (nm)",
       y = "Absorbance") +
  theme_minimal()

Plotting the noisy spectrum

ggplot(data = data.frame(wavelength = my_wavelengths, absorbance = my_noisy_spc),
       aes(x = wavelength, y = absorbance)) +
  geom_line(color = "blue") +
  labs(title = "Spectrum 1 with added Noise",
       x = "Wavelength (nm)",
       y = "Absorbance") +
  theme_minimal()

Removing the noise

No let´s apply to this spectrum two methods used to remove noise:

Moving Window Average

library(prospectr)
# Replace ‘my_noisy_spc’ with your actual data frame or matrix and ‘w’ with the desired
# window size
mwa_result <- movav(my_noisy_spc, w = 11)

Savitzky-Golay Filter

# Replace ‘my_noisy_spc’ with your actual data frame or matrix, ‘m’ with the window size,
# and ‘p’ with the polynomial order
sg_result <- savitzkyGolay(
my_noisy_spc, # the noisy spectrum
m = 0, # this m is for the derivative order, 0 is for no derivative
p = 3, # the polynomial order
w = 11 # the window size
)

# Create a long-format data frame for the MVA spectrum
long_mwa_spc <- data.frame(
wavelength = as.numeric(names(mwa_result)), # Extract the wavelengths
spectral_value = mwa_result, # Store the spectral values from the MVA
method = 'mva' # Label the data
)
# Create a long-format data frame for the Savitzky-Golay (SG) spectrum
long_sg_spc <- data.frame(
wavelength = as.numeric(names(sg_result)), # Extract the wavelengths
spectral_value = sg_result, # Store the spectral values from the SG
method = 'sg' # Label the data
)
# Combine the long-format MVA and SG data frames into one
denoised_long <- rbind(long_mwa_spc, long_sg_spc)

Now we create a data frame for the original noisy spectrum and some labels.

# Create a data frame for the original noisy spectrum
original_noisy <- data.frame(
wavelength = my_wavelengths, # Use the original wavelengths
spectral_value = my_noisy_spc # Use the original noisy spectrum values
)

my_labels <- c(
'mva' = 'Moving average smoothing',
'sg' = 'Savitzky-Golay filtering'
)

Now we want to overplot the corrected spectra over the noisy one:

# Plot the spectra using ggplot2
ggplot() +
# Plot the original noisy spectrum as a red line in the background
geom_line(
data = original_noisy,
aes(x = wavelength, y = spectral_value),
color = 'black'
) +
# Plot the MVA and SG processed spectra on top of the original noisy spectrum
geom_line(
data = denoised_long,
aes(x = wavelength, y = spectral_value, color = method, group = method)
) +
# Add titles and labels for the axes
labs(
x = 'Wavelength (nm)',
y = 'Reflectance'
) +
#coord_cartesian(xlim = c(1500, 1800), ylim = c(0.425, 0.450)) +
# Manually set the colors for the MVA and SG lines
scale_color_manual(
values = c('mva' = '#00AFBB', 'sg' = '#E7B800'),
labels = my_labels
) +
# Facet the plot by the method, creating separate plots for MVA and SG
facet_grid(. ~ method, labeller = labeller(method = my_labels)) +
theme(legend.position = 'top') + guides(color = guide_legend(title = NULL))

We can zoom the spectra to see the differences more clearly:

ggplot() +
# Plot the original noisy spectrum as a red line in the background
geom_line(
data = original_noisy,
aes(x = wavelength, y = spectral_value),
color = 'black'
) +
# Plot the MVA and SG processed spectra on top of the original noisy spectrum
geom_line(
data = denoised_long,
aes(x = wavelength, y = spectral_value, color = method, group = method), size = 1
) +
# Add titles and labels for the axes
labs(
x = 'Wavelength (nm)',
y = 'Reflectance'
) +
coord_cartesian(xlim = c(1500, 1800), ylim = c(0.425, 0.450)) +
# Manually set the colors for the MVA and SG lines
scale_color_manual(
values = c('mva' = '#00AFBB', 'sg' = '#E7B800'),
labels = my_labels
) +
# Facet the plot by the method, creating separate plots for MVA and SG
facet_grid(. ~ method, labeller = labeller(method = my_labels)) +
theme(legend.position = 'top') + guides(color = guide_legend(title = NULL))

Bibliography:

Removing the artifacts in the spectra

2025-06-09T12:08:00.000+02:00

To follow this post, remember than in the first post we build a matrix with the raw NIR spectra, and that matrix is the one we are going to treat in this post to remove the spectra so we will give to the new matrix (without the artifacts). This will be the first of several conversions of the spsectra matris where we continue applying other math treatments.

First we load the libraries we are going to use and the data we prepared in the previous post:

library(tidyverse)
library(prospectr)
load("C:/BLOG/Workspaces/NIR Soil Tutorial/post2.RData")

The function we are going to use is spliceCorrection from the prospectr package. This function is used to correct the spectra for artifacts due to detector changes or other reasons. The function takes as input a matrix of spectra and returns a corrected matrix of spectra. To find information:

?prospectr::spliceCorrection

dat$spc <- spliceCorrection(dat$spc_raw, my_wavelengths, splice = c(1000, 1830))

And now we create a long version of the dataframe to plot it with ggplot2, as we have done in the second post:

spc_long <- data.frame(
sample = rep(1:nrow(dat), each = ncol(dat$spc)),
oc = rep(dat$Organic_Carbon, each = ncol(dat$spc)),
clay = rep(dat$Clay, each = ncol(dat$spc)),
silt = rep(dat$Silt, each = ncol(dat$spc)),
sand = rep(dat$Sand, each = ncol(dat$spc)),
wavelength = rep(my_wavelengths, nrow(dat)),
absorbance = as.vector(t(log(1 / dat$spc, 10))))


ggplot(
spc_long,
aes(x = wavelength, y = absorbance, group = sample, color = oc)) +
geom_line(alpha = 0.5) + # Set alpha to 0.5 for transparency
scale_color_gradient(low = 'blue', high = 'red') +
theme_minimal() +
labs(
x = 'Wavelength (nm)',
y = 'Reflectance',
color = 'Organic Carbon (%)') +
theme(legend.position = 'top')

Now we can check if the artifacts are removed:

ggplot(spc_long,
       aes(x = wavelength, y = absorbance, group = sample)) +
  geom_line(alpha = 0.5) +
  theme_minimal() +
  labs(
    title = "Artifacts at 1000 nm removed",
    x = 'Wavelength (nm)',
    y = 'Reflectance'
  ) +
  coord_cartesian(xlim = c(990, 1010), ylim = c(0.4, 0.45)) +
  geom_vline(xintercept = 1000, linetype = "dashed", color = "red")

Looking for artifacts in the spectra

2025-05-31T20:21:00.006+02:00

In this post we are going to look to the NIR spectra in detail looking for artifacts in the sample:

Plotting the spectra :

Looking to the NIR spectra we can see that there are some artifacts in the spectra at 1000 nm and 1830 nm. These artifacts are due to the change of detector in the spectrometer. We can plot these artifacts using the ggplot2 package.

In this case we are going to vave the plots and use the patchwork package to combine them.

library(tidyverse)
library(patchwork)
load("C:/BLOG/Workspaces/NIR Soil Tutorial/post2.RData")

artifact_1000nm <- ggplot(vnir_long,
       aes(x = wavelength, y = absorbance, group = sample)) +
  geom_line(alpha = 0.5) +
  theme_minimal() +
  labs(
    title = "Artifacts at 1000 nm",
    subtitle = "Due to detector change",
    x = 'Wavelength (nm)',
    y = 'Calculated absorbance'
  ) +
  coord_cartesian(xlim = c(990, 1010), ylim = c(0.4, 0.45)) +
  geom_vline(xintercept = 1000, linetype = "dashed", color = "red")

artifact_1830nm <- ggplot(vnir_long,
       aes(x = wavelength, y = absorbance, group = sample)) +
  geom_line(alpha = 0.5) +
  theme_minimal() +
  labs(
    title = "Artifacts at 1830 nm",
    subtitle = "Less visible than at 1000 nm",
    x = 'Wavelength (nm)',
    y = 'Calculated absorbance'
  ) +
  coord_cartesian(xlim = c(1820, 1840), ylim = c(0.2, 0.4)) +
  geom_vline(xintercept = 1830, linetype = "dashed", color = "red")

Combining the plots

We can combine both plots:

artifact_1000nm + artifact_1830nm + plot_layout(ncol = 2)

Bibliography:

Follow the posts on Netlify

Plotting MIR spectra with ggplot2

2025-05-31T20:18:00.003+02:00

Loading and preparing the MIR data.

The data is available in the data folder of the repository. Is in a CSV format and can be loaded using the read_csv function from the readr package.

library(readr)

url_mir <- "https://raw.githubusercontent.com/FAO-SID/SoilFER-Spec/main/data/dat2MIR.csv"
dat2MIR <- read_csv(url_mir)

Now we will follow the instructions of the Soil spectroscopy training material changing the name of the first column (index of the samples) to “sample”.

colnames(dat2MIR)[1] <- "sample"
head(dat2MIR, c(10, 7))

# A tibble: 10 × 7
   sample `4001.65608` `3999.72758` `3997.79907` `3995.87056` `3993.94205`
    <dbl>        <dbl>        <dbl>        <dbl>        <dbl>        <dbl>
 1      1       0.180        0.18         0.180        0.180        0.181 
 2      2       0.150        0.151        0.151        0.151        0.151 
 3      3       0.0484       0.0487       0.0489       0.0491       0.0492
 4      4       0.133        0.133        0.133        0.134        0.134 
 5      5       0.227        0.228        0.228        0.228        0.228 
 6      6       0.120        0.120        0.121        0.121        0.121 
 7      7       0.312        0.313        0.313        0.313        0.314 
 8      8       0.225        0.226        0.226        0.226        0.227 
 9      9       0.218        0.219        0.220        0.220        0.221 
10     10       0.113        0.113        0.114        0.114        0.114 
# ℹ 1 more variable: `3992.01354` <dbl>

Looking to the CSV file, we see that the first column is the sample number, and the last one the organic carbon parameter. All the columns in the middle ate the wavenumbers.

Let´s create a vector with the wavenumbers. The wavenumbers are in the column names from the second to the penultimate column, so we can extract them using the following code:

wavenumbers <- as.numeric(colnames(dat2MIR)[2:(ncol(dat2MIR) - 1)])

Let´s prepare the wavenumbers vector to be used as column names in the matrix of spectra.

wavelengths_ir <- round(10000000 / wavenumbers)
wavenumbers_ir <- round(10000000/wavelengths_ir) # Convert to cm-1

my_spectra_ir <- as.matrix(dat2MIR[, 2:(ncol(dat2MIR) - 1)])
colnames(my_spectra_ir) <- wavenumbers_ir

Let´s prepare the MIR dataframe to work with it:

dat_mir <- dat2MIR[, c(1, 1767)]
dat_mir$spc_raw_ir <- my_spectra_ir
rm(my_spectra_ir)

library(tidyverse)

my_wavenumbers <- as.numeric(colnames(dat_mir$spc_raw_ir))

#creating the long dataframe
ir_long <- data.frame(
sample = rep(1:nrow(dat_mir), each = ncol(dat_mir$spc_raw_ir)),
oc = rep(dat_mir$Organic_Carbon, each = ncol(dat_mir$spc_raw_ir)),
wavenumber = rep(my_wavenumbers, nrow(dat_mir)),
absorbance = as.vector(t(dat_mir$spc_raw_ir))
)

head(ir_long)

  sample        oc wavenumber absorbance
1      1 0.7950426       4002    0.17972
2      1 0.7950426       4000    0.18000
3      1 0.7950426       3998    0.18025
4      1 0.7950426       3995    0.18048
5      1 0.7950426       3994    0.18076
6      1 0.7950426       3992    0.18107

We can use the ggplot function to create a plot.

ggplot(ir_long,
       aes(x = wavenumber, y = absorbance, group = sample, color = oc)) +
  geom_line(alpha = 0.5) +
  scale_color_gradient(low = 'blue', high = 'red') +
  scale_x_reverse() +  # Invierte el eje X
  theme_minimal() +
  labs(x = 'Wavenumber (cm⁻¹)',
       y = 'Calculated absorbance',
       color = 'OC (%)')

Bibliography:

Follow the NIR-Chemometrics blog on Netlify

Plotting NIR spectra with ggplot2

2025-05-31T20:04:00.001+02:00

This is the post number two of the Soil spectroscopy training material. We are following the paper Soil spectroscopy training material and at the same time we make small changes to the code.

In the previous post, we plot the NIR spectra using the classical R function matplot. Now we want to use the ggplot2 package to plot the spectra but for that we have to prepare the data in a different way, creating a long format indeed a wide format.

library(tidyverse)
load("dat.RData")
load("C:/BLOG/Workspaces/NIR Soil Tutorial/post1.RData")
#converting the wavelengths into numeric values
my_wavelengths <- as.numeric(colnames(dat$spc_raw))

#creating the long dataframe
vnir_long <- data.frame(
sample = rep(1:nrow(dat), each = ncol(dat$spc_raw)),
oc = rep(dat$Organic_Carbon, each = ncol(dat$spc_raw)),
clay = rep(dat$Clay, each = ncol(dat$spc_raw)),
silt = rep(dat$Silt, each = ncol(dat$spc_raw)),
sand = rep(dat$Sand, each = ncol(dat$spc_raw)),
wavelength = rep(my_wavelengths, nrow(dat)),
absorbance = as.vector(t(log(1 / dat$spc_raw, 10))))

head(vnir_long)

  sample        oc     clay silt sand wavelength absorbance
1      1 0.7950426 14.95713 40.1 44.9        350   1.206405
2      1 0.7950426 14.95713 40.1 44.9        351   1.204785
3      1 0.7950426 14.95713 40.1 44.9        352   1.225360
4      1 0.7950426 14.95713 40.1 44.9        353   1.239876
5      1 0.7950426 14.95713 40.1 44.9        354   1.238837
6      1 0.7950426 14.95713 40.1 44.9        355   1.235567

We can use the ggplot function to create a plot.

ggplot(vnir_long,
aes(x = wavelength, y = absorbance, group = sample, color = clay)) +
geom_line(alpha = 0.5) + # Set alpha to 0.5 for transparency
scale_color_gradient(low = 'blue', high = 'red') +
theme_minimal() +
  labs(x = 'Wavelength (nm)',y = 'Calculated absorbance',color = 'Clay (%)')

This plot shows the spectra of the samples with the clay content as a color gradient (the paper use the organic carbon). The x-axis represents the wavelength in nano-meters, and the y-axis represents the calculated absorbance. The color gradient indicates the clay content, with blue representing the samples with the lower clay content and the red ones with the higher clay content.

Due to the scatter we don´t see any patters clearly in the spectra yet.

The paper also show us how to work with the MIR (Middle Infrared Spectra), so in the next post we will do the same as we did in these first two posts, but with the MIR dataframe.

Bibliography:

New "Soil Near Infrared" Training Material (part 1)

2025-04-25T20:42:00.001+02:00

Load and prepare the data.

The data is available in the data folder of the repository. Is in a CSV format and can be loaded using the read_csv function from the readr package.

library(readr)

url <- "https://raw.githubusercontent.com/FAO-SID/SoilFER-Spec/main/data/dat2.csv"
dat2 <- read_csv(url)

Now we will follow the instructions of the Soil spectroscopy training material changing the name of the first column (index of the samples) to “sample”.

colnames(dat2)[1] <- "sample"
head(dat2, c(10, 7))

# A tibble: 10 × 7
   sample Organic_Carbon  Clay  Silt  Sand   X350   X351
    <dbl>          <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
 1      1          0.795 15.0   40.1  44.9 0.0622 0.0624
 2      2          0.696  8.79  52.5  38.7 0.0735 0.0725
 3      3          1.46  13.8   41.8  44.4 0.0685 0.0699
 4      4          3.36  31.1   49.5  19.4 0.0590 0.0615
 5      5          3.71  33.1   62.3   4.6 0.0624 0.0626
 6      6          1.39  33.3   60.6   6.1 0.0997 0.0988
 7      7          3.38  53.8   44.7   1.5 0.0920 0.0936
 8      8          3.16  27.3   66.3   6.4 0.0703 0.0673
 9      9          3.49  30.5   48.5  21   0.0473 0.0440
10     10          2.34  21.5   61.5  17   0.0666 0.0691

Now we have the sample number, the four parameters and the wavelengths from 350 to 2500 nm in steps of 1 nm. In total we have 2151 wavelengths also called datapoints. Our next step is to remove the X from the wavelengths column names. There are several ways to do this. One of the most common ways is to use the gsub function which replaces all occurrences of a pattern in a string with a replacement string. In this case, we want to replace “X” with ““ (an empty string). We can do this for all the columns except the first four (sample number and parameters) using the following code:

colnames(dat2)[-c(1:5)] <- gsub("X", "", colnames(dat2)[-c(1:5)])
head(dat2, c(10, 7))

# A tibble: 10 × 7
   sample Organic_Carbon  Clay  Silt  Sand  `350`  `351`
    <dbl>          <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
 1      1          0.795 15.0   40.1  44.9 0.0622 0.0624
 2      2          0.696  8.79  52.5  38.7 0.0735 0.0725
 3      3          1.46  13.8   41.8  44.4 0.0685 0.0699
 4      4          3.36  31.1   49.5  19.4 0.0590 0.0615
 5      5          3.71  33.1   62.3   4.6 0.0624 0.0626
 6      6          1.39  33.3   60.6   6.1 0.0997 0.0988
 7      7          3.38  53.8   44.7   1.5 0.0920 0.0936
 8      8          3.16  27.3   66.3   6.4 0.0703 0.0673
 9      9          3.49  30.5   48.5  21   0.0473 0.0440
10     10          2.34  21.5   61.5  17   0.0666 0.0691

Now we can isolate the wavelengths into a matrix. We can do this by using the as.matrix function. The as.matrix function converts a data frame to a matrix. In this case, we want to convert all columns except the first four columns (sample number and parameters) to a matrix. We can do this using the following code:

my_spectra <- as.matrix(dat2[, -c(1:5)])

Now we remove the index and parameter columns from dat2:

dat <- dat2[, c(1:5)]

and reassign the spectra to dat as a single variable called spc_raw, and remove my_spectra from the environment.

dat$spc_raw <- my_spectra
rm(my_spectra)

As a first check we can use matplot to visualize the spectra. We can use the matplot function to plot the spectra. The matplot function is used to plot matrices. In this case, we want to plot the spectra (the matrix) against the wavelengths (the column names of the matrix). We can do this using the following code:

Plot the spectra (using classical R)

matplot(colnames(dat$spc_raw), t(dat$spc_raw), type = "l", lty = 1, col = "grey", xlab = "Wavelength (nm)", ylab = "Reflectance (%)", main = "NIR Spectra of the samples")

This is my first approach to this training material. As you see in the paper the authors use this time the ggplot2 package to plot the spectra, and I will use it in the next post.

save.image("C:/BLOG/Workspaces/NIR Soil Tutorial/post1.RData")

Bibliography:

New Soil "Near Infrared Spectroscopy" Training Material

2025-04-13T00:00:00.001+02:00

Load and prepare the data.

The data is available in the data folder of the repository. Is in a CSV format and can be loaded using the read_csv function from the readr package.

library(readr)

url <- "https://raw.githubusercontent.com/FAO-SID/SoilFER-Spec/main/data/dat2.csv"
dat2 <- read_csv(url)

Now we will follow the instructions of the Soil spectroscopy training material changing the name of the first column (index of the samples) to “sample”.

colnames(dat2)[1] <- "sample"
head(dat2, c(10, 7))

# A tibble: 10 × 7
   sample Organic_Carbon  Clay  Silt  Sand   X350   X351
    <dbl>          <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
 1      1          0.795 15.0   40.1  44.9 0.0622 0.0624
 2      2          0.696  8.79  52.5  38.7 0.0735 0.0725
 3      3          1.46  13.8   41.8  44.4 0.0685 0.0699
 4      4          3.36  31.1   49.5  19.4 0.0590 0.0615
 5      5          3.71  33.1   62.3   4.6 0.0624 0.0626
 6      6          1.39  33.3   60.6   6.1 0.0997 0.0988
 7      7          3.38  53.8   44.7   1.5 0.0920 0.0936
 8      8          3.16  27.3   66.3   6.4 0.0703 0.0673
 9      9          3.49  30.5   48.5  21   0.0473 0.0440
10     10          2.34  21.5   61.5  17   0.0666 0.0691

colnames(dat2)[-c(1:5)] <- gsub("X", "", colnames(dat2)[-c(1:5)])
head(dat2, c(10, 7))

# A tibble: 10 × 7
   sample Organic_Carbon  Clay  Silt  Sand  `350`  `351`
    <dbl>          <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
 1      1          0.795 15.0   40.1  44.9 0.0622 0.0624
 2      2          0.696  8.79  52.5  38.7 0.0735 0.0725
 3      3          1.46  13.8   41.8  44.4 0.0685 0.0699
 4      4          3.36  31.1   49.5  19.4 0.0590 0.0615
 5      5          3.71  33.1   62.3   4.6 0.0624 0.0626
 6      6          1.39  33.3   60.6   6.1 0.0997 0.0988
 7      7          3.38  53.8   44.7   1.5 0.0920 0.0936
 8      8          3.16  27.3   66.3   6.4 0.0703 0.0673
 9      9          3.49  30.5   48.5  21   0.0473 0.0440
10     10          2.34  21.5   61.5  17   0.0666 0.0691

my_spectra <- as.matrix(dat2[, -c(1:5)])

Now we remove the index and parameter columns from dat2:

dat <- dat2[, c(1:5)]

and reassign the spectra to dat as a single variable called spc_raw, and remove my_spectra from the environment.

dat$spc_raw <- my_spectra
rm(my_spectra)

Plot the spectra (using classical R)

matplot(colnames(dat$spc_raw), t(dat$spc_raw), type = "l", lty = 1, col = "grey", xlab = "Wavelength (nm)", ylab = "Reflectance (%)", main = "NIR Spectra of the samples")

This is my first approach to this training material. As you see in the paper the authors use this time the ggplot2 package to plot the spectra, and I will use it in the next post.

save.image("C:/BLOG/Workspaces/NIR Soil Tutorial/post1.RData")

Bibliography:

Estadísticas de tráfico marítimo en el puerto de Avilés (parte 2)

2025-03-19T12:05:00.002+01:00

Una de las preguntas que nos podemos hacer con los datos que proporciona la web del puerto de Avilés, es:

¿Que mercancías se exportaron principalmente del Puerto de Avilés en 2024?.

La respuesta es clara:

No es sorprendente, ya que en Avilés se encuentra una fábrica de ArcelorMittal.

podemos seguir haciéndonos preguntas:

Cual es el destino de los buques con estas boninas embarcadas?.

Vemos que el principal destino de las bobinas es Ravenna (Italia).

¿Qué tonelaje tenía cada envío a lo largo del 2024?

Podemos comprobar como las herramientas de los paquetes de R (ggplot2,...) nos permiten realizar gráficos muy claros para interpretar datos de hojas Excel o bases de datos que dispongamos.

Estadísticas de trafico marítimo en el puerto de Avilés (parte 1)

2025-03-18T12:58:00.000+01:00

Con el fin de coger práctica con R y las herramientas de tidyverse, durante los próximos post trabajaré con los datos que aporta la página web del puerto de Avilés para analizar diferentes estadísticas sobre el uso del puerto para embarque y desembarque de mercancías en el año 2024. Para ello es interesante acceder a la página web de dicho puerto y descargar la hoja Excel con dichos datos e importar dicha hoja a R, aunque es preferible que arregléis previamente dicha Excel antes de la importación a RStudio.

Este es el enlace de la página web del puerto: puerto de Avilés

En está página hay un enlace a "Buques en el puerto" y desde ahí podéis seleccionar la opción "Histórico de movimientos" y descargaros la hoja excel que incluye los datos desde el 1 de enero de 2024 hasta la fecha que hagáis la descarga.

He estado trabajando con la hoja Excel y pronto la tendré disponible si queréis ahorraros los pasos, aunque os los recomiendo para que vayáis cogiendo soltura con el manejo de Excel previo a la importación a RStudio.

Recientemente he estado haciendo unas visitas al puerto y es muy interesante el tema del tráfico marítimo, así como toda la industria que se genera entorno a los puertos.

Incluyo en este post una foto de un buque de más de 190 metros de eslora entrando al puerto.