使用`sp`和`sf`使R脚本适应管理区域的问题

如何解决使用`sp`和`sf`使R脚本适应管理区域的问题

我正在尝试在 Qgis 3.14 上运行 R 脚本，该脚本在 Qgis 2.18 版中运行良好。它根据 DEM 地形高程栅格文件划分管理区域。 但是我遇到了一些问题。 这是我的脚本

##Determina_zonas_de_manejo=group
##Raster_de_EnTrada=raster
##Quantidade_de_Zonas_de_Manejo=number 2
##Poligono_Unidade_Produtiva=optional vector
##Ponto_Inicial=optional vector
##Inicial_do_identificador=string A
##Porcentagem_de_vertices_a_serem_mantidos_0_a_100=number 100
##Zonas_de_Manejo=output vector

library(raster)
library(sp)
library(rgdal)
library(e1071)
library(igraph)
library(rgeos)
library(TSP)

raster_in  <- Raster_de_EnTrada[[1]]
groups <- Quantidade_de_Zonas_de_Manejo
pol_mask <- Poligono_Unidade_Produtiva
pts_inicial <- Ponto_Inicial
iniID <- Inicial_do_identificador

#Porcentagem_de_vertices_a_serem_mantidos_0_a_100 <- 100
porc_vert <- Porcentagem_de_vertices_a_serem_mantidos_0_a_100/100

#raster_in <- brick('G:\\Meu Drive\\Drive\\RemoteSensing\\UGD_v2\\correcoes14_02_2020\\ArquivosRenomearZM\\NDVI_Final_01.tif')[[1]]
#groups <- 5
#pol_file <- 'G:\\Meu Drive\\Drive\\RemoteSensing\\UGD_v2\\correcoes14_02_2020\\ArquivosRenomearZM\\01.shp'
#pol_mask <- readOGR(dirname(pol_file),strsplit(basename(pol_file),'[.]')[[1]][1] )
#iniID <- "B"
#pts_file <- 'G:\\Meu Drive\\Drive\\RemoteSensing\\UGD_v2\\correcoes14_02_2020\\ArquivosRenomearZM\\ponto_inicial_teste.shp'
#pts_inicial <- readOGR(dirname(pts_file),strsplit(basename(pts_file),'[.]')[[1]][1] )

#set.seed(1)
raster_in_pts <- as.data.frame(rasterToPoints(raster_in))
pc_calc <- prcomp(~ raster_in_pts[,1] + raster_in_pts[,2] + raster_in_pts[,3],scale = TRUE)
pc_data <- as.data.frame(pc_calc$x)
pc_kmeans <- cmeans(pc_data,groups,500,verbose=TRUE,method="cmeans")

raster_in_pts <- data.frame(x = raster_in_pts$x,y = raster_in_pts$y,cluster = pc_kmeans$cluster)
cluster_raster <- rasterFromXYZ(raster_in_pts)
cluster_raster@crs <- raster_in@crs
#plot(cluster_raster,col = topo.colors(30))

# Converter falsas zonas de manejo em nodata
cluster_ref <- setValues(cluster_raster,NA)
for(i in 1:groups){
mask_i <- setValues(cluster_raster,NA)
mask_i[cluster_raster == i ] <- i
clump_i <- clump(mask_i,direction = 4)
clump_i_freq <- na.exclude(as.data.frame(freq(clump_i)))
valueID <- clump_i_freq$value[which.max(clump_i_freq$count)]
maxi <- max(clump_i_freq$count)
excludeID <- clump_i_freq$value[which(clump_i_freq$count < maxi)]
clump_i_mask <- clump_i
clump_i_mask[clump_i %in% excludeID] <- NA
clump_i_mask[clump_i == valueID] <- i
cluster_ref[clump_i_mask == i ] <- i
}
#plot(cluster_ref,col = topo.colors(30))

# Função para preencher nodata
func_viz <- function(x){
feq_class <- na.exclude(as.data.frame(table(x)))
if(nrow(feq_class)!=0) {
colnames(feq_class) <- c("Class","Freq")
value_noData <- as.numeric(as.character(feq_class$Class[feq_class$Freq==max(feq_class$Freq)]))
return( value_noData ) }
else {
return( NA ) }
}

