[R-sig-Geo] how to calculate areas of shape files using R
Roger Bivand
Roger.Bivand at nhh.no
Mon Jan 31 09:59:20 CET 2011
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
More information about the R-sig-Geo
mailing list