1. Example from https://github.com/r-hyperspec/unmixR

Both methods produce very similar results:

data("demo_data")
# Reduce data dimensionality with PCA
pca <- prcomp(demo_data)
x <- pca$x[, 1:2]
# Perform N-FINDR in reduced space
nf <- nfindr(x, p = 3)
# Get endmembers in original space
ems <- endmembers(nf, demo_data)
# Get endmembers in reduced space
ems_pca <- endmembers(nf, x)
# Calculate abundances using barycentric coordinates
# it works only in the reduced space
ab_bary <- abundances(nf, x, method = "bary")
# Calculate abundances using NNLS
# it works in both spaces but it is better to be applied in the original space
ab_nnls <- abundances(nf, demo_data, method = "nnls", normalize = TRUE)

print(summary(ab_bary))
       V1                   V2               V3        
 Min.   :-0.0008319   Min.   :0.0000   Min.   :0.0000  
 1st Qu.: 0.0135290   1st Qu.:0.1641   1st Qu.:0.2073  
 Median : 0.2637398   Median :0.2929   Median :0.2771  
 Mean   : 0.3200710   Mean   :0.3818   Mean   :0.2981  
 3rd Qu.: 0.4971472   3rd Qu.:0.6297   3rd Qu.:0.3234  
 Max.   : 1.0000000   Max.   :1.0000   Max.   :1.0000  
print(summary(ab_nnls))
       V1                V2               V3        
 Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.01664   1st Qu.:0.1449   1st Qu.:0.2062  
 Median :0.25512   Median :0.2946   Median :0.2784  
 Mean   :0.31918   Mean   :0.3810   Mean   :0.2998  
 3rd Qu.:0.48584   3rd Qu.:0.6305   3rd Qu.:0.3291  
 Max.   :1.00000   Max.   :1.0000   Max.   :1.0000  
ggplot(data=data.frame(bary=ab_bary[,1], nnls=ab_nnls[,1])) +
  geom_point(aes(x=bary, y=nnls))

2. ASTER image

“bary” and “nnls” methods produce very different result.

#remotes::install_github("r-hyperspec/unmixR")
library(terra)  
library(RStoolbox)
library(unmixR)
library(ggplot2)

rdir <- "../Tazoult_AST/AST-07XT_20020524"

2.1 Data input

ast <- rast(file.path(rdir, "AST_07XT_VNS2_aoi_sh.tif"))
#ast <- rast(file.path(rdir, "AST_07XT_VNS2_aoi_pangfFRccsm20.tif")) #COMPTE: NA values!!!
ast
names(ast) <- paste0("B0",1:9)
print(paste0("NA present? ", any(is.na(values(ast)))))
#table(values(ast==0))
#range(values(ast))
wl1 <- c(0.52, 0.63, 0.76, 1.6, 2.145, 2.185, 2.235, 2.295, 2.360)
wl2 <- c(0.6, 0.69, 0.86, 1.7, 2.185, 2.225, 2.285, 2.365, 2.430)
wl <- (wl2-wl1)/2 + wl1

2.2 PCA

X <- values(ast)
print(summary(prcomp(X)))
Importance of components:
                           PC1      PC2      PC3     PC4    PC5     PC6     PC7     PC8     PC9
Standard deviation     88.2258 24.53153 16.43535 7.18504 5.4601 3.34969 2.92804 2.24463 2.07628
Proportion of Variance  0.8879  0.06865  0.03081 0.00589 0.0034 0.00128 0.00098 0.00057 0.00049
Cumulative Proportion   0.8879  0.95657  0.98739 0.99327 0.9967 0.99796 0.99893 0.99951 1.00000
X_pc_demo <- prcomp(X)$x

2.3. Perform N-FINDR in reduced space

p=7
Xnfindrfinal<- nfindr(X_pc_demo[,1:p], p = p+1, n_init=50)
e <- endmembers(Xnfindrfinal, X)
df <- reshape2::melt(e)
colnames(df) <- c("EM", "Band", "Reflectance")
df$EM <- as.character(df$EM)
dict <- data.frame(Band=unique(df$Band), Wavelength=wl)
df$Wavelength <- plyr::mapvalues(df$Band, from=dict$Band, to=dict$Wavelength)
df$Wavelength <- as.numeric(as.character(df$Wavelength))
ggplot(data=df) + 
  geom_line(aes(x=Wavelength, y=Reflectance, group=EM,color=EM)) +
  facet_wrap(~EM,ncol=3)

2.4. Abundances

abnd_bary <- abundances(Xnfindrfinal, X_pc_demo[,1:7], method = "bary")
# Calculate abundances using NNLS
# it works in both spaces but it is better to be applied in the original space
abnd_nnls <- abundances(Xnfindrfinal, X, method = "nnls", normalize = TRUE)

2.5 Compare “bary” to “nnls”

options(width = 200)
print(summary(abnd_bary[,1:5])) #select 1:5 for a cleaner output
       V1                 V2                 V3                  V4                V5          
 Min.   :-0.39599   Min.   :-0.52147   Min.   :-0.534884   Min.   :-0.3615   Min.   :-0.52115  
 1st Qu.:-0.08633   1st Qu.: 0.07177   1st Qu.: 0.006923   1st Qu.: 0.2647   1st Qu.:-0.04271  
 Median :-0.03632   Median : 0.15139   Median : 0.095025   Median : 0.3428   Median : 0.01783  
 Mean   :-0.02787   Mean   : 0.14933   Mean   : 0.096283   Mean   : 0.3474   Mean   : 0.02192  
 3rd Qu.: 0.02484   3rd Qu.: 0.22475   3rd Qu.: 0.182124   3rd Qu.: 0.4295   3rd Qu.: 0.07937  
 Max.   : 1.00000   Max.   : 1.00000   Max.   : 1.000000   Max.   : 1.0000   Max.   : 1.00000  
print(summary(abnd_nnls[,1:5]))
       V1                V2                V3                V4               V5         
 Min.   :0.00000   Min.   :0.00000   Min.   :0.00000   Min.   :0.0000   Min.   :0.00000  
 1st Qu.:0.00000   1st Qu.:0.00000   1st Qu.:0.02470   1st Qu.:0.1975   1st Qu.:0.00000  
 Median :0.00000   Median :0.04996   Median :0.09918   Median :0.2943   Median :0.02202  
 Mean   :0.02089   Mean   :0.07052   Mean   :0.12047   Mean   :0.2899   Mean   :0.06056  
 3rd Qu.:0.03304   3rd Qu.:0.12889   3rd Qu.:0.18378   3rd Qu.:0.3891   3rd Qu.:0.09744  
 Max.   :1.00000   Max.   :1.00000   Max.   :1.00000   Max.   :1.0000   Max.   :1.00000  
df <- data.frame(bary=abnd_bary[,1], nnls=abnd_nnls[,1]) 
ggplot(data=df) +
  geom_bin2d(aes(x=bary, y=nnls),bins=70)

  • Abundances <0 with bary method
  • Resulting Abundances are very different
LS0tCnRpdGxlOiAiTkZJTkRSX2xvZy5SbWQiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCiMjIDEuIEV4YW1wbGUgZnJvbSBodHRwczovL2dpdGh1Yi5jb20vci1oeXBlcnNwZWMvdW5taXhSCgpCb3RoIG1ldGhvZHMgcHJvZHVjZSB2ZXJ5IHNpbWlsYXIgcmVzdWx0czoKCmBgYHtyLGZpZy53aWR0aD00LCBmaWcuaGVpZ2h0PTMsICBmaWcuc2hvdz0naG9sZCd9CmRhdGEoImRlbW9fZGF0YSIpCiMgUmVkdWNlIGRhdGEgZGltZW5zaW9uYWxpdHkgd2l0aCBQQ0EKcGNhIDwtIHByY29tcChkZW1vX2RhdGEpCnggPC0gcGNhJHhbLCAxOjJdCiMgUGVyZm9ybSBOLUZJTkRSIGluIHJlZHVjZWQgc3BhY2UKbmYgPC0gbmZpbmRyKHgsIHAgPSAzKQojIEdldCBlbmRtZW1iZXJzIGluIG9yaWdpbmFsIHNwYWNlCmVtcyA8LSBlbmRtZW1iZXJzKG5mLCBkZW1vX2RhdGEpCiMgR2V0IGVuZG1lbWJlcnMgaW4gcmVkdWNlZCBzcGFjZQplbXNfcGNhIDwtIGVuZG1lbWJlcnMobmYsIHgpCiMgQ2FsY3VsYXRlIGFidW5kYW5jZXMgdXNpbmcgYmFyeWNlbnRyaWMgY29vcmRpbmF0ZXMKIyBpdCB3b3JrcyBvbmx5IGluIHRoZSByZWR1Y2VkIHNwYWNlCmFiX2JhcnkgPC0gYWJ1bmRhbmNlcyhuZiwgeCwgbWV0aG9kID0gImJhcnkiKQojIENhbGN1bGF0ZSBhYnVuZGFuY2VzIHVzaW5nIE5OTFMKIyBpdCB3b3JrcyBpbiBib3RoIHNwYWNlcyBidXQgaXQgaXMgYmV0dGVyIHRvIGJlIGFwcGxpZWQgaW4gdGhlIG9yaWdpbmFsIHNwYWNlCmFiX25ubHMgPC0gYWJ1bmRhbmNlcyhuZiwgZGVtb19kYXRhLCBtZXRob2QgPSAibm5scyIsIG5vcm1hbGl6ZSA9IFRSVUUpCgpwcmludChzdW1tYXJ5KGFiX2JhcnkpKQpwcmludChzdW1tYXJ5KGFiX25ubHMpKQpnZ3Bsb3QoZGF0YT1kYXRhLmZyYW1lKGJhcnk9YWJfYmFyeVssMV0sIG5ubHM9YWJfbm5sc1ssMV0pKSArCiAgZ2VvbV9wb2ludChhZXMoeD1iYXJ5LCB5PW5ubHMpKQpgYGAKCiMgMi4gQVNURVIgaW1hZ2UKCiJiYXJ5IiBhbmQgIm5ubHMiIG1ldGhvZHMgcHJvZHVjZSB2ZXJ5IGRpZmZlcmVudCByZXN1bHQuCgpgYGB7ciwgIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIHJlc3VsdHM9J2hpZGUnfQojcmVtb3Rlczo6aW5zdGFsbF9naXRodWIoInItaHlwZXJzcGVjL3VubWl4UiIpCmxpYnJhcnkodGVycmEpICAKbGlicmFyeShSU3Rvb2xib3gpCmxpYnJhcnkodW5taXhSKQpsaWJyYXJ5KGdncGxvdDIpCgpyZGlyIDwtICIuLi9UYXpvdWx0X0FTVC9BU1QtMDdYVF8yMDAyMDUyNCIKYGBgCgojIyMgMi4xIERhdGEgaW5wdXQKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCByZXN1bHRzPSdoaWRlJ30KYXN0IDwtIHJhc3QoZmlsZS5wYXRoKHJkaXIsICJBU1RfMDdYVF9WTlMyX2FvaV9zaC50aWYiKSkKI2FzdCA8LSByYXN0KGZpbGUucGF0aChyZGlyLCAiQVNUXzA3WFRfVk5TMl9hb2lfcGFuZ2ZGUmNjc20yMC50aWYiKSkgI0NPTVBURTogTkEgdmFsdWVzISEhCmFzdApuYW1lcyhhc3QpIDwtIHBhc3RlMCgiQjAiLDE6OSkKcHJpbnQocGFzdGUwKCJOQSBwcmVzZW50PyAiLCBhbnkoaXMubmEodmFsdWVzKGFzdCkpKSkpCiN0YWJsZSh2YWx1ZXMoYXN0PT0wKSkKI3JhbmdlKHZhbHVlcyhhc3QpKQp3bDEgPC0gYygwLjUyLCAwLjYzLCAwLjc2LCAxLjYsIDIuMTQ1LCAyLjE4NSwgMi4yMzUsIDIuMjk1LCAyLjM2MCkKd2wyIDwtIGMoMC42LCAwLjY5LCAwLjg2LCAxLjcsIDIuMTg1LCAyLjIyNSwgMi4yODUsIDIuMzY1LCAyLjQzMCkKd2wgPC0gKHdsMi13bDEpLzIgKyB3bDEKYGBgCgojIyMgMi4yIFBDQQoKYGBge3J9ClggPC0gdmFsdWVzKGFzdCkKcHJpbnQoc3VtbWFyeShwcmNvbXAoWCkpKQpYX3BjX2RlbW8gPC0gcHJjb21wKFgpJHgKYGBgCgojIyMgMi4zLiBQZXJmb3JtIE4tRklORFIgaW4gcmVkdWNlZCBzcGFjZQoKYGBge3J9Cm5jcD03ClhuZmluZHJmaW5hbDwtIG5maW5kcihYX3BjX2RlbW9bLDE6bmNwXSwgcCA9IG5jcCsxLCBuX2luaXQ9NTApCmBgYAoKYGBge3J9CmUgPC0gZW5kbWVtYmVycyhYbmZpbmRyZmluYWwsIFgpCmRmIDwtIHJlc2hhcGUyOjptZWx0KGUpCmNvbG5hbWVzKGRmKSA8LSBjKCJFTSIsICJCYW5kIiwgIlJlZmxlY3RhbmNlIikKZGYkRU0gPC0gYXMuY2hhcmFjdGVyKGRmJEVNKQpkaWN0IDwtIGRhdGEuZnJhbWUoQmFuZD11bmlxdWUoZGYkQmFuZCksIFdhdmVsZW5ndGg9d2wpCmRmJFdhdmVsZW5ndGggPC0gcGx5cjo6bWFwdmFsdWVzKGRmJEJhbmQsIGZyb209ZGljdCRCYW5kLCB0bz1kaWN0JFdhdmVsZW5ndGgpCmRmJFdhdmVsZW5ndGggPC0gYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIoZGYkV2F2ZWxlbmd0aCkpCmdncGxvdChkYXRhPWRmKSArIAogIGdlb21fbGluZShhZXMoeD1XYXZlbGVuZ3RoLCB5PVJlZmxlY3RhbmNlLCBncm91cD1FTSxjb2xvcj1FTSkpICsKICBmYWNldF93cmFwKH5FTSxuY29sPTMpCmBgYAoKIyMjIDIuNC4gQWJ1bmRhbmNlcwoKYGBge3J9CmFibmRfYmFyeSA8LSBhYnVuZGFuY2VzKFhuZmluZHJmaW5hbCwgWF9wY19kZW1vWywxOjddLCBtZXRob2QgPSAiYmFyeSIpCiMgQ2FsY3VsYXRlIGFidW5kYW5jZXMgdXNpbmcgTk5MUwojIGl0IHdvcmtzIGluIGJvdGggc3BhY2VzIGJ1dCBpdCBpcyBiZXR0ZXIgdG8gYmUgYXBwbGllZCBpbiB0aGUgb3JpZ2luYWwgc3BhY2UKYWJuZF9ubmxzIDwtIGFidW5kYW5jZXMoWG5maW5kcmZpbmFsLCBYLCBtZXRob2QgPSAibm5scyIsIG5vcm1hbGl6ZSA9IFRSVUUpCmBgYAojIyMgMi41IENvbXBhcmUgImJhcnkiIHRvICJubmxzIgoKYGBge3IsZmlnLndpZHRoPTQsIGZpZy5oZWlnaHQ9MywgIGZpZy5zaG93PSdob2xkJ30Kb3B0aW9ucyh3aWR0aCA9IDIwMCkKcHJpbnQoc3VtbWFyeShhYm5kX2JhcnlbLDE6NV0pKSAjc2VsZWN0IDE6NSBmb3IgYSBjbGVhbmVyIG91dHB1dApwcmludChzdW1tYXJ5KGFibmRfbm5sc1ssMTo1XSkpCmRmIDwtIGRhdGEuZnJhbWUoYmFyeT1hYm5kX2JhcnlbLDFdLCBubmxzPWFibmRfbm5sc1ssMV0pIApnZ3Bsb3QoZGF0YT1kZikgKwogIGdlb21fYmluMmQoYWVzKHg9YmFyeSwgeT1ubmxzKSxiaW5zPTcwKQoKYGBgCiogQWJ1bmRhbmNlcyA8MCB3aXRoIGJhcnkgbWV0aG9kCiogUmVzdWx0aW5nIEFidW5kYW5jZXMgYXJlIHZlcnkgZGlmZmVyZW50IAo=