calculate_pixel_area

## Original JavaScrpit code by Bill Hube and Rich's anwser in StackExchange
## https://gis.stackexchange.com/questions/127165/more-accurate-way-to-calculate-area-of-rasters
## May available in other CRSs by defining the right "a" and "b"
## Calculate m^2 area of a wgs84 square pixel.
## Adapted from: https://gis.stackexchange.com/a/127327/2397
## Parameters:
##      pixel_size (float): length of side of pixel in degrees.
##      center_lat (float): latitude of the center of the pixel. Note this
##      value +/- half the `pixel-size` must not exceed 90/-90 degrees
##      latitude or an invalid area will be calculated.
## Returns:
## Area of square pixel of side length `pixel_size` centered at
## `center_lat` in m^2

## "mosaic" package needed
## for users haven't download package "mosaic",
## please download the package by run following code in console :
## 'install.packages("mosaic")'

library(mosaic)
area_of_pixel <- function(center_lat){
  a = 6378137 # unit:meters
  b = 6356752.3142 # unit:meters
  c = sqrt(1 - (b/a)^2)
  
  zm_a <- 1 - c*sin(deg2rad(center_lat+1/2))
  zp_a <- 1 + c*sin(deg2rad(center_lat+1/2))
  area_a <- 
    pi * b^2 * (
      log(zp_a/zm_a) / (2*c) +
        sin(deg2rad(center_lat+1/2)) / (zp_a*zm_a))

  zm_b <- 1 - c*sin(deg2rad(center_lat-1/2))
  zp_b <- 1 + c*sin(deg2rad(center_lat-1/2))
  area_b <- 
  pi * b^2 * (
    log(zp_b/zm_b) / (2*c) +
      sin(deg2rad(center_lat-1/2)) / (zp_b*zm_b))

    return(1 / 360 * (area_a - area_b))
}