diff --git a/code/A98_PredictorsExtract.R b/code/A98_PredictorsExtract.R index ace4faa6258737425edd4cfd2f9adf43ecbf4dc2..6b3e2f3d4d52c8e28373c42789bced46865fbab7 100644 --- a/code/A98_PredictorsExtract.R +++ b/code/A98_PredictorsExtract.R @@ -18,13 +18,15 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA, require(rgeos) require(rgdal) require(geosphere) + require(spatialEco) + require(sf) print(paste("Loading", typp, "data :", toextract)) print(paste("output will be:", output)) if(typp=="raster"){ mypredictor <- raster(toextract)} else mypredictor <- readOGR(toextract) colnames(mypredictor@data) - header.shp <-readOGR(x.shp) + header.shp <- readOGR(x.shp) crs(mypredictor) <- crs(header.shp) #should be verified beforehand! @@ -54,21 +56,37 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA, if(!is.na(output)) {write.csv(out, file = output)} else{return(out)} } else { out <- sp::over(x=header.shp, y=mypredictor) - toassign <- header.shp[which(is.na(out[,1])),] + i.toassign <- which(is.na(out[,1])) + toassign <- header.shp[i.toassign,] print(paste(length(toassign), "plots not matched directly - seek for NN")) if(length(toassign)>0){ + ## explode multipolygon + mypredictor.expl <- explode(mypredictor) %>% + as_Spatial() nearest <- foreach(i=1:length(toassign), .packages=c('raster'), .combine=rbind) %dopar% { print(i) - nearest.tmp <- tryCatch(geosphere::dist2Line(header.shp[i,], mypredictor), - error = function(e){data.frame(distance=NA, lon=NA, lat=NA,ID=NA)} + ## create a subset of geoentities based on a 5° buffer radius around each target plot. + tmp.buff <- gBuffer(toassign[i,], width=5) + tryCatch( + tmp.mypredictor <- spatialEco::spatial.select( + x = tmp.buff, + y = mypredictor.expl, + distance = 0.1, + predicate = "intersect" + ), + error = function(e) { + print(paste("Nothing close enough for plot", toassign@data$PlotObservationID[i])) + } + ) + # find nearest neighbour + nearest.tmp <- tryCatch(geosphere::dist2Line(toassign[i,], tmp.mypredictor), + error = function(e){data.frame(i=i, distance=NA, lon=NA, lat=NA,ID=NA)} ) - #nearest.tmp <- geosphere::dist2Line(toassign[i,], mypredictor) return(nearest.tmp) } - - out[which(is.na(out[,1])),] <- mypredictor@data[nearest[,"ID"],] + out[i.toassign,] <- mypredictor@data[nearest[,"ID"],] } if(!is.na(output)) {write.csv(out, file = output)} else{return(out)} }