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
