Boa tarde colegas listeiros, como vão? Em relação a minha ultima dúvida, consegui bolar uma lógica que me retornasse o que eu estava buscando no e-mail encaminhado acima e gostaria de compartilha-la com vocês. Logo após fazer a clusterização e gerar uma lista com valores de cluster, adicionei essa camada de informação diretamente ao meu shapefile de pontos original:
shapefile_dados_das_camadas_utilizadas@dados$Cluster = lista_clusteres
Em seguida, gerei um arquivo raster com as informações de cluster e de dimensões 20x20 m:
modelo_interpolador = gstat(formula = shapefile_dados_das_camadas_utilizadas$Cluster ~ 1, locations = shapefile_dados_das_camadas_utilizadas, nmax = 1, set = list(idp = 0)) raster_Cluster = raster(predict( modelo_interpolador, as(celulas_dentro_do_talhao, 'SpatialPixels')))
Por fim, transformei o arquivo raster em um shape de poligonos e passei TRUE no argumento dissolve, para unir os polígonos de acordo com o ID 'Cluster'
shapefile_poligonos = rasterToPolygons(raster_Cluster , n = 4, dissolve = TRUE)
E essa foi a forma com que contornei o problema proposto no inicio da mensagem. Caso alguém tenha bolado um outro caminho e quiser compartilhar, tenho muito interesse em ver logicas diferentes de se trabalhar com dados espacializados no R! Aproveito para deixar aqui uma dúvida que surgiu ao final dessa tarefa. Seria possível suavizar os contornos dos polígonos gerados, assim como no ArcGis, ao utilizar as funções spline ou smooth, sem que buracos e/ou sobreposições de polígonos ocorram? No R, ao utilizar a função smooth(), do pacote smoothr, acontece sempre de aparecerem buracos entre polígonos e também sobreposições de polígonos no meu resultado final, o que não é desejado. Mais uma vez, agradeço pela atenção de todos! Yury Duarte Engenheiro Agrônomo - ESALQ/USP Em seg, 3 de set de 2018 às 14:20, Yury Duarte <yurynepomuceno@gmail.com> escreveu:
Boa tarde colegas listeiros, como vão?
Há algum tempo atras recorri a lista com algumas duvidas de manipulação do SIG com o R e me levaram bem longe, mas ainda sigo com algumas dificuldades. Estou fazendo uma analise de agrupamento para um local com base em informações georreferenciadas e espacializadas numa grade de 20x20 metros. Inicialmente tenho essas informações num shape de pontos e também tenho o shape do contorno da minha área. O resultado da minha analise de agrupamento é uma lista indicando o cluster de cada conjunto lat;long. Após esse passo, gostaria de criar um filtro que fosse capaz de identificar e substituir o valor de um pixel caso ele estiver cercado por valores de cluster diferente. existe alguma ferramenta que faça algo do tipo?
Paralelo a isso, gerei um shape de polígonos com a mesma dimensão da grade de pontos (20x20 metros). Nesse momento ocorre a quebra do script, pois o número de polígonos gerado é menor do que o número de informações que eu tenho no meu shape de pontos original e na minha lista de clusters. Minha ideia foi de fazer um merge (ou coisa do tipo) entre minha lista de cluster e uma lista com os centroides de cada polígono, assim, eu eliminaria os pontos "excedentes" em relação aos polígonos, mas não consegui fazer com que funcionasse. Segue script ao final da mensagem.
Desde já, agradeço pela ajuda de todos!
library(spatstat) library(maptools) library(shapefiles) library(foreign) library(magrittr) library(formattable) library(e1071) library(rgeos) library(rgdal) library(rgeos) library(raster) library(sp) library(RcppCNPy) library(deldir) library(dismo) library(dplyr) library(ggplot2) library(gstat)
#------Ler shape de contorno e shape de dados + transformar coordenadas em UTM------# #-----Lendo shape de dados------# shapefile_dados_das_camadas_utilizadas = readOGR(paste0('diretorio', 'dados.shp'))
#-----Lendo shape de contorno------# shapefile_contorno_do_talhao = readOGR(paste0('diretorio', 'contorno.shp'))
#------Clusterizacao------# shapefile_dados_padronizados = scale(shapefile_dados_das_camadas_utilizadas@data) shapefile_dados_geolocalizacao = shapefile_dados_das_camadas_utilizadas@coords wss = (nrow(shapefile_dados_padronizados)-1)*sum(apply(shapefile_dados_padronizados,2,var)) cluster_maximo = 30 for (i in 1:cluster_maximo) wss[i] = sum(kmeans(shapefile_dados_padronizados, centers=i)$withinss) plot(1:cluster_maximo, wss, type="b", xlab="Número of Clusters", ylab="Soma dos quadrados dentro dos clusteres")
#dendograma = shapefile_dados_padronizados %>% dist %>% hclust #plot(dendograma) #rect.hclust(dendograma, k = 3, border = "blue") #rect.hclust(dendograma, k = 4, border = "red")
escolher_num_cluster = 6 lista_clusteres = data.frame((kmeans(shapefile_dados_padronizados, centers = escolher_num_cluster)$cluster),shapefile_dados_geolocalizacao)
#------Transformando em poligono------# x = shapefile_dados_das_camadas_utilizadas@coords[,1] y = shapefile_dados_das_camadas_utilizadas@coords[,2] df_xy = data.frame(shapefile_dados_das_camadas_utilizadas@coords [,1],shapefile_dados_das_camadas_utilizadas@coords[,2]) #coordinates(df_xy) <- ~x+y z = deldir(x,y,plot=TRUE,wl='triang',wp="none") w = tile.list(z) polys = vector(mode="list", length=length(w)) for (i in seq(along=polys)) { pcrds = cbind(w[[i]]$x, w[[i]]$y) pcrds = rbind(pcrds, pcrds[1,]) polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i)) } SP = SpatialPolygons(polys) plot(SP) celulas_dentro_do_talhao = raster::intersect(SP, shapefile_contorno_do_talhao) plot(celulas_dentro_do_talhao)
#------Adicionando Informacao de cluster ao poligono------# length(celulas_dentro_do_talhao) confirmando_pontos_dentro_da_area = data.frame(coordinates(celulas_dentro_do_talhao)) length(lista_clusteres[,1]) colnames(lista_clusteres) = c('Cluster', 'Long', 'Lat') colnames(confirmando_pontos_dentro_da_area) = c('Long', 'Latitude') #confirmando_pontos_dentro_da_area$Long = round(confirmando_pontos_dentro_da_area$Long, 5) #lista_clusteres$Long = round(lista_clusteres$Long, 5) #teste = merge(confirmando_pontos_dentro_da_area, lista_final, all.x = TRUE) #celulas_dentro_do_talhao@data$Cluster = '???' plot(celulas_dentro_do_talhao) writeOGR(celulas_dentro_do_talhao, 'C:\\Users\\Yury\\Desktop\\InCeres\\UGD', layer = 'sss', driver = 'ESRI Shapefile')
Yury Duarte Engenheiro Agrônomo - ESALQ/USP