[R-sig-Geo] rgeos::gArea() to all polygons within sp polygons DF
Roger Bivand
Roger.Bivand at nhh.no
Mon Apr 16 16:52:40 CEST 2012
This is a request to try out development versions of sp, rgdal, and rgeos,
related to Agus' question.
I've commited sp SVN rev. 1246 (in rspatial), rgdal SVN rev. 325, and
rgeos SVN rev. 332 to R-Forge.
The objects defined in sp for polygon geometries do not follow OGC SFS
specifications, but can be mapped to them by pointers assigning interior
rings to their containing exterior rings, stored in comments assigned to
Polygons objects. These comments may be created in maptools or rgeos. It
also turns out that the required assignments can be extracted from input
OGR data.
By adding a top-level comment to SpatialPolygons* objects set to TRUE if
all member Polygons objects have comments set, we can avoid the very
costly comment assignment loop. If in addition the SpatialPolygons*
objects are read with readOGR() in rgdal, the comments can be set directly
while reading, as OGR seems to expect the first ring in an OGC SFS Polygon
object to be exterior, the rest interior.
This means that for Agus' example, there is no extra cost involved in
setting the comments, and his gArea call is 0.25s rather than over 100s
with comment setting in rgeos, and more with full containment analysis via
maptools (shapelib in maptools does not support OGC SFS geometries).
So I would like to receive feedback on reading polygon data with readOGR()
and on whether it seems reasonable to accept the assignment of interior
rings and exterior rings for Polygons objects. If we can believe the
Polygons coming from OGR, we can save a lot of time in preparing them for
further handling in rgeos.
I do not want to release these revisions until they have been tried out
thoroughly.
Roger
On Tue, 10 Apr 2012, Roger Bivand wrote:
> On Mon, 9 Apr 2012, Agustin Lobo wrote:
>
>> Thanks Roger.
>> Apologies for the missing output of sessionInfo() .
>> Confirmed that the results are correct with rgeos 0.2-5, although it's
>> now very slow.
>
> You can say:
>
> bn <- createSPComment(bn)
>
> first, which takes some time as it recreates the object. If the comments
> assigning holes (interior rings) to exterior rings are present, they are not
> recreated. Shapefiles do not store this information, and it is flattened out
> in making Polygons objects, which were (perhaps unfortunately) conceived of
> as matching shapefiles.
>
> Arguably, this comment creation, needed to map sp Polygons objects to OGC SFS
> Polygon or Multipolygon objects, should be done in sp - but the computational
> geometry functions neede to find our which hole is in which exterior ring are
> in GEOS, so this would involve a circularity in dependencies.
>
> So the way to ensure that polygon objects work as expected in rgeos and when
> we believe the input hole status is:
>
> library(rgdal)
> pls <- readOGR(...)
> library(rgeos)
> pls <- createSPComment(pls)
> # I'll look at making this faster in C
> # subsequent operations (which still check that comments are present in #
> each Polygons object)
>
> For your case, which has many very complex polygons and as such is a good
> test case, I have:
>
>> system.time(bn <- readOGR(".", "BorneoLCv2google"))
> OGR data source with driver: ESRI Shapefile
> Source: ".", layer: "BorneoLCv2google"
> with 32036 features and 6 fields
> Feature type: wkbPolygon with 2 dimensions
> user system elapsed
> 9.928 0.178 10.388
>> system.time(b2 <- gArea(bn, byid=TRUE))
> user system elapsed
> 154.242 10.033 169.401
>> system.time(bn <- createSPComment(bn))
> user system elapsed
> 152.045 0.085 156.501
>> set_do_poly_check(FALSE)
> # without this, it redoes the comment checking
>> system.time(b2 <- gArea(bn, byid=TRUE))
> user system elapsed
> 0.251 0.000 0.252
>
> In this case, it is also possible to get the areas without using rgeos:
>
> pls <- slot(bn, "polygons")
> exter <- sapply(pls, function(x) {
> Pls <- slot(x, "Polygons")
> sum(sapply(Pls, function(xx)
> ifelse(!slot(xx, "hole"), slot(xx, "area"), 0)))
> })
> inter <- sapply(pls, function(x) {
> Pls <- slot(x, "Polygons")
> sum(sapply(Pls, function(xx)
> ifelse(slot(xx, "hole"), slot(xx, "area"), 0)))
> })
> area <- exter - inter
>
>> print(area[25967], digits=16)
> [1] 8550196769.930017
>
> since the net area of each Polygons object is the gross area minus the holes;
> this is also fast.
>
> Hope this clarifies,
>
> Roger
>
>> Ticket in QGIS updated (and closed).
>>
>> Agus
>>
>> El día 9 de abril de 2012 17:36, Roger Bivand <Roger.Bivand at nhh.no>
>> escribió:
>>> On Mon, 9 Apr 2012, Agustin Lobo wrote:
>>>
>>>> Thanks Edzer,
>>>>
>>>> b2 <- gArea(BorneoLCg, byid=T)
>>>>
>>>> I get many warnings such as:
>>>> 40: In RGEOSMiscFunc(spgeom, byid, "rgeos_area") :
>>>> Polygons object missing comment attribute ignoring hole(s). See
>>>> function createSPComment.
>>>>
>>>> is this normal?
>>>>
>>>> The result is identical to
>>>> b <- sapply(BorneoLCg at polygons, function(x) x at area)
>>>>
>>>> while it is different than the one provided by QGIS (known to be wrong
>>>> for
>>>> very
>>>> large polygons only because of a format problem) and postgis (assumed
>>>> to be correct). Please give a look to:
>>>> http://hub.qgis.org/issues/5332
>>>
>>>
>>> Please always provide sessionInfo() output, and a copy of the rgeos
>>> startup
>>> messages (showing which version of GEOS is used). With rgeos 0.2-5 and
>>> GEOS
>>> runtime version: 3.3.2-CAPI-1.7.2, I have:
>>>
>>> bn <- readOGR(".", "BorneoLCv2google")
>>> b2 <- gArea(bn, byid=TRUE)
>>>>
>>>> print(b2[25967], digits=16)
>>>
>>> 25966
>>> 8550196769.93012
>>>>
>>>> print(b2[26570], digits=16)
>>>
>>> 26569
>>> 173647473841.3217
>>>
>>> so the hole warnings were the reason for the difference. SVN revision 327
>>> of
>>> 16 March enforced the use of createSPComment() in gArea() - the user was
>>> expected to denote holes properly themselves before this change. Update
>>> your
>>> rgeos and try again.
>>>
>>> Please update the QGIS #5332 to state that from rgeos 0.2-5, the values
>>> are
>>> correct. The area slot in the Polygons object is a gross area used to
>>> order
>>> plotting (because until recently R polygons simply overplotted one
>>> another),
>>> so it is not guaranteed to be a planar area measure for polygons with
>>> holes.
>>>
>>> Roger
>>>
>>>
>>>>
>>>> Agus
>>>>
>>>>
>>>> El día 8 de abril de 2012 17:26, Edzer Pebesma
>>>> <edzer.pebesma at uni-muenster.de> escribió:
>>>>>
>>>>> Agus, have you tried
>>>>>
>>>>> gArea(BorneoLCg, byid = TRUE)
>>>>>
>>>>> ?
>>>>>
>>>>> On 04/08/2012 04:52 PM, Agustin Lobo wrote:
>>>>>>
>>>>>> Hi!
>>>>>>
>>>>>> How do I apply gArea()) to all polygons within an spatial pol. DF?
>>>>>> I've tried
>>>>>> areas <- apply(BorneoLCg,gArea)
>>>>>> and
>>>>>> areas <- sapply(BorneoLCg,gArea)
>>>>>>
>>>>>> but this is not correct.
>>>>>>
>>>>>> I want to compare to
>>>>>> areas <- sapply(BorneoLCg at polygons, function(x) x at area)
>>>>>>
>>>>>> Thanks
>>>>>>
>>>>>> Agus
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> R-sig-Geo at r-project.org
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>>
>>>>> --
>>>>> Edzer Pebesma
>>>>> Institute for Geoinformatics (ifgi), University of Münster
>>>>> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
>>>>> 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
>>>>> http://www.52north.org/geostatistics e.pebesma at wwu.de
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo at r-project.org
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>> _______________________________________________
>>>> 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
>>>
>>> _______________________________________________
>>> 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
More information about the R-sig-Geo
mailing list