[R-sig-Geo] rgeos::gArea() to all polygons within sp polygons DF
Roger Bivand
Roger.Bivand at nhh.no
Tue Apr 10 10:48:27 CEST 2012
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