Calculating areas on Earth’s surface

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Calculating areas on Earth’s surface

Karl Ove Hufthammer
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@polygons, function(x) x@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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: 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@polygons, function(x) x@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
>