[R-sig-Geo] Intersection between two SpatialPolygonsDataFrame using gIntersection function in R

[Roger Bivand]

On Mon, 17 Dec 2012, Sophie Padié wrote:

> Dear subscribers,
>
> After a quick search on the net and on this mailing list, I didn't find any
> solution, so I hope my question won't be that trivial.
>
> I want to select a subset of a map (a SpatialPolygonsDataFrame) to create a
> new vector map = a part of the initial map included inside a given polygon
> (eventually, inside several polygons). I have a large study area, under the
> object "parcel", and would like to select the part of the map inside the
> polygon "hr_kud". This is basically the part of my study area used by an
> animal, so I want to compute basic statistics on this new map and to use it
> for new analyses (of habitat selection). I have to repeat this operation
> for about 100 animals...
>
> Here are my data : http://dl.free.fr/vX4TA8nsE
>
> *parcel* = my study area map. and *hr_kud* = an example of homerange used
> to reduce the study area
>
> Here is a part of my R code, used under Rstudio (R version 2.15.2
> (2012-10-26)) :

Never run examples under any IDE, it may only obscure, made worse here by 
posting in HTML (which is not allowed).

>
>> library(rgeos)> load("Data.RData")
>

The problem arises earlier than this. If I tabulate the number of 
characters in the OGC SFS comments of the Polygons objects, I see:

res <- sapply(slot(parcel, "polygons"), comment)
any(is.null(res))

# [1] FALSE

table(sapply(res, nchar))

#    0     1     3     5     7     9    11    13
#   31 19673   299    49    12     2     2     1

and:

parcel1 <- parcel[-odd,]
gIsValid(parcel1)

# Error in createPolygonsComment(p) :
#  rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon 
#  for hole at index 3

So there are problems in the creation of parcel, and these are leading to 
the crash in rgeos (empty assigned comment strings; they are NULL if not 
assigned). So you need to specify the antecedents of parcel - as well, of 
course, as reporting the output of sessionInfo() - including all the 
packages used for creating the damaged object, and the complete startup 
message of rgeos (and rgdal if implicated).

I'll add code to guard against attempts to use empty comment strings, but 
it would be very helpful to establish where they came from (rgdal?).

In addition to these issues, parcel has further topology issues.

A fix is to use the hole analysis wrapper in maptools, which uses rgeos to 
try to correct wrong hole assignments, also correcting the damaged 
comments:

library(maptools)
pls <- slot(parcel, "polygons")
pls1 <- sapply(pls, checkPolygonsHoles)
slot(parcel, "polygons") <- pls1
res <- sapply(slot(parcel, "polygons"), comment)
table(sapply(res, nchar))

#     1     3     5     7     9    11    13
# 19701   302    49    12     2     2     1

gIsValid(parcel)
# [1] TRUE
hr_map <- gIntersection(parcel, hr_kud, byid=TRUE)
# Warning message:
# In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") :
#   spgeom1 and spgeom2 have different proj4 strings

which runs to completion with no problems; you could remove the warning in 
the obvious way, the CRS are:

[1] " +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666667 
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356514.999904194 
+units=m +no_defs"

[1] " +init=epsg:27572 +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 
+towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"

which are almost the same, but are not string-identical.

Hope this helps,

Roger

>
>> class(parcel)[1] "SpatialPolygonsDataFrame"
> attr(,"package")
> [1] "sp"
>
>> head(parcel)  AREAHA RECNO      AREA AREASQKM        OS_06 GRD_CAT_06
> 0 0.0519     0  519.1941        0 chemin prive        way
> 1 0.0127     1  127.0963        0 chemin prive        way
> 2 0.0596     2  596.2204        0 chemin prive        way
> 3 0.0248     3  248.0452        0    route D90        way
> 4 0.0375     4  375.0141        0 route privee        way
> 5 0.1260     5 1260.3941        0 chemin prive        way
>
>> dim(parcel) [1] 20069 7 > class(hr_kud) [1] "SpatialPolygonsDataFrame"
> attr(,"package") [1] "sp" > object.size(hr_kud) 9536 bytes >
> object.size(parcel)
> 88250824 bytes
>
>
>
> And when I run my final command, it causes Rstudio to crash :
>
>> hr_map <- gIntersection(parcel, hr_kud, byid=TRUE)
>
>
>
> It is not the first crash reported, but none of the answers find on the net
> are efficient for me... It's exactly the same using R or Rstudio and all my
> packages are updated.
>
> So If you have any idea of why it's crashing, or of an alternative method
> to create my subset map or to intersect two SpatialPolygonsDataFrame in R,
> It will be really appreciated !
>
> thanks in advance,
> Sophie
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
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