[R-sig-Geo] Overlaying SpatialLines by SpatialPolygons

[Roger Bivand]

On Thu, 9 Jun 2011, Mathieu Rajerison wrote:

> Hi,
>
>
> Thanks a lot for taking my problem into account.
>
> I've put a data subset as attached file. It contains two files: roads.shp
> and cities.shp.
>
> routes is the result of an intersection between initial complete roads file
> and cities. Features in actual roads have been split up, therefore gContains
> test can directly be executed.

This works, but could be tidied up in a loop:

library(rgeos)
library(rgdal)
cities <- readOGR(".", "cities")
routes <- readOGR(".", "routes")
candidates <- gIntersects(routes, cities, byid=TRUE)
i1 <- gIntersection(routes[candidates[1,],], cities[1,])
i2 <- gIntersection(routes[candidates[2,],], cities[2,])
i3 <- gIntersection(routes[candidates[3,],], cities[3,])
ii <- rbind(i1, i2, i3, makeUniqueIDs=TRUE)
row.names(ii) <- row.names(cities)
SpatialLinesLengths(ii)

With a loop:

library(rgdal)
cities <- readOGR(".", "cities")
routes <- readOGR(".", "routes")
candidates <- gIntersects(routes, cities, byid=TRUE)
nc <- nrow(candidates)
out <- vector(mode="list", length=nc)
for (i in seq(along=out)) {
   out[[i]] <- gIntersection(routes[candidates[i,],], cities[i,])
   row.names(out[[i]]) <- as.character(i)
}
ii <- do.call("rbind", out)
row.names(ii) <- row.names(cities)
SpatialLinesLengths(ii)

This gives the total lengths by city that you need. Using a candidate 
matrix avoids trying to do intersections on objects that do not intersect.

Hope this helps,

Roger

>
> In roads, a column named "contained" indicates whether the road is
> considered contained or not after my gContains test
>
>
> Mathieu
>
> 2011/6/8 Roger Bivand <Roger.Bivand at nhh.no>
>
>> On Wed, 8 Jun 2011, Mathieu Rajerison wrote:
>>
>>  Hi,
>>>
>>>
>>> First, thanks to all. It's the first time I use rgeos.
>>>
>>> I finally managed to get a layer of roads (routes.bdc object) that are the
>>> result of the intersection between roads and cities (bdc object).
>>>
>>> But the problem now is that some intersected roads are not considered
>>> being
>>> contained in any of the cities, although they should be, as the result of
>>> an
>>> intersection. I put two images of them as attached files (city contours in
>>> light orange, roads in blue, and some selected roads in yellow)
>>>
>>> http://hpics.li/750b6bf
>>> http://hpics.li/578ab3d
>>>
>>
>> What is needed is a reproducible example, so that the classes of the output
>> objects become evident. Can you provide such an example with publically
>> available data, or put a relevant data subset on a link?
>>
>> Roger
>>
>>
>>
>>> Here is the code, for the moment:
>>> ####
>>> library(rgeos)
>>>
>>> pls <- slot(bdc, "polygons")
>>> pls1 <- lapply(pls, checkPolygonsHoles)
>>> slot(bdc, "polygons") <- pls1
>>>
>>> routes.bdc = gIntersection(routes, bdc)
>>>
>>> routes.bdc = slot(routes.bdc, "lineobj")
>>>
>>> bdc.contains = gContains(bdc, routes.bdc, byid=TRUE)
>>>
>>> nbdc = apply(bdc.contains, 1, function(x) as.numeric(which(x)))
>>>
>>> apply(bdc.contains, 1, function(x) any(x))
>>> ###
>>>
>>> apply(bdc.contains, 1, function(x) any(x))
>>> mentions many FALSE values which indicate that some roads are not
>>> contained
>>> in any city...(?)
>>>
>>>
>>> Any help would be greatly appreciated!
>>>
>>> Thanks again
>>>
>>> 2011/6/7 Roger Bivand <Roger.Bivand at nhh.no>
>>>
>>>  On Tue, 7 Jun 2011, Mathieu Rajerison wrote:
>>>>
>>>>  Thanks for the answer!
>>>>
>>>>>
>>>>> I think I surely have to use gIntersection then gLength for this purpose
>>>>>
>>>>> I've tried gIntersection like this:
>>>>> routes2 = gIntersection(routes, bdc)
>>>>>
>>>>> But I get the following message:
>>>>> Erreur dans createPolygonsComment(p) :
>>>>>  rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon
>>>>> for
>>>>> hole at index 1
>>>>>
>>>>> Do you have an idea why I get this message and how to correct it?
>>>>>
>>>>>
>>>> The polygon object is malformed, and at least one polygon claims only to
>>>> be
>>>> an interior ring, which is impossible. You need to look at how you are
>>>> creating your SpatialPolygons object, and possibly use
>>>> checkPolygonsHoles()
>>>> in maptools to repair damaged objects. Typically running:
>>>>
>>>> library(maptools)
>>>> pls <- slot(bdc, "polygons")
>>>> pls1 <- lapply(pls, checkPolygonsHoles)
>>>> slot(bdc, "polygons") <- pls1
>>>>
>>>> may help. As Colin says, posting a link to the polygons object, and to
>>>> code
>>>> used to create it, may be helpful.
>>>>
>>>> Roger
>>>>
>>>>
>>>>
>>>>>
>>>>> 2011/6/7 Roman Lu??trik <roman.lustrik at gmail.com>
>>>>>
>>>>>  Package rgeos has all sorts of functions that might come handy.
>>>>>
>>>>>>
>>>>>> Cheers,
>>>>>> Roman
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Tue, Jun 7, 2011 at 6:13 PM, Mathieu Rajerison <
>>>>>> mathieu.rajerison at gmail.com> wrote:
>>>>>>
>>>>>>  Hi,
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> I have a SpatialLines object representing roads and a SpatialPolygons
>>>>>>> object
>>>>>>> containing cities.
>>>>>>>
>>>>>>> I'd like to know how to overlay a SpatialLines object by
>>>>>>> aSpatialPolygons
>>>>>>> object.
>>>>>>>
>>>>>>> I'd like to know the total distance of roads that cover each of my
>>>>>>> cities.
>>>>>>>
>>>>>>> Is it possible?
>>>>>>>
>>>>>>>
>>>>>>> Thanks,
>>>>>>>
>>>>>>> Mathieu
>>>>>>>
>>>>>>>      [[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
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> In God we trust, all others bring data.
>>>>>>
>>>>>>
>>>>>>        [[alternative HTML version deleted]]
>>>>>
>>>>>
>>>>>
>>>>>  --
>>>> 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
>>>>
>>>>
>>>>
>>>        [[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
>> 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