[R-sig-Geo] Calculating areas on Earth’s surface

[Robert Hijmans]

Dear Karl,

For the area of a spherical polygon, you could do

library(geosphere)
areaPolygon(xy)

I do not know how accurate this is relative to your alternative approaches. 

Robert


> Calculating areas on Earth’s surface 

> I’m interested in the distortion in apparent area from various
> projections, 
> and need to calculate the *real* area of various polygons (or grid cells). 
> In an earlier post on this list, it was suggested to use a equal-area 
> projection. This partly works, but the accuracy is not good enough for 
> measuring the distortion for my purpose (I guess an accuracy < 0.05% for 
> moderately large areas would be OK). 
>  
> Here’s a simple example: 
>
> library(sp) 
> library(rgdal) 
>
> unproj=CRS("+init=epsg:4326")                  # Unprojected, i.e.
> longitudes and latitudes 
> xy=GridTopology(c(0,55),c(8, 4), c(5,5))       # Small grid (covering
> Norway) 
> xy=as.SpatialPolygons.GridTopology(xy, unproj) # Grid as polygons 
>
> # Calculate the area of a polygons using the 'proj' projection 
> calcArea=function(xy, proj) { 
>  xy.trans=spTransform(xy, proj) 
>  sapply( xy.trans at polygons, function(x) x at area)/1e6 
> } 
>
> # Various equal-area projections 
> proj1=CRS("+proj=moll +lat_0=65 +lon_0=10") # Mollweide 
> proj2=CRS("+proj=sinu +lat_0=65 +lon_0=10") # Sinusoidal 
> proj3=CRS("+proj=tcea +lat_0=65 +lon_0=10") # Transverse Cylindrical 
> proj4=CRS("+init=epsg:32633")               # UTM 33N (*not* equal area) 
>
> # Areas according to the different projections all differ 
> a1=calcArea(xy, proj1) 
> a2=calcArea(xy, proj2) 
> a3=calcArea(xy, proj3) 
> a4=calcArea(xy, proj4) 
>
> # All areas calculated using Mollweide are small than the areas 
> # calculated using the sinusoidal projection: 
> all(a1<a2) 
>
> Since the areas from the various projections all differ, I’m not sure
> which 
> one is most accurate. Does anyone know of a ‘real’ area function available 
> somewhere in R, which can calculate the area of simple polygons on the WGS
> 84 
> ellipsoid (or perhaps even a spherical approximation would be good
> enough). 
> It would be very useful for drawing maps showing the distortion in area
> for 
> various projections. 
>
> -- 
> Karl Ove Hufthammer 
>


--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Calculating-areas-on-Earth-s-surface-tp6606454p6704251.html
Sent from the R-sig-geo mailing list archive at Nabble.com.