# 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))