diff --git a/code/A98_PredictorsExtract.R b/code/A98_PredictorsExtract.R index 19e18c0367976ff4d5ed536da691834840ef67bd..b630aacc60b021478b88d2b4187f4203b135e7fe 100644 --- a/code/A98_PredictorsExtract.R +++ b/code/A98_PredictorsExtract.R @@ -25,11 +25,8 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA, print(paste("output will be:", output)) if(typp=="raster"){ mypredictor <- raster(toextract)} else mypredictor <- readOGR(toextract) - colnames(mypredictor@data) header.shp <- readOGR(x.shp) crs(mypredictor) <- crs(header.shp) #should be verified beforehand! - - ## Divide in chunks if requested if(chunkn>1 & !is.na(chunk.i)){ @@ -47,14 +44,16 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA, cl <- makeForkCluster(ncores, outfile="" ) registerDoParallel(cl) - print("start main foreach loop") + if(typp=="raster"){ - out <- foreach(i=1:length(header.shp), .packages=c('raster'), .combine=rbind) %dopar% { + print("start main foreach loop") + out <- foreach(i=1:length(header.shp), .packages=c('raster'), .combine=rbind) %dopar% { tmp <- raster::extract(mypredictor, header.shp[i,], buffer=min(10000, max(250, header.shp@data$lc_ncrt[i])), fun=myfunction) } if(!is.na(output)) {write.csv(out, file = output)} else{return(out)} } else { + print("Match sequentially") out <- sp::over(x=header.shp, y=mypredictor) i.toassign <- which(is.na(out[,1])) toassign <- header.shp[i.toassign,] @@ -67,7 +66,7 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA, as_Spatial() } else { mypredictor.expl <- mypredictor} - + print("start main foreach loop") nearest <- foreach(i=1:length(toassign), .packages=c('raster'), .combine=rbind) %dopar% { print(i) ## create a subset of geoentities based on a 5° buffer radius around each target plot. @@ -84,13 +83,16 @@ PredExtr <- function(x.shp, myfunction=NA, output=NA, } ) # 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 <- tryCatch(tmp.mypredictor@data[geosphere::dist2Line(toassign[i,], tmp.mypredictor)[,"ID"],], + error = function(e){ + ee <- myptedictor@data[1,] + ee[1,] <- rep(NA, ncol(mypredictor)) + } ) return(nearest.tmp) } - out[i.toassign,] <- tmp.mypredictor@data[nearest[,"ID"],] + out[i.toassign,] <- nearest } if(!is.na(output)) {write.csv(out, file = output)} else{return(out)} }