# Converter poligono da unidade produtiva em raster para usar como mascara
if(class(pol_mask) == "SpatialpolygonsDataFrame"){
raster_in_copy <- raster_in
extent(raster_in_copy) <- extent(pol_mask)
res(raster_in_copy) <- res(raster_in)[1]
rp <- rasterize(pol_mask,raster_in_copy,1)
raster_in_01 <- !is.na(rp)
}else{
raster_in_01 <- !is.na(raster_in)
}

# Preenchimento de nodata
repeat {
cluster_ref <- extend(cluster_ref,c(1,1))
raster_in_01 <- extend(raster_in_01,1))
cluster_ref_foc <- focal(cluster_ref,w=matrix(1,nrow=3,ncol=3),fun=func_viz,pad=T)
in_s_out_01 <- is.na(cluster_ref_foc)
in_s_out_01[raster_in_01 != 1] <- NA
in_s_out_01[is.na(raster_in_01)] <- NA
#plot(in_s_out_01)
soma_vazio <- sum(in_s_out_01[raster_in_01 == 1],na.rm = TRUE)
print(soma_vazio)
#cluster_ref_foc[raster_in_01 != 1] <- NA
#plot(cluster_ref_foc)
#print(cluster_ref_foc)
if (soma_vazio == 0){
cluster_ref_foc <- focal(cluster_ref_foc,pad=T)
cluster_ref_foc <- focal(cluster_ref_foc,pad=T)
break
} else{
cluster_ref_foc[raster_in_01 != 1] <- NA
#plot(cluster_ref_foc)
cluster_ref <- cluster_ref_foc
}
}

# Converter raster para poligono
pol_clust_ref <- rasterTopolygons(cluster_ref_foc,dissolve=TRUE)
#spplot(pol_clust_ref)
if(Porcentagem_de_vertices_a_serem_mantidos_0_a_100 != 100){
library(rmapshaper)
pol_clust_ref <- ms_simplify(pol_clust_ref,keep = porc_vert)#,keep_shapes = TRUE,explode = TRUE)
pol_clust_ref@data$rmapshaperid <- NULL
}

# Cortar poligono
if(class(pol_mask) == "SpatialpolygonsDataFrame"){
row.names(pol_mask) <- "1"
clip_mask <- gIntersection(pol_clust_ref,pol_mask,byid = T,drop_lower_td = TRUE)
row.names(clip_mask) <- gsub(" 1","",row.names(clip_mask))
keep <- row.names(clip_mask)
clip_mask <- spChFIDs(clip_mask,keep)
pol_clust_ref <- SpatialpolygonsDataFrame(clip_mask,data.frame(layer = pol_clust_ref@data[keep,]))
}else{
raster_in_1 <- raster_in
raster_in_1[!is.na(raster_in_1)] <- 1
pol_1 <- rasterTopolygons(raster_in_1,dissolve=TRUE)
row.names(pol_1) <- "1"
clip_mask <- gIntersection(pol_clust_ref,pol_1,]))
}

# ur.area<-sapply(slot(pol_clust_ref,"polygons"),function(x) sapply(slot(x,"polygons"),slot,"area"))
# mult_part <- unlist(lapply(ur.area,length))
# if(max(mult_part) > 1){
#     single_pol <- multitoone(pol_clust_ref)
#     #single_pol <- SpatialpolygonsDataFrame(single_pol,data.frame(layer = row.names(single_pol)))
#
# }

#caminho minimo

centroids <- gCentroid(pol_clust_ref,byid=TRUE)
centroids_pts <- SpatialPointsDataFrame(centroids,as.data.frame(centroids))
proj4string(centroids_pts) <- pol_clust_ref@proj4string

pts_xy <- as.data.frame(centroids@coords)
pts_xy$id_old <- 1:nrow(pts_xy)
dist_xy <- as.matrix(dist(pts_xy))
adj_matriz <- gtouches(pol_clust_ref,byid = TRUE)
dist_xy[!adj_matriz] <- dist_xy[!adj_matriz] + max(dist_xy) * 1000 # dist_xy[!adj_matriz] + max(dist_xy)
diag(dist_xy) <- 0.0
tsp <- TSP(dist_xy)

if(class(pts_inicial) == "SpatialPointsDataFrame"){
dists_viz <- sqrt((pts_xy$x - pts_inicial@coords[1])^2 + (pts_xy$y - pts_inicial@coords[2])^2)
start_pts_id <- which.min(dists_viz)
} else{
start_pts_id <- which.min(pts_xy$y) # miny
}

ptm <- proc.time()
dist_total <- 1e8
true_n <- 1
Numero_de_Repeticoes <- 1
while(true_n <= Numero_de_Repeticoes) {
tour <- solve_TSP(tsp,method =  "farthest_insertion",rep=20,two_opt =TRUE,two_opt_repetitions = 10,start = start_pts_id)
tour_pts_temp <- pts_xy[tour,]
dist_total_temp <- sum(sqrt(diff(tour_pts_temp$x)^2+diff(tour_pts_temp$y)^2))
if(dist_total_temp < dist_total){
tour_pts <- tour_pts_temp
dist_total <- dist_total_temp
}
true_n <- true_n + 1
print(dist_total)
}
print("Tempo minimizar caminho:")
proc.time() - ptm

# Organizar poligono de saida
pol_resul <- pol_clust_ref[tour_pts$id_old,]
pol_resul@data$layer <- NULL
pol_resul@data$NomeID <- 1:nrow(pts_xy)
#width_str <- length(strsplit(as.character(nrow(pts_xy)),"")[[1]])
width_str <- 3
pol_resul@data$Nome <- paste(iniID,formatC(1:nrow(pts_xy),width = width_str,flag = "0"),sep = "")

mean_r <- extract(raster_in,pol_resul,fun = mean,na.rm = T,small = T)
sd_r <- extract(raster_in,fun = sd,small = T)
cv <- sd_r/mean_r * 100
#len_r <- extract(raster_in,fun=function(x,...)length(x),small = T)
#area_cell <- res(raster_in)[1]^2
#area_r <- len_r * area_cell
#area_total <- sum(area_r)

if(groups > 5) {
apt_kmeans <- kmeans(mean_r,5)$cluster
apt_mean_cluster <- as.data.frame(cbind(mean_r,apt_kmeans))
apt_mean_cluster$ordem <- 1:nrow(apt_mean_cluster)
agreg <- aggregate( V1 ~ apt_kmeans,FUN = mean,data = apt_mean_cluster)
agreg <- agreg[order(agreg$V1),]
agreg$aptidao <- c("MUITO BAIXA","BAIXA","MEDIA","ALTA","MUITO ALTA")
agreg_merg <- merge(apt_mean_cluster,agreg,'apt_kmeans',sort = FALSE)
ord <- agreg_merg$ordem
field_apt <- agreg_merg[order(ord),]$aptidao
} else if(groups > 3){
apt_kmeans <- kmeans(mean_r,3)$cluster
apt_mean_cluster <- as.data.frame(cbind(mean_r,]
agreg$aptidao <- c("BAIXA","ALTA")
agreg_merg <- merge(apt_mean_cluster,]$aptidao

}else{
field_apt <- "MEDIA"
}

pol_resul@data$Media <- as.vector(mean_r)
pol_resul@data$Desv_Pad <- as.vector(sd_r)
pol_resul@data$CV <- as.vector(cv)
pol_resul@data$Area_ha <- area(pol_resul) / 10000
pol_resul@data$Aptidao <- field_apt
pol_resul@data$Constante <- 1

Zonas_de_Manejo =  pol_resul

报告的错误如下：

 Warning message:
    In if (class(pol_mask) == "SpatialpolygonsDataFrame") { :
    the condition has length> 1 and only the first element will be used Error in loadNamespace (i,c 
    (lib.loc,.libPaths ()),versionCheck = vI [[i]])
    vI object not found
    Calls: ms_simplify ... tryCatch -> tryCatchList -> tryCatchOne -> <Anonymous>
    Execution interrupted