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

Fri Nov 5 14:42:13 CET 2010

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