Commit 6650dcbc authored by Leila Schuh's avatar Leila Schuh
Browse files

structural diversity entropy

parent e75563ba
.Rproj.user
.Rhistory
.RData
.Ruserdata
#=================================================================================================
# Bin NDVI Data
#=================================================================================================
library(raster)
library(radiomics)
#=================================================================================================
# NDVI Data has been downloaded from google earth engine, representing yearly average NDVI values
# over the growing season of 2018, in a study region in north eastern Eurasia (NEE).
# The few NA gaps were filled with a local neighborhood average.
#=================================================================================================
nee <- raster("Data/NDVI.NAfill.2018.nee.tif")
neemat <- matrix(nee, nrow=nrow(nee), ncol =ncol(nee), byrow = TRUE)
#==============================================================================================
# The following function was copied from thee radiomics package, which was available on CRAN
# at the time this experiment was carried out, but has since been removed from CRAN.
#==============================================================================================
discretizeImage <- function(data, n_grey=32, verbose=TRUE){
#Error checking
#Check validity of input
if (!is.matrix(data)) {
stop(paste0("Object of class ", class(data), ". is.matrix(object) must evaluate TRUE."))
}
if (any(data < 0, na.rm=TRUE)) {
stop("Object contains negative values. All values must be greater than 0.")
}
#Not a perfect solution. Makes n_grey breaks, but doesn't necessarily populate all of them
# eg. n_gey could be 100, but only 75 of the levels are used by pixels
l_unique <- length(unique(c(data)))
#Make sure discretization is valid
if(l_unique == 0 ) stop("Function not valid for empty input")
if(sum(is.na(data)) == length(data)){
if(verbose) warning("Matrix composed entirely of NA's\nReturning Matrix")
return(data)
}
#If we don't need to do anything, we don't do anything
if(n_grey == l_unique){
return(data)
} else if(n_grey > l_unique){
if(verbose) message(sprintf("n_grey (%d) cannot be larger than the number of gray levels in the image (%d)", n_grey, l_unique))
if(verbose) message(sprintf("n_grey set to %d", l_unique))
return(data)
}
discretized <- cut(data, breaks=seq(min(data, na.rm=TRUE), max(data, na.rm=TRUE), length.out=(n_grey + 1)),
labels = seq(1, n_grey, 1),
include.lowest=TRUE, right=FALSE)
return(matrix(as.numeric(discretized), nrow=nrow(data)))
}
#==============================================================================================
# Bin NDVI data into different equal-size bins and save data as .tif files
#==============================================================================================
nee.1000GL <- discretizeImage(neemat, n_grey = 1000)
nee.1000GL <- raster(nee.100GL0, crs = crs(nee))
extent(nee.1000GL) <- extent(nee)
writeRaster(nee.1000GL, "Data/nee.1000GL.tif")
nee.500GL <- discretizeImage(neemat, n_grey = 500)
nee.500GL <- raster(nee.500GL, crs = crs(nee))
extent(nee.500GL) <- extent(nee)
writeRaster(nee.500GL, "Data/nee.500GL.tif")
nee.100GL <- discretizeImage(neemat, n_grey = 100)
nee.100GL <- raster(nee.100GL, crs = crs(nee))
extent(nee.100GL) <- extent(nee)
writeRaster(nee.100GL, "Data/nee.100GL.tif")
nee.50GL <- discretizeImage(neemat, n_grey = 50)
nee.50GL <- raster(nee.50GL, crs = crs(nee))
extent(nee.50GL) <- extent(nee)
writeRaster(nee.50GL, "Data/nee.50GL.tif")
nee.25GL <- discretizeImage(neemat, n_grey = 25)
nee.25GL <- raster(nee.25GL, crs = crs(nee))
extent(nee.25GL) <- extent(nee)
writeRaster(nee.25GL, "Data/nee.25GL.tif")
nee.15GL <- discretizeImage(neemat, n_grey = 15)
nee.15GL <- raster(nee.15GL, crs = crs(nee))
extent(nee.15GL) <- extent(nee)
writeRaster(nee.15GL, "Data/nee.15GL.tif")
nee.10GL <- discretizeImage(neemat, n_grey = 10)
nee.10GL <- raster(nee.10GL, crs = crs(nee))
extent(nee.10GL) <- extent(nee)
writeRaster(nee.10GL, "Data/nee.10GL.tif")
nee.9GL <- discretizeImage(neemat, n_grey = 9)
nee.9GL <- raster(nee.9GL, crs = crs(nee))
extent(nee.9GL) <- extent(nee)
writeRaster(nee.9GL, "Data/nee.9GL.tif")
nee.8GL <- discretizeImage(neemat, n_grey = 8)
nee.8GL <- raster(nee.8GL, crs = crs(nee))
extent(nee.8GL) <- extent(nee)
writeRaster(nee.8GL, "Data/nee.8GL.tif")
nee.7GL <- discretizeImage(neemat, n_grey = 7)
nee.7GL <- raster(nee.7GL, crs = crs(nee))
extent(nee.7GL) <- extent(nee)
writeRaster(nee.7GL, "Data/nee.7GL.tif")
nee.6GL <- discretizeImage(neemat, n_grey = 6)
nee.6GL <- raster(nee.6GL, crs = crs(nee))
extent(nee.6GL) <- extent(nee)
writeRaster(nee.6GL, "Data/nee.6GL.tif")
nee.5GL <- discretizeImage(neemat, n_grey = 5)
nee.5GL <- raster(nee.5GL, crs = crs(nee))
extent(nee.5GL) <- extent(nee)
writeRaster(nee.5GL, "Data/nee.5GL.tif")
nee.4GL <- discretizeImage(neemat, n_grey = 4)
nee.4GL <- raster(nee.4GL, crs = crs(nee))
extent(nee.4GL) <- extent(nee)
writeRaster(nee.4GL, "Data/nee.4GL.tif")
nee.3GL <- discretizeImage(neemat, n_grey = 3)
nee.3GL <- raster(nee.3GL, crs = crs(nee))
extent(nee.3GL) <- extent(nee)
writeRaster(nee.3GL, "Data/nee.3GL.tif")
nee.2GL <- discretizeImage(neemat, n_grey = 2)
nee.2GL <- raster(nee.2GL, crs = crs(nee))
extent(nee.2GL) <- extent(nee)
writeRaster(nee.2GL, "Data/nee.2GL.tif")
#=================================================================================================
\ No newline at end of file
#=================================================================================================
# Structural diversity of NDVI data
#=================================================================================================
# This is the template script for structural diversity estimation in NDVI data, using
# window side length (WSL) = 3
# For all other WLS ∈ W1={3, 5, 7, 9, 11, 13, 15, 17, 19, 25, 35, 45}, please replace arguments
# "NgbProb" and "NgbHet" with the respective WSL, and adapt file names.
# For all binned data, please load the respective file instead, and adapt file names.
# All model settings remain the same
# Metrics tested: entropy, weighted entropies, contrast, dissimilarity, homogeneity
# "Spectral" means "value", and "Position" means "rank" in the function settings.
# Distance ("Dist") 1, all directions ("angle") specifies the spatial relationship of pixels.
# Please use the .tar.gz file to install the R package StrucDivDevel.
# If any NAs occur in the window, NA is returned ("na.rm = FALSE").
#=================================================================================================
library(raster)
library(StructDivDevel)
#=================================================================================================
dir.create("Methods/StructuralDiversity/tif")
#=================================================================================================
nee <- raster("Data/NDVI.NAfill.2018.nee.tif")
#=================================================================================================
# WSL 3
#=================================================================================================
# Entropy
nee.Ent.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "Entropy",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
# Weighted Entropy
nee.AbsValEnt.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "WeightedEntropyAbsSpectral",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
nee.SqrValEnt.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "WeightedEntropySqrSpectral",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
nee.AbsRankEnt.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "WeightedEntropyAbsPosition",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
nee.SqrRankEnt.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "WeightedEntropySqrPosition",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
# Contrast
nee.ValCon.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "ContrastSpectral",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
nee.RankCon.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "ContrastPosition",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
# Dissimilarity
nee.ValDis.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "DissimilaritySpectral",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
nee.RankDis.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "DissimilarityPosition",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
# Homogeneity
nee.Hom.3 <- FocalRcppNew(x = nee, NgbProb = 3, NgbHet = 3, fun = "HomogeneityPosition",
Dist = 1, na.rm = FALSE, angle = "all", parallelize = TRUE)
#=================================================================================================
writeRaster(nee.Ent.3, "Methods/StructuralDiversity/tif/nee.Ent.3.tif")
writeRaster(nee.AbsValEnt.3, "Methods/StructuralDiversity/tif/nee.AbsValEnt.3.tif")
writeRaster(nee.SqrValEnt.3, "Methods/StructuralDiversity/tif/nee.SqrValEnt.3.tif")
writeRaster(nee.AbsRankEnt.3, "Methods/StructuralDiversity/tif/nee.AbsRankEnt.3.tif")
writeRaster(nee.SqrRankEnt.3, "Methods/StructuralDiversity/tif/nee.SqrRankEnt.3.tif")
writeRaster(nee.RankCon.3, "Methods/StructuralDiversity/tif/nee.RankCon.3.tif")
writeRaster(nee.ValDis.3, "Methods/StructuralDiversity/tif/nee.ValDis.3.tif")
writeRaster(nee.RankDis.3, "Methods/StructuralDiversity/tif/nee.RankDis.3.tif")
writeRaster(nee.Hom.3, "Methods/StructuralDiversity/tif/nee.Hom.3.tif")
#=================================================================================================
#=================================================================================================
# Resampling simulation
#=================================================================================================
# This is the template script for resampling simulation with function entropy
# For all other metrics, please run the script with other function names, as can bee seen
# at the bottom of this script
#=================================================================================================
library(raster)
library(radiomics)
#=================================================================================================
dir.create("Methods/ResamplingSimulation/csv")
set.seed(42)
#=================================================================================================
# Sampling from NDVI data
#=================================================================================================
nee <- raster("Data/NDVI.NAfill.2018.nee.tif")
nee.mat <- matrix(nee, nrow=nrow(nee), ncol =ncol(nee), byrow = TRUE)
#=================================================================================================
# Function to calculate co-occurence matrix
#=================================================================================================
# x - data/ sample
# WSL - window side length
# d - distance between pixels, here 1, direct pixel neighbors
#=================================================================================================
# Pixel pairs:
# Nrp - number of pixel pairs, we consider all directions.
# We consider pixel pairs within B_k. Therefore, ncol(B_k) = nrow(B_k) = WSL =>
# Nrp = ((WSL - d) * WSL + (WSL -d )* WSL + (WSL - d) * (WSL - d) * 2) * 2
# Nrp = 4 * (WSL - d) * (2 * WSL - d)
#=================================================================================================
# The following function returns the gray level co-occurrence matrix (GLCM)
#=================================================================================================
PMat <- function(x, WSL, d){
xMat <- matrix(x, nrow = WSL, ncol =WSL, byrow=T)
Vals <- unique(sort((xMat)))
SpaDe <- matrix(0, nrow = length(Vals), ncol = length(Vals))
colnames(SpaDe) <- rownames(SpaDe) <- Vals
for(i in 1:nrow(SpaDe)){
for(j in 1:ncol(SpaDe)){
for(a in 1:nrow(xMat)){
for(b in 1:ncol(xMat)){
if(b <= (ncol(xMat) - d)){
if(Vals[[i]] == xMat[[a,b]] & Vals[[j]] == xMat[[a,b+d]]){
SpaDe[[i,j]] <- SpaDe[[i,j]] + 1
}
}
if(b > d & a <= (nrow(xMat) - d)){
if(Vals[[i]] == xMat[[a,b]] & Vals[[j]] == xMat[[a+d,b-d]]){
SpaDe[[i,j]] <- SpaDe[[i,j]] + 1
}
}
if(a <= (nrow(xMat) - d)){
if(Vals[[i]] == xMat[[a,b]] & Vals[[j]] == xMat[[a+d,b]]){
SpaDe[[i,j]] <- SpaDe[[i,j]] + 1
}
}
if(b <= (ncol(xMat) - d) & a <= (nrow(xMat) - d)){
if(Vals[[i]] == xMat[[a,b]] & Vals[[j]] == xMat[[a+d,b+d]]){
SpaDe[[i,j]] <- SpaDe[[i,j]] + 1
}
}
}
}
}
}
SpaDe <- t(SpaDe) + SpaDe
Nrp <- 4 * (WSL - d) * (2*WSL - d)
SpatMatP <- SpaDe/Nrp
return(SpatMatP)
}
# test. If sum(p) = 1, everything is good
# p <- PMat(nee.mat, 3, 1)
# sum(p)
#=================================================================================================
# Diversity functions
#=================================================================================================
# The following functions calculate the diversity metric for every element of the co-occurrence matrix
# (enters the function as pmat, which is the GLCM), which results in the metric matrix. The sum of all
# elements of the metric matrix is returned.
#=================================================================================================
# Entropy
#=================================================================================================
Entropy <- function(pmat){
Hmat <- pmat * (-log(pmat))
IndivEntropy<- c(Hmat)
Entropy <- (sum(Hmat, na.rm = TRUE))
out <- as.data.frame(cbind(IndivEntropy, Entropy))
return(out)
}
#=================================================================================================
# Rank dissimilarity
#=================================================================================================
RankDissim <- function(pmat){
Dmat <- matrix(0, nrow = nrow(pmat), ncol = ncol(pmat))
rownames(Dmat) <- colnames(Dmat) <- rownames(pmat)
for (r in 1:nrow(pmat)){
for (s in 1:ncol(pmat)){
Dmat[r,s] <- pmat[r,s] * abs(r-s)
}
}
IndivRankDissim<- c(Dmat)
RankDissim <- (sum(Dmat, na.rm = TRUE))
out <- as.data.frame(cbind(IndivRankDissim, RankDissim))
return(out)
}
#=================================================================================================
# Rank contrast
#=================================================================================================
RankContrast <- function(pmat){
Cmat <- matrix(0, nrow = nrow(pmat), ncol = ncol(pmat))
rownames(Cmat) <- colnames(Cmat) <- rownames(pmat)
for (r in 1:nrow(pmat)){
for (s in 1:ncol(pmat)){
Cmat[r,s] <- pmat[r,s] * (r-s)^2
}
}
IndivRankContrast <- c(Cmat)
RankContrast <- (sum(Cmat, na.rm = TRUE))
out <- as.data.frame(cbind(IndivRankContrast, RankContrast))
return(out)
}
#=================================================================================================
# Homogeneity
#=================================================================================================
Homogen <- function(pmat){
HoMat <- matrix(0, nrow = nrow(pmat), ncol = ncol(pmat))
rownames(HoMat) <- colnames(HoMat) <- rownames(HoMat)
for (r in 1:nrow(pmat)){
for (s in 1:ncol(pmat)){
HoMat[r,s] <- pmat[r,s] / (1 + (r-s)^2)
}
}
IndivHomogen <- c(HoMat)
Homogen <- (sum(HoMat, na.rm = TRUE))
out <- as.data.frame(cbind(IndivHomogen, Homogen))
return(out)
}
#=================================================================================================
# Square rank weighted entropy
#=================================================================================================
SqrRankEnt <- function(pmat){
SqrRankEntMat <- matrix(0, nrow = nrow(pmat), ncol = ncol(pmat))
rownames(SqrRankEntMat) <- colnames(SqrRankEntMat) <- rownames(pmat)
for (r in 1:nrow(pmat)){
for (s in 1:ncol(pmat)){
SqrRankEntMat[r,s] <- pmat[r,s] * (-log(pmat[r,s])) * (r-s)^2
}
}
IndivSqrRankEnt <- c(SqrRankEntMat)
SqrRankEnt <- (sum(SqrRankEntMat, na.rm = TRUE))
out <- as.data.frame(cbind(IndivSqrRankEnt, SqrRankEnt))
return(out)
}
#=================================================================================================
# Absolute rank weighted entropy
#=================================================================================================
AbsRankEnt <- function(pmat){
AbsRankEntMat <- matrix(0, nrow = nrow(pmat), ncol = ncol(pmat))
rownames(AbsRankEntMat) <- colnames(AbsRankEntMat) <- rownames(pmat)
for (r in 1:nrow(pmat)){
for (s in 1:ncol(pmat)){
AbsRankEntMat[r,s] <- pmat[r,s] * (-log(pmat[r,s])) * abs(r-s)
}
}
IndivAbsRankEnt <- c(AbsRankEntMat)
AbsRankEnt <- (sum(AbsRankEntMat, na.rm = TRUE))
out <- as.data.frame(cbind(IndivAbsRankEnt, AbsRankEnt))
return(out)
}
#=================================================================================================
# Simulate windows B_k
#=================================================================================================
# rep - nr. of times B_k should be sampled
# tab = table of raster matrix, columns renamed accordingly
# wsl = window side length
# d = distance between pixels (1)
# PFun - PMat (function to calculate GLCM)
# m = number of GL data has been binned to
# div = diversity function
#=================================================================================================
# The following function samples the moving window B_k rep times and calculates the specified
# diversity metric
#=================================================================================================
SimB <- function(rep, tab, wsl, d, PFun, m, div){
B_out <- data.frame()
m_i <- replicate(rep, sample(x = tab$M_i, size = wsl^2, replace = T, prob = tab$P_i))
for(i in 1:rep){
pmat <- PFun(m_i[,i], wsl, d)
metric <- div(pmat)
P_ij <- as.vector(pmat)
df <- as.data.frame(cbind(P_ij, metric, m, wsl))
B_out <- rbind(B_out, df)
}
return(B_out)
}
#=================================================================================================
# Complete resampling simulation
#=================================================================================================
# M - vector. GL to be sampled
# WSL - vector. WSLs to be sampled
# D - distance between pixel pairs (here 1)
# Mat - data to be sampled from as matrix (hier nee.mat)
# BFun - SimB, function to sample B_k
# Div - Diversity function
# Rep - number of times B_k shall be sampled
#=================================================================================================
# the following function runs the whole resampling simulation
#=================================================================================================
Sim <- function(M, WSL, D, Mat, BFun, Pmat, Div, Rep){
out <- data.frame()
for (mm in M){
A <- discretizeImage(Mat, n_grey = mm)
tabA <- as.data.frame(table(A))
tabA <- cbind(tabA, tabA$Freq/sum(tabA$Freq))
colnames(tabA) <- c("M_i","Freq_i","P_i")
for (ws in WSL){
HB <- SimB(rep = Rep, tab = tabA, wsl = ws, d = D, PFun = Pmat, m = mm, div = Div)
out <- rbind(out, HB)
}
}
return(out)
}
#=================================================================================================
# Set parameters
#=================================================================================================
Ms <- c(2, 4, 6, 8, 10, 12, 14, 16, 18, 20,
25, 50, 100, 250, 500, 1000)
WSLs <- c(3, 7, 11, 15, 19, 25, 45)
times <- 10
distance <- 1
#=================================================================================================
# Run simulations
#=================================================================================================
# Entropy
#=================================================================================================
Entropy <- Sim(M = Ms, WSL = WSLs, D = distance, Mat = nee.mat, BFun = SimB,
Pmat = PMat, Div = Entropy, Rep = times)
write.table(Entropy, "Methods/ResamplingSimulation/csv/Entropy.csv", row.names = FALSE, col.names = TRUE)
#=================================================================================================
# For all other metrics run the respective function as seen below:
##=================================================================================================
# # Absolute rank entropy
#=================================================================================================
# AbsRankEntSim <- Sim(M = Ms, WSL = WSLs, D = distance, Mat = nee.mat, BFun = SimB,
# Pmat = PMat, Div = AbsRankEnt, Rep = times)
#
# write.table(AbsRankEntSim, "Methods/ResamplingSimulation/csv/AbsRankEnt.csv", row.names = FALSE, col.names = TRUE)
#
#=================================================================================================
# # Square rank entropy
#=================================================================================================
#
# SqrRankEntSim <- Sim(M = Ms, WSL = WSLs, D = distance, Mat = nee.mat, BFun = SimB,
# Pmat = PMat, Div = SqrRankEnt, Rep = times)
#
# write.table(SqrRankEntSim, "Methods/ResamplingSimulation/csv/SqrRankEnt.csv", row.names = FALSE, col.names = TRUE)
#
#=================================================================================================
# # Rank contrast
#=================================================================================================
#
# RankContrastSim <- Sim(M = Ms, WSL = WSLs, D = distance, Mat = nee.mat, BFun = SimB,
# Pmat = PMat, Div = RankContrast, Rep = times)
#
# write.table(RankContrastSim, "Methods/ResamplingSimulation/csv/RankContrast.csv", row.names = FALSE, col.names = TRUE)
#
#=================================================================================================
# # Rank Dissimilarity
# #=============================================================================================================================
#
# RankDisSim <- Sim(M = Ms, WSL = WSLs, D = distance, Mat = nee.mat, BFun = SimB,
# Pmat = PMat, Div = RankDissim, Rep = times)
#
# write.table(RankDisSim, "Methods/ResamplingSimulation/csv/RankDissim.csv", row.names = FALSE, col.names = TRUE)
#
#=================================================================================================
# # Homogeneity
#=================================================================================================
#
# HomogenSim <- Sim(M = Ms, WSL = WSLs, D = distance, Mat = nee.mat, BFun = SimB,
# Pmat = PMat, Div = Homogen, Rep = times)
#
# write.table(HomogenSim, "Methods/ResamplingSimulation/csv/Homogen.csv", row.names = FALSE, col.names = TRUE)
#
#=================================================================================================
\ No newline at end of file