[R-sig-Geo] Defining Holes Inside a Polygon

Fri Nov 5 15:23:18 CET 2010

Thank you very much Roger. But , unfortunately it's was almost there. Like I
told before, I'll plot this Shapefiles over a raster map, and even if I set
the pbg='transparent' the color of the layer below (another polygon) will
still appear over my raster image. That's why I'm trying to turn these
'islands' (that will be represented by the background raster image) into
'holes' of my new shapefile.

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.Bivand at 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.Bivand at 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.aluizio at 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.Bivand at 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.Bivand at nhh.no