Fwd: Representação de gráfico tipo torta em um shapefile
Boa tarde Walmes, Estou tentando postar minha dúvida abaixo e aparece um erro no servidor, saberia me dizer o porque? Obrigado, Alexandre -------- Mensagem original -------- Assunto: Representação de gráfico tipo torta em um shapefile Data: Sat, 24 May 2014 11:24:58 -0400 De: ASANTOS <alexandresantosbr@yahoo.com.br> Para: r-br@listas.c3sl.ufpr.br Bom dia Pessoal, Estou querendo inserir um gráfico do tipo torta em um shapefile de bacias hidrográficas da América do Sul (bacias), fazendo a inserção de um gráfico torta para cada bacia , onde localizam-se meus pontos de interesse (pontos$LEVEL2), porém os valores do gráfico gostaria que fosse a frequência de pontos$CT em cada pontos$LEVEL2, mas estou tento dificuldades e vou explicar passo a passo: # Start ------------------------------------------------------------------------------------------------------ #Pacotes require (maptools) require(shapefiles) require(rgdal) require(mapplots) require(RColorBrewer) # #Shapefile e dados 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") tokens <- gsub("^.*/s/","",dirname(links)) fileNames <- basename(links) newLinks <- file.path("http://dl.dropbox.com/s", tokens, fileNames); newLinks for (a in newLinks) { tryCatch(download.file(a, dest=basename(a), mode='wb'), error=function(...) print("Falha no download!"))} # #Abrindo o shapefile # CRS.new <- CRS("+proj=longlat +datum=WGS84") bacias <- readOGR(".", "sa_bas_ll_r500m") proj4string(bacias) <- CRS.new plot(bacias) # # #Coordenadas dos pontos + atributos pontos<-read.csv("indv_atributos.csv", sep=";",h=T) coordinates(pontos) <- c("x_long_dec","y_lat_dec") proj4string(pontos) <- CRS.new # # Selecionos os pontos que estão dentro de cada bacia hidrográfica inside.sa <- !is.na(over(pontos, as(bacias, "SpatialPolygons")))## inside.sa # # Pontos pertencentes a cada bacia nível 2, tenho varios niveis de classificação e escolhi o LEVEL2 bacias@bbox pontos@bbox pontos$LEVEL2 <- over(pontos, bacias)$LEVEL2 ## Adiciono o gráfico torta para cada pontos$LEVEL2 a frequencia de pontos$CT add.pie(z=pontos$CT, x=pontos$x_long_dec, y=pontos$y_lat_dec, radius=1, col=c(alpha("orange", 0.6), alpha("blue", 0.6), alpha("black", 0.6), alpha("red", 0.6), alpha("grey", 0.6)), labels=" ") #END---------------------------------------------------------------------------------------------------------- Mas não obtive sucesso, alguém poderia me dar um help? Obrigado, -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Alexandre, boa tarde! Sugiro que você trabalhe com o mapa de bacias "dissolvidas" para o nível 2. Para isso você pode fazer uso de maptools::unionSpatialPolygons(). Para fins de visualização você pode ainda utilizar outras funções de simplificação. Feito isso, além de visualizar mais rápido, você poderá obter o centróide dos polígonos dissolvidos pra locar os gráficos de torta. Depois de rodar o sp::over() você precisa tabular as frequências por bacia e adicionar a coordenada de cada uma. A entrada do add.pie() é individual, então você poderá pensar num laço ou uma função da família apply pra entradas múltiplas. Atte., Éder Comunello <c <comunello.eder@gmail.com>omunello.eder@gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]
Boa noite Éder, Tentei seguir seu roteiro usando maptools::unionSpatialPolygons() e não deu certo a dissolução de bacias, fiz: require(shapefiles) library(maptools) library(gpclib) library(rgdal) library(PBSmapping) #Shapefile e dados links <- c( "https://www.dropbox.com/s/3ph2630xicpexo0/dem.crop.final.tif", "https://www.dropbox.com/s/3ph2630xicpexo0/dem.crop.final.tif", "https://www.dropbox.com/s/3ph2630xicpexo0/dem.crop.final.tif", "https://www.dropbox.com/s/3ph2630xicpexo0/dem.crop.final.tif") if (sum(sapply(basename(links), file.exists))<4) { lapply(links, function(a) tryCatch(download.file(a, dest=basename(a), mode='wb'), error=function(...) message("ERRO!")))} else message('OK!') # #abrido shapefile # bacias <- readOGR(".", "sa_bas_ll_r500m") CRS.new <- CRS("+proj=longlat +datum=WGS84") proj4string(bacias) <- CRS.new plot(bacias) str(bacias) ## Separando os ID dos atributos do LEVEL2, pois tenho 6 níveis de bacias IDlevel2<- cut(bacias@data$LEVEL2, range(bacias@data$LEVEL2), include.lowest=TRUE) ## Bacias "dissolvidas" para o nível 2 Dissolve_bacia <- unionSpatialPolygons(bacias ,Dissolve_bacia) Dissolve_bacia_2 <- SpatialPolygons2PolySet(DissolveResult) ## Representação gráfica plotPolys(Dissolve_bacia_2 , proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude") E a figura resultante em plotPolys não são minhas bacias no nível 2. Poderia me dizer o que estou fazendo de errado, Obrigado, -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ====================================================================== Em 26/05/2014 15:12, Éder Comunello escreveu: Alexandre, boa tarde! Sugiro que você trabalhe com o mapa de bacias "dissolvidas" para o nível 2. Para isso você pode fazer uso de maptools::unionSpatialPolygons(). Para fins de visualização você pode ainda utilizar outras funções de simplificação. Feito isso, além de visualizar mais rápido, você poderá obter o centróide dos polígonos dissolvidos pra locar os gráficos de torta. Depois de rodar o sp::over() você precisa tabular as frequências por bacia e adicionar a coordenada de cada uma. A entrada do add.pie() é individual, então você poderá pensar num laço ou uma função da família apply pra entradas múltiplas. Atte., Éder Comunello <c <mailto:comunello.eder@gmail.com>omunello.eder@gmail.com <mailto:omunello.eder@gmail.com>> Dourados, MS - [22 16.5'S, 54 49'W]
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]
Boa tarde Éder, GENIAL!!!!! Era isso mesmo que eu precisava, muito clara e limpa sua programação, tinha conseguido achar o erro e usar bacias@data$LEVEL2, mas na hora de desenhar as tortas, estava tentando calcular os centroides das bacias com trueCentroids = gCentroid(bacias2,byid=TRUE) do pacote GISTools e convertendo para um novo shape e depois novamente para polígonos espaciais, tava ficando uma gambiarra para desenhar as tortas, fico com a sua opção, Muito Obrigado!!! Alexandre Em 30/05/2014 08:14, Éder Comunello escreveu:
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 <http://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 <http://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 <mailto:comunello.eder@gmail.com>omunello.eder@gmail.com <mailto:omunello.eder@gmail.com>> Dourados, MS - [22 16.5'S, 54 49'W]
_______________________________________________ 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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
participantes (2)
-
ASANTOS -
Éder Comunello