Alexandre, bom dia! Basicamente o erro está no vetor de identificadores que você definiu para unionSpatialPolygons(). Ao invés de 'Dissolve_bacia', bastaria entrar com " bacias@data$LEVEL2". Resolvi rodar o procedimento todo e obtive a figura abaixo. Obtive 'tortas' sem subdivisão, provavelmente porque tem bacias com a ocorrência de apenas um valor do atributo $CT. Disponibilizo o código que produzi e se for de ajuda, peço que disponibilize o produto final pra vermos como ficou. Comentei na medida do possível... ### <code r> setwd("C:/LAB/RGIS/maps/bacias") pkgs <- c("maptools", "shapefiles", "rgdal", "mapplots", "RColorBrewer") sapply(pkgs, require, character=T) ### Download ---------------------------------------------------------------- # links <- c( "https://www.dropbox.com/s/enrdm739ts4eba0/sa_bas_ll_r500m.shp", "https://www.dropbox.com/s/fnnumaix0v88jo1/sa_bas_ll_r500m.shx", "https://www.dropbox.com/s/o152vc0n15w54em/sa_bas_ll_r500m.dbf", "https://www.dropbox.com/s/2l24sso20645511/indv_atributos.csv") newLinks <- gsub("^.*/s/", "http://dl.dropbox.com/s/", links); newLinks vDown <- sum(sapply(newLinks, function(x) file.exists(basename(x)))) if (vDown<4) { sapply(newLinks, function(a) { tryCatch(download.file(a, dest=basename(a), mode='wb'), error=function(...) print("Falha no download!"))})} sapply(basename(newLinks), file.info)[1,] ### sizes ### Carregando shapefile ---------------------------------------------------- # CRS.new <- CRS("+proj=longlat +datum=WGS84") bacias.ori <- readOGR(".", "sa_bas_ll_r500m") proj4string(bacias.ori) <- CRS.new ### Visualização excessivamente demorada! Melhor utilizar o vetor simplificado! ### Os tempos anotados (s) foram tomados num netbook de configuração modesta. # system.time(plot(bacias.ori, col=topo.colors(32), border=NA)) ### 44s ### Entendendo os dados ----------------------------------------------------- # t(names(bacias.ori@data)) ### LEVELS c(5:10) sapply(bacias.ori@data[,c(3:10)], function(x) length(unique(x))) # SA_BAS_ SA_BAS_ID LEVEL1 LEVEL2 LEVEL3 LEVEL4 LEVEL5 LEVEL6 # 5339 5339 11 92 569 1775 2374 2430 # ^ ^ ^ ^ ^ ^ ^ ^ # | | | | | | | | # Elementos Elementos MegaBacias Bacias Sub-Bacias Mini-Bacias Micro-Bacias Nano-Bacias # Lembrando que a escala acima é só uma brincadeira... :) # # São 5339 elementos (possivelmente todos polígonos) que podem ser retratados como 11 Grandes # Bacias ou 92 Bacias. # Quando você usa o cut() pra atribuir cor, na verdade o R desenha os 5339 elementos e # 'pinta' os que tem o mesmo rótulo com a mesma cor. Usando o unionSpatialPolygons() com # $LEVEL2 você passará a ter 92 polígonos. # # Na hora de locar as 'tortas' é outro problema, pois cada bacia do LEVEL2 é representada por # diversos polígonos, por vezes centenas, e você teria que escolher um desses polígonos ou # # obter a coordenda média desses. ### Simplificando ----------------------------------------------------------- # bacias.mod <- unionSpatialPolygons(bacias.ori, bacias.ori$LEVEL2) ### modIFICADO bacias <- rgeos::gSimplify(bacias.mod, .4, topologyPreserve=TRUE) # .4 degrees length(bacias.ori); length(bacias.mod); length(bacias) # system.time(plot(bacias.ori, col=topo.colors(32), border=NA)) ### 44s # system.time(plot(bacias.mod, col=topo.colors(32), border=NA)) ### 30s system.time(plot(bacias, col=topo.colors(32), border=NA)) ### 18s ### Info -------------------------------------------------------------------- # str(bacias, max=2) str(bacias@data) ### "SpatialPolygons" - Não tem @data!!! ### Podemos utilizar o @data original (bacias.ori) str(bacias.ori@data) ### "SpatialPolygonsDataFrame" - DATA!!! ### Continuaremos a utilizar bacia.ori para quantificações! ### Remover Polígonos menores ----------------------------------------------- # ### Esse trecho é para eliminar polígonos menores para fins de agilizar a visualização ### Deixo apenas pra referência. pois não teve muito efeito # subPolys <-sapply(bacias@polygons, function(a) sum(sapply(a@Polygons, length))) # table(subPolys) # subAreas <- lapply(bacias@polygons, function(a) sapply(a@Polygons, function(b) b@area)) # sapply(subAreas, sum) # sum(sapply(subAreas, sum)<2) # str(bacias, max.level=3) # bigPolys <- which(sapply(subAreas, sum) >= 2); length(bigPolys) # BIG <- lapply(bigPolys, function(a) {bacias@polygons[[a]]}); length(BIG) # bacias2 <- SpatialPolygons(BIG) # system.time(plot(bacias2, col=topo.colors(32), border=NA)) ### 18s ### Coordenadas dos pontos + atributos pontos <- read.csv("indv_atributos.csv", sep=";", head=T) names(pontos) <- c("num_loc", "LAT", "LON", "CT") coordinates(pontos) <- ~LON+LAT proj4string(pontos) <- CRS.new ### Pontos no mapa system.time(plot(bacias, col=topo.colors(32), border=NA)) ### 18s points(pontos) ### Pontos pertencentes a cada bacia nível 2 (LEVEL2) bacias@bbox; pontos@bbox id2 <- over(pontos, bacias.ori)$LEVEL2 id2i <- id2[!is.na(id2)] ### Inside Points; if outside == NA newData <- reshape2::dcast(pontos@data, id2i~CT, length, value.var='CT'); newData coords <- coordinates(bacias) coords.df <- data.frame(id2i=as.numeric(row.names(coords)), coords) newData2 <- merge(newData, coords.df); newData2 ### Observe que em vários casos a torta terá apenas um valor!!! ### Adiciono o gráfico torta para cada pontos$LEVEL2 a frequencia de pontos$CT x11() system.time(plot(bacias, col=topo.colors(32), border=NA, axes=T)) ### 18s lapply(1:nrow(newData2), function(x) add.pie(z=as.numeric(newData[x,2:7]), x=newData2$X1[x], y=newData2$X2[x], radius=3, labels=" ")) ### </code> Éder Comunello <c <comunello.eder@gmail.com>omunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
