Alexandre, bom dia!
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...
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.
### <code r>
setwd("C:/LAB/RGIS/maps/bacias")
pkgs <-
c("maptools", "shapefiles", "rgdal", "mapplots",
"RColorBrewer")
sapply(pkgs, require,
character=T)
### Download
----------------------------------------------------------------
#
links <- c(
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>
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.