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

Roger Bivand Roger.Bivand at nhh.no
Mon Apr 16 20:45:32 CEST 2012 

On Mon, 16 Apr 2012, Agustin Lobo wrote:

> 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?

Yes, but from source, I'm afraid, only sp is in current build there, while 
rgdal and rgeos are not built.

Roger

>
> 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
>>
>

-- 
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