Anyway, thank you very much for the effort. I won't be able to figure
something like out that so fast.
Rodrigo.
-----Mensagem original-----
De: Roger Bivand [mailto:roger.biv...@nhh.no]
Enviada em: sexta-feira, 5 de novembro de 2010 11:42
Para: Rodrigo Aluizio
Cc: 'R Help'
Assunto: Re: RES: [R-sig-Geo] Defining Holes Inside a Polygon
On Fri, 5 Nov 2010, Rodrigo Aluizio wrote:
Very sorry about the attached file, I'm not used to mailing lists,
won't happen again.
Well, I've already tried the checkPolygonsHoles function as specified
in a book of Yours, and even then nothing changes in the final
Shapefiles, all the polygons remain with FALSE "hole" slots.
Below are the code and sessionInfo().
# Required Packages
library(rgdal)
library(RODBC)
library(maptools)
gpclibPermit()
#Loading Data
Data<-odbcConnectExcel2007('Sub-Provincias-Forams-Dados.xlsx',readOnly
=T)
PCdata<-sqlFetch(Data,'PC',rownames='ID')
odbcCloseAll()
# Loading Lines
PC<-MapGen2SL('PCForams.dat',proj4string=CRS('+proj=longlat
+ellps=WGS84'))
# Verifing closed Lines
PCP<-sapply(slot(PC,'lines'),function (x) {
crds<-slot(slot(x,'Lines')[[1]],'coords')
identical(crds[1,],crds[nrow(crds),])
})
# Preparing the Poligons
PCP2<-PC[PCP==T]
lista_de_pols<-slot(PCP2,'lines')
PCSP<-SpatialPolygons(lapply(lista_de_pols, function(x) {
Polygons(list(Polygon(slot(slot(x,'Lines')[[1]],
'coords'))),ID=slot(x,'ID'))
}),proj4string=CRS('+proj=longlat +ellps=WGS84'))
The problem is that you had 39 Lines objects, and now have 39 Polygons
objects, each with a single outer ring by definition. The
checkPolygonsHoles() function checks multiple rings in a Polygons object.
What you need to do from here will depend on whether you know which inner
rings belong to which outer rings or not. If you do not, you could try
something like:
crds <- coordinates(PCSP)
pips <- lapply(slot(PCSP, "polygons"), function(x)
sp:::pointsInPolygons(crds, x))
ipips <-sapply(pips, which)
which(sapply(ipips, length) > 1)
nipips <- which(sapply(ipips, length) > 1) npips <- ipips[nipips] ors <-
integer(length(slot(PCSP, "polygons"))) ors[npips[[1]]] <- 1 ors[npips[[2]]]
<- 2 ors ors[2] <- 3 ors[3] <- 4 ors
to get a candidate vector by seeing which rings contain more than one
centroid.
In one case, we can use unionSpatialPolygons() to combine three external
"islands", which are not holes, and in two cases there are no multi-rings:
pls <- vector(mode="list", length=4)
oo2 <- PCSP[which(ors == 2)]
gpclibPermit()
pls[2] <- slot(unionSpatialPolygons(oo2, rep("2", 3)), "polygons")[[1]]
pls[3] <- slot(PCSP[2], "polygons")[[1]] pls[4] <- slot(PCSP[3],
"polygons")[[1]]
The final collection is harder:
oo1 <- PCSP[which(ors == 1)]
oo1pls <- slot(oo1, "polygons")
areas <- sapply(oo1pls, function(x) sapply(slot(x, "Polygons"), slot,
"area"))
marea <- which.max(areas)
o1lst <- vector(mode="list", length=length(oo1pls)) for (i in
1:length(oo1pls)) {
oi <- slot(slot(oo1pls[[i]], "Polygons")[[1]], "coords")
hole <- ifelse(i == marea, FALSE, TRUE)
o1lst[[i]] <- Polygon(oi, hole=hole)
}
pls[1] <- Polygons(o1lst, ID="1")
PCSPH <- SpatialPolygons(pls)
plot(PCSPH, axes=TRUE, col="grey", pbg="white")
now works. You'll need to revisit the data, however, as you only have four
Polygons objects now, so will need a PCdata object with 4 rows. You'll
probably also need to set the row.names of both PCSPH and PCdata to match
exactly too, as they may no longer do so. You'll also need to make sure that
the proj4string is properly set - I've ignored it here.
Hope this helps,
Roger
# Defining the Holes
pls_out<-lapply(slot(PCSP,"polygons"),checkPolygonsHoles)
PCSPH<-SpatialPolygons(pls_out,proj4string=CRS('+proj=longlat
+ellps=WGS84'))
PCSPDF<-SpatialPolygonsDataFrame(PCSPH,PCdata)
writeOGR(PCSPDF,dsn='.',layer='PCForams',driver='ESRI Shapefile')
# Session Info
R version 2.12.0 (2010-10-15)
Platform: x86_64-pc-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=Portuguese_Brazil.1252
[2] LC_CTYPE=Portuguese_Brazil.1252
[3] LC_MONETARY=Portuguese_Brazil.1252
[4] LC_NUMERIC=C
[5] LC_TIME=Portuguese_Brazil.1252
attached base packages:
[1] grDevices datasets splines graphics stats
[6] tcltk utils methods base
other attached packages:
[1] gpclib_1.5-1 maptools_0.7-38 lattice_0.19-13
[4] foreign_0.8-41 RODBC_1.3-2 rgdal_0.6-28
[7] sp_0.9-72 svSocket_0.9-50 TinnR_1.0.3
[10] R2HTML_2.2 Hmisc_3.8-3 survival_2.35-8
loaded via a namespace (and not attached):
[1] cluster_1.13.1 grid_2.12.0 svMisc_0.9-60
[4] tools_2.12.0
Thanks.
Rodrigo.
-----Mensagem original-----
De: Roger Bivand [mailto:roger.biv...@nhh.no] Enviada em: sexta-feira,
5 de novembro de 2010 08:01
Para: Rodrigo Aluizio
Cc: R Help
Assunto: Re: [R-sig-Geo] Defining Holes Inside a Polygon
On Fri, 5 Nov 2010, Rodrigo Aluizio wrote:
Hi List.
I?m trying to solve an issue here, but it?s getting tricky.
I have an MapGen .dat file (attached) that contains about 39 closed
lines.
Three of them (the largest ones) will be my main polygons and the
rest of them must become holes inside one of the three major areas.
I?m able to transform all the closed lines into SpatialPolygons and
even create shapefiles with them, but how can I modify the ?hole?
slot of the smaller polygons so I can keep them transparent when
overlaying another map (those holes will overlay islands of a raster
image)?
Next time avoid sending 200K to over 1700 people, please, we have a
planet to protect! Put any desired attachments on a website and provide a
link.
You should also have provided the code (in your message) to get from:
library(maptools)
l1 <- MapGen2SL("PCForams.dat")
to your problem - as those choices may be the cause of the problem;
you have not provided the output of sessionInfo() either. Once things
are lined up, you should - if you may - say:
gpclibPermit() # very restrictive license pls_out <-
lapply(slot(SpPols0, "polygons"), checkPolygonsHoles)
SpPols1 <- SpatialPolygons(pls_out)
Or if you can install rgeos from source on R-Forge, equivalently:
library(rgeos)
SpPols1 <- createSPComment(SpPols0)
Hope this helps,
Roger
Thank you for the attention and help.
Regards.
-------------------------------------------------------------------
MSc. <mailto:r.alui...@gmail.com> Rodrigo Aluizio
Centro de Estudos do Mar/UFPR
Laboratório de Foraminíferos e Micropaleontologia Ambiental Avenida
Beira Mar s/n - CEP 83255-971 Pontal do Paraná - PR - Brazil
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School
of Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: roger.biv...@nhh.no
--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: roger.biv...@nhh.no
