area Subroutine

public subroutine area(a, f, lats, lons, n, AA, PP)

Determine the area of a geodesic polygon

Arbitrarily complex polygons are allowed. In the case of self-intersecting polygons the area is accumulated "algebraically", e.g., the areas of the 2 loops in a figure-8 polygon will partially cancel. There's no need to "close" the polygon by repeating the first vertex. The area returned is signed with counter-clockwise traversal being treated as positive.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: a

the equatorial radius (meters).
real(kind=wp), intent(in) :: f

the flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
real(kind=wp), intent(in) :: lats(n)

an array of the latitudes of the vertices (degrees). lats should be in the range [-90 deg, 90 deg].
real(kind=wp), intent(in) :: lons(n)

an array of the longitudes of the vertices (degrees).
integer, intent(in) :: n

number of points
real(kind=wp), intent(out) :: AA

the (signed) area of the polygon ().
real(kind=wp), intent(out) :: PP

the perimeter of the polygon.

Source Code

subroutine area(a, f, lats, lons, n, AA, PP)

  integer, intent(in) :: n !! number of points
  real(wp), intent(in) :: a !! the equatorial radius (meters).
  real(wp), intent(in) :: f !! the flattening of the ellipsoid.  Setting `f = 0` gives
                            !! a sphere.  Negative `f` gives a prolate ellipsoid.
  real(wp), intent(in) :: lats(n) !! an array of the latitudes of the vertices (degrees).
                                  !! lats should be in the range [-90 deg, 90 deg].
  real(wp), intent(in) :: lons(n) !! an array of the longitudes of the vertices (degrees).
  real(wp), intent(out) :: AA !! the (signed) area of the polygon (\(m^2\)).
  real(wp), intent(out) :: PP !! the perimeter of the polygon.

integer :: i, outmask, cross
real(wp) :: s12, azi1, azi2, dummy1,dummy2,dummy3,dummy4, SS12, b, e2, c2, area0, &
            Aacc(2), Pacc(2)

outmask = 8
Aacc = 0.0_wp
Pacc = 0.0_wp
cross = 0
do i = 0, n-1
  call inverse(a, f, lats(i+1), lons(i+1), &
      lats(mod(i+1,n)+1), lons(mod(i+1,n)+1), &
      s12, azi1, azi2, outmask, dummy1, dummy2, dummy3, dummy4, SS12)
  call accadd(Pacc, s12)
  call accadd(Aacc, -SS12)
  cross = cross + transit(lons(i+1), lons(mod(i+1,n)+1))
end do
PP = Pacc(1)
b = a * (1.0_wp - f)
e2 = f * (2.0_wp - f)
if (e2 == 0.0_wp) then
  c2 = a**2
else if (e2 > 0.0_wp) then
  c2 = (a**2 + b**2 * atanh(sqrt(e2)) / sqrt(e2)) / 2.0_wp
else
  c2 = (a**2 + b**2 * atan(sqrt(abs(e2))) / sqrt(abs(e2))) / 2.0_wp
end if
area0 = 4.0_wp * pi * c2
call accrem(Aacc, area0)
if (mod(abs(cross), 2) == 1) then
  if (Aacc(1) < 0.0_wp) then
    call accadd(Aacc, +area0/2.0_wp)
  else
    call accadd(Aacc, -area0/2.0_wp)
  end if
end if
if (Aacc(1) > area0/2.0_wp) then
  call accadd(Aacc, -area0)
else if (Aacc(1) <= -area0/2.0_wp) then
  call accadd(Aacc, +area0)
end if
AA = Aacc(1)

end subroutine area