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

Fri Nov 5 15:45:05 CET 2010

On Fri, 5 Nov 2010, Rodrigo Aluizio wrote:

> 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.

The SpatialPolygons object will be exported with the "hole" statuses set 
correctly to the shapefile. If you are plotting externally, you depend on 
the facilities there. In R, you are using 2.12.0, which is the first 
release to handle polygon holes repectfully. Try:

plot(as(PCSPH, "Spatial"), axes=TRUE, bg="azure")
plot(PCSPH, add=TRUE, col="grey", usePolypath=TRUE)

to "see through" the holes. You are right that with usePolypath=FALSE, and 
in earlier versions of R, inner ring holes were painted over the filled 
outer rings with a background colour. We probably should make usePolypath 
TRUE by default now, and check the R version in use.

Roger

>
> 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
>
>

-- 
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