Spectral Unmixing Functionality in RStoolbox

Spectral Unmixing Functionality in RStoolbox

written by Martin Wegmann

March 20, 2018

Please find below a blog post by one of our EAGLE M.Sc. students, Jakob Schwalb-Willmann, about his implementation of spectral unmixing into the RStoolbox R package, developed at our department by Benjamin Leutner:

 

“Recently, in January, I finished the development of the first version of a spectral unmixing function being part of RStoolbox, an R package offering numerous tools for remote sensing analysis written by Benjamin Leutner. The multiple endmember spectral mixture analysis (mesma) function makes it possible to unmix multi- and hyper-spectral imagery by sets of spectral endmember profiles.

For this, a non-negative least squares (NNLS) solver was implemented. NNLS is a statistical approach to fit model parameters to data, assuming that the model parameters are always expressed linearly to those not expressed by the model (unknown parameters) and that the model parameters can never be negative. There are different approaches to solve the NNLS problem. A popular one had been introduced by Lawson & Hanson (1974), which can be considered as fundamental work on practical Least Square Problems solving. It was originally published as FORTRAN code, which is still widely used, e.g. by the R NNLS package or by the Python scipy library. However, compared to newer developing frameworks and languages emerged, the FORTRAN implementation is relatively slow. Thus and to be independent from existing solutions, I wrote a C++ NNLS solver for mesma(), based on a sequential coordinate-wise algorithm (SCA) introduced by Franc et al. (2005). The latter method inherits strong control about the solver’s iteration stopping conditions.

Apart from NNLS, we are planning to add further solver methods (which is why I am currently looking for other practicable unmixing methods).

In this post, I want to demonstrate, how to use the unmixing function RStoolbox::mesma() by giving an example. It can be reproduced on any device running R and having installed the current RStoolbox beta. To do so, execute devtools::install_github(“bleutner/RStoolbox”).

First, we are going to load the Landsat example imagery delivered with RStoolbox:

1
2
3
4
5
6
7
8
9
10
#if not already done, install RStoolbox beta:
devtools::install_github("bleutner/RStoolbox")
#required caret 6.0-79, which is not on CRAN yet
#load packages
library(raster)
library(RStoolbox)
#load an example dataset
data(lsat)

 

To create some endmember spectra, we simply collect the spectral profiles of “water” and “land” from our imagery. We keep it simple here – instead, you could use spectra from a spectral library.

1
2
3
4
5
#make up some endmember spectra: water and land
em_names <- c("water", "land")
pts <- data.frame(class=em_names, cell = c(47916,5294))
em <- lsat[pts$cell]
rownames(em) <- em_names

 

That’s all for the pre-processing! Now you have an image and two endmembers. Note that you need to have at least two endmembers (as we have in this example) to unmix an image. Also, take a look at “em” and “lsat” before continuing: Both have the same spectral resolution (band number), which is a prerequisite. This means that if you want to use data from a spectral library for unmixing, you simply need to resample the data to the same spectral resolution of your imagery to use mesma().

Now, just call mesma(). It returns a probability raster, each layer representing one class of your endmembers, except for the last layer, which gives you the RMSE for each pixel.

1
2
#unmix the image for water and land
probs <- mesma(lsat, em, method = "NNLS")

 

Since we decided to develop the NNLS solver from scratch in C++ and not to use existing older FORTRAN implementations, mesma() is quiet fast. It supports parallel processing, if you want to unmix large amounts of data and have multiple CPU cores available. To do so, just create a cluster before calling mesma() using raster::beginCluster() and stop it afterwards using raster::endCluster().

Now, let us look at the output’s:

1
2
3
#take a look
raster::plot(probs$water, col = c("white","blue"))
raster::plot(probs$land, col = c("white","brown"))

 

Here are the probabilities for water…

and here they are for land:

You clearly see, how mesma() could easily differentiate between these two endmembers with only one spectrum per class. If you have ideas or find a problem, please report us via GitHub!

 

References:

Franc, V., Hlaváč, V., & Navara, M. (2005). Sequential coordinate-wise algorithm for the non-negative least squares problem. In: International Conference on Computer Analysis of Images and Patterns (pp. 407-414). Berlin, Heidelberg.

Lawson, C. L., & Hanson, R. J. (1974). Solving least squares problems (Vol. 15). Siam.

you may also like:

UAV Lidar mapping for rockfall Risk Estimation

UAV Lidar mapping for rockfall Risk Estimation

The area of the rockfall event above vineyards within a nature protection zone Due to a recent rockfall event and in order to analyze the landscape conditions within this area did we conduct a UAV Lidar mapping campaign in Thüngersheim, close to Würzburg. The main aim...

innovative urban climate in-situ measurements for Earth Observation

innovative urban climate in-situ measurements for Earth Observation

Bikair is a project aiming at measuring urban climate parameters with in-situ and Earth Observation. It focuses on testing low-cost Arduino-based sensors in an urban environment such as the city of Würzburg. Eventually, the project aims to correlate in-situ data with...

most recent news:

UAV Lidar mapping for rockfall Risk Estimation

UAV Lidar mapping for rockfall Risk Estimation

The area of the rockfall event above vineyards within a nature protection zone Due to a recent rockfall event and in order to analyze the landscape conditions within this area did we conduct a UAV Lidar mapping campaign in Thüngersheim, close to Würzburg. The main aim...

innovative urban climate in-situ measurements for Earth Observation

innovative urban climate in-situ measurements for Earth Observation

Bikair is a project aiming at measuring urban climate parameters with in-situ and Earth Observation. It focuses on testing low-cost Arduino-based sensors in an urban environment such as the city of Würzburg. Eventually, the project aims to correlate in-situ data with...