[R-sig-Geo] rgeos::gArea() to all polygons within sp polygons DF

Agustin Lobo alobolistas at gmail.com
Mon Apr 16 19:25:49 CEST 2012 

Roger,

In order to help testing, should we just install
sp, rgdal and rgeos from r-forge.r-project.org?
or is there an special way to force that we install the SVN revs that
you mention?

I will not be able to do it until early in Wednesday, as I'm away of
the office with just a little macbook air.

Agus

El día 16 de abril de 2012 16:52, Roger Bivand <Roger.Bivand at nhh.no> escribió:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

More information about the R-sig-Geo mailing list