# R code from www.graemetlloyd.com
#
# The following implements the function described by Bevis and Cambareri, 1987
# (Mathematical Geology, 19(4), 335-346). It additionally assumes the 'inside'
# area of the polygon is the target of interest.
# 
# Vectors of latitudes and longitudes should be fed to the functions 'vlat'
# and 'vlon' respectively and the area will be returned.
#
# Note: N latitudes and E longitudes should be entered as positives, S latitudes 
# and W longitudes as negatives.

sphpolyarea<-function(vlat, vlon)
{ 
  points<-chull(vlon,vlat)
  vlon<-vlon[points]
  vlat<-vlat[points]
  sum<-0
  nv<-length(vlat)
  dtr<-pi/180
  rad<-6371 # Radius of Earth given in kilometres
  for (i in 1:nv) {
    if (i == 1) {
      flat<-vlat[2]
      flon<-vlon[2]
      blat<-vlat[nv]
      blon<-vlon[nv]
    }
    if (i != 1 && i < nv) {
      flat<-vlat[i+1]
      flon<-vlon[i+1]
      blat<-vlat[i-1]
      blon<-vlon[i-1]
    }
    if (i == nv) {
      flat<-vlat[1]
      flon<-vlon[1]
      blat<-vlat[nv-1]
      blon<-vlon[nv-1]
    }
    plat<-vlat[i]
    plon<-vlon[i]
    qlat<-flat
    qlon<-flon
    t<-sin((qlon-plon)*dtr)*cos(qlat*dtr)
    b<-sin(dtr*qlat)*cos(plat*dtr)-cos(qlat*dtr)*sin(plat*dtr)*cos((qlon-plon)*dtr)
    fang<-atan2(t,b)
    plat<-vlat[i]
    plon<-vlon[i]
    qlat<-blat
    qlon<-blon
    t<-sin((qlon-plon)*dtr)*cos(qlat*dtr)
    b<-sin(dtr*qlat)*cos(plat*dtr)-cos(qlat*dtr)*sin(plat*dtr)*cos((qlon-plon)*dtr)
    bang<-atan2(t,b)
    fvb<-bang-fang
    if (fvb < 0) fvb<-fvb+(2*pi)
    sum<-sum+fvb
  }
  area<-(sum-pi*(nv-2))*rad^2
  area
}