# load data
library(openxlsx)

lf.tab = read.xlsx("LANDFORM/Export_ZP_landform.xlsx")
head(lf.tab)
colnames(lf.tab)[2] = "prop_NA"

library(rgdal)
zp = readOGR("Extract_raster_zp(20_02_2020)", "CEPAZ_ZP")
head(zp@data)

library(raster)
expoeo = raster("Extract_raster_zp(20_02_2020)/expoeo.tif")
projection(expoeo)
expons = raster("Extract_raster_zp(20_02_2020)/expons.tif")
projection(expons)
tpi = raster("Extract_raster_zp(20_02_2020)/extrtpi.tif")
projection(tpi)
pente = raster("Extract_raster_zp(20_02_2020)/pentezp2.tif")
projection(pente)

# ray = raster("Extract_raster_zp(20_02_2020)/rayextr.tif")
# projection(ray) == projection(pente)
# ray2 = resample(ray, pente, method = "ngb")
# saveRDS(ray2, "ray2.rds")
ray2 = readRDS("ray2.rds")
ray2

lf = raster("LANDFORM/landformsTPI_saga.tif")
projection(lf)
lf_proj = projectRaster(lf, pente, method = "ngb")
lf.clip = crop(lf_proj, extent(pente))
# saveRDS(lf.clip, file = "lf_clip.rds")
# lf.clip = readRDS("lf_clip.rds")
lf.clip
projection(zp) == projection(expoeo)

library(sp)
zp_proj = spTransform(zp, crs(projection(pente)))

#=========================

zp_x = zp_proj[1,]
expoeo.vec = extract(expoeo, zp_x)[[1]]
expons.vec = extract(expons, zp_x)[[1]]

projection(lf.clip) == projection(zp_x)
lf.vec = extract(lf.clip, zp_x)[[1]]


tab = data.frame( expoeo = expoeo.vec, 
                  expons = expons.vec,
                  lf = ,
                  ray = extract(ray2, zp_x)[[1]],
                  pente = extract(pente, zp_x)[[1]],
                  tpi = extract(tpi, zp_x)[[1]] )
tab$lf = as.factor(tab$lf)

library(FD)
Dij = gowdis(tab)

res_x = mean(as.vector(Dij))