### Function to calculate regional area-weighted average from a matrix X (lon x lat)
### from regular cartesian grid
### Inputs: data matrix (lon x lat), lat centers vector, lon centers vector

calc_awa = function(xx, lon, lat){ # lon in deg_east, lat in deg_north (S-N)
  R = 6371 # Earth's radius (km)
  dy = abs(lat[2]-lat[1]); dx = abs(lon[2] - lon[1]);
  lon.edges = lon-(dx/2); lon.edges[length(lon.edges)+1] = max(lon)+(dx/2)
  lat.edges = lat-(dy/2); lat.edges[length(lat.edges)+1] = max(lat)+(dy/2)
  lat.radians = lat.edges*(pi/180) # Need to convert to radians for sin/cos
  total.area = (pi/180)*(R^2)*abs(sin(max(lat.radians))-sin(min(lat.radians)))*abs(max(lon.edges)-min(lon.edges))
  grid.area = array(NA, dim=dim(xx))
  for (i in 1:length(lon)){
    for (j in 1:length(lat)){
      grid.area[i,j] = (pi/180)*(R^2)*abs(sin(lat.radians[j+1])-sin(lat.radians[j]))*abs(lon.edges[i+1]-lon.edges[i])
    }
  }
  res = sum(xx*grid.area, na.rm=TRUE)/total.area
  return(res) 
}

#===================================================================================================================