[R-sig-Geo] how to calculate areas of shape files using R

[Roger Bivand]

On Mon, 31 Jan 2011, Roman Luštrik wrote:

> When you convert your polygons to Spatial object, the area will be
> accessible through the slots. This can be very tedious to grasp at first
> (but try str()). For instance, for a SpatialPolygons of one Polygon

A small caveat - if the coordinates are planar, the explanation holds, if 
geographical (degrees), it doesn't. The naive area values have been used 
internally to plot larger Polygon objects before smaller Polygon objects, 
to avoid overplotting in earlier versions of R that did not "identify" 
holes automatically. So first make sure that you have planar coordinates, 
otherwise the "area" values have no metric. If need be, project with 
spTransform() methods in rgdal. If I recall, Elaine was working with an 
area covering much of the Western Pacific, hence the remark - hopefully, 
this question is about a project with a smaller extent.

In addition, if the Polygons objects are simple (one Polygon object per 
Polygons object), the "area" slot of the Polygons object is equal to that 
of the member Polygon, and:

sapply(slot(SP, "polygons"), slot, "area")

will return a vector of areas. Check with:

sapply(slot(SP, "polygons"), function(x) length(slot(x, "Polygons")))

If there are non-simple Polygons objects (with more than one Polygon 
member), we also need to check whether there are holes (or should be - try 
plot() too) - see below. If there are real holes not declared as holes, 
the total area will be inflated.

Example:

library(maptools)
fn <- system.file("shapes/sids.shp", package="maptools")[1]
prj <- CRS("+proj=longlat +ellps=clrk66")
xx <- readShapeSpatial(fn, IDvar="FIPSNO", proj4string=prj) 
library(rgdal)
xxx <- spTransform(xx, CRS("+init=epsg:3631"))
# NC State Plane latest
sapply(slot(xxx, "polygons"), function(x) length(slot(x, "Polygons")))
# some non-simple Polygons objects
sapply(slot(xxx, "polygons"), function(x) {
   xi <- slot(x, "Polygons")
   any(sapply(xi, slot, "hole"))
})
# no declared holes
sapply(slot(xxx, "polygons"), slot, "area")
# areas in square metres

Here the non-simple Polygons objects are counties made up of more than one 
part/island, with no holes. checkPolygonsHoles() in maptools may be used 
to check hole status if required.

Hope this helps,

Roger

>
>> my.object
> An object of class "SpatialPolygons"
> Slot "polygons":
> [[1]]
> An object of class "Polygons"
> Slot "Polygons":
> [[1]]
> An object of class "Polygon"
> Slot "labpt":
> [1] -38.80103  29.16010
>
> Slot "area":
> [1] 14371.66
>
> Slot "hole":
> [1] FALSE
>
> Slot "ringDir":
> [1] 1
>
> Slot "coords":
>           [,1]      [,2]
> [1,] -109.76555  61.91851
> [2,]   21.09385 105.52317
> [3,]   21.09385 -29.65128
> [4,] -105.85930 -25.29081
> [5,] -109.76555  61.91851
>
>
>
> Slot "plotOrder":
> [1] 1
>
> Slot "labpt":
> [1] -38.80103  29.16010
>
> Slot "ID":
> [1] "1"
>
> Slot "area":
> [1] 14371.66
>
>
>
> Slot "plotOrder":
> [1] 1
>
> Slot "bbox":
>         min       max
> x -109.76555  21.09385
> y  -29.65128 105.52317
>
> Slot "proj4string":
> CRS arguments: NA
>
> And for this particular case, you can access the area slot through
>
>
> my.object at polygons[[1]]@Polygons[[1]]@area
>
>
> Cheers,
> Roman
>
>
>
>
> On Mon, Jan 31, 2011 at 12:42 AM, elaine kuo <elaine.kuo.tw at gmail.com>wrote:
>
>> Dear List,
>>
>> I am dealing with hundreds of shape files. (generated by ArcGIS 9.3)
>> They are composed of multiple polygons.
>> The polygons are classified as 1, 2, 3, and 4 according to their seasonal
>> status.
>>
>> I want to calculate the total areas of each shape files.
>> Furthermore, the areas of 1, 2, 3, and 4 are also needed to be computed for
>> each shape file.
>>
>> Please kindly advise if it is possible to use R packages and functions are
>> suitable for the tasks.
>> (may be adehabitat)
>> Thank you.
>>
>> Elaine
>>
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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