[R-sig-Geo] Computing polygon area with decimal degree coordinates

[Roger Bivand]

On Wed, 31 May 2006, Patrick Giraudoux wrote:

> Well, thinking about it, I realize this question was a bit weird/silly 
> in the real world... Decimal degree coordinates does not tell us 
> anything about the ellipsoid/datum they refer to and this must be done 
> within a CRS framework. If we want a reasonable precision on small 
> areas, sound mathematics to get reasonable projections in an euclidean 
> space are likely unavoidable (which is easy within rgdal and 
> sptranform()). Except if one think that the earth is a perfect sphere. 
> Would even be easier if still flat...

The Great Circle distance in sp assumes WGS84, and is probably OK for 
provinces, and maybe larger counties where the area should just be an 
order of magnitude, and should preserve the rank order of the polygons. 
But it cannot be accurate anyway, because it is limited to the straight 
line segments of the points on the polygon boundary - if they are thinned, 
you will get a different area, the same if they increase in detail, such 
as along a river bank. 

But it could be done (roughly), if the spTransform() route seems
unhelpful, or if the data are at a continental scale, where projection
will lead to big errors away from the map centre.

Roger

> 
> 
> Roger Bivand a écrit :
> > On Wed, 31 May 2006, Patrick Giraudoux wrote:
> >
> >   
> >> Dear Listers,
> >>
> >> The function areapl() of the package splancs computes polygon areas  in 
> >> the coordinate units. This means in square meters when using UTM or 
> >> Lambert projections but in "square degrees" when using longlat degrees.
> >>
> >> Does anybody knows a R function (or an algorithm) to compute the area 
> >> directly (without projecting or projecting internally), eg in square 
> >> meters, of a polygon whose node coordinates are given in decimal degrees?
> >>     
> >
> > I think it could be done internally in the sp/src/Rcentroid.c function if 
> > Area2() was rewritten to call sp_gcdist() on each of the four segments in:
> >
> >   area = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]);
> >
> > but it would only be acceptable for larger areas. As it is, the area of 
> > SpatialPolygons objects is really only used to create the plot order (to 
> > plot larger polygons before smaller ones). In fact we only have the 
> > measurement at the Polygon object level sensibly - above that we don't 
> > really know whether multiple polygons in a Polygons object are holes or 
> > not, so summing the component Polygon object areas in a Polygons object 
> > may be wrong.
> >
> > Roger
> >
> >   
> >> Thanks for any hint,
> >>
> >> Patrick
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> 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