rgeos::gArea() to all polygons within sp polygons DF

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

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

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

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

print(b2[25967], digits=16)

print(b2[26570], digits=16)

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

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

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

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

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

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

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

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

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.

system.time(bn <- readOGR(".", "BorneoLCv2google"))

system.time(b2 <- gArea(bn, byid=TRUE))

system.time(bn <- createSPComment(bn))

set_do_poly_check(FALSE)

system.time(b2 <- gArea(bn, byid=TRUE))

print(area[25967], digits=16)

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

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

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

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

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

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

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

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

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

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

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

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

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

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