How do I calculate the proportion of smaller squares covered by a larger circle? I would like to determine the proportion covered by the circle (say r = 2800m) of each cell (say 1000 x 1000 m). The circle is centered around the center of the focal cell of a rectangular neighborhood. I provide a rough visual estimation in the attached picture to clarify the idea. I.e. "1" means fully included, and "0" means no overlap.
The problem occurs in spatial modeling and coordinates are in meters (epsg:6933). This rules out sampling approaches as extremely inefficient (thanks for the suggestion though!). 
Ideally the solution should should be implemented easily in R.
Any suggestion would be highly welcome. Thanks.
[Yes , I know the dot is missing in two "0.09"s]

 A: This question asks how to construct a rasterized two-dimensional kernel for a circular uniform distribution.
The underlying mathematical problem is this: given a planar circle $R$ and a rectangle $X,$ find the area of $X\cap R$.  The kernel's values will be these areas divided by the area of $R$ (so that it integrates to unity).
This is an elementary integration problem but it's messy.  Borrowing statistical methods for describing probability distributions in two dimensions, we can construct the solution out of simpler ones.  These are:


*

*Shift and scale.  Choose the circle's center as the origin and choose units of measurement in which it has a unit radius.

*Create a "distribution function."  This function $F$ of two arguments $x$ and $y$ returns the area of the circle located to the left of $x$ and beneath $y$.  Using it we may express the area of the circle covered by the rectangle with lower left corner $(x_0,y_0)$ and upper right corner $(x_1,y_1)$ as
$$f(x_0,y_0; x_1,y_1) = F(x_1,y_1) - F(x_0,y_1) - F(x_1,y_0) + F(x_0,y_0).$$

*Exploit symmetries. By negating $x$ and or $y$ and swapping $x$ with $y$ as necessary we can relate $F(x,y)$ to the value of $F(u,v)$ where $u \ge v\ge 0$.

*The computation requires the "marginal" distributions $F_X(x)$, the area of the circle to the left of $x$, and $F_Y(y)$, the area beneath $y$. By symmetry $F_X=F_Y$.
By integration (or elementary trigonometry) we find that for $-1 \le x \le 1$,
$$F_X(x) = \pi/2  + \arcsin(x) + x\sqrt{1-x^2}.$$
(The pieces are twice the areas of a quarter circle, a sector, and a triangle, respectively.) This function rises from $0$ to $\pi$.
I leave the details of (1) - (3) to the interested reader (they can be found in the implementation of pucirc below).  They can all be worked out by carving $R$ along the vertical lines at $0$ and $x$ and the horizontal lines at $0$ and $y$.  There are really only two distinct cases to consider: $(x,y)$ lies inside the circle or outside it (that is, $x^2+y^2$ is less than $1$ or greater than $1$).  All the resulting pieces can be expressed in terms of $F_X$, simple constants, and simple functions (like $xy$, which is the area of the rectangle from the origin to the point $xy$).
To illustrate, here are the overlap areas in a unit grid for a central circle of radius $2.8$, as requested in the question.

The circle's boundary is overlaid in orange.  The areas of each cell are posted (rounded to $0.001$).  They were found by applying $f$ repeatedly to each of the cells: see the "Application" section in the code for details.
As a quick check, the last line of this code sums all the cell areas and compares them to the theoretical area of the circle.  The resulting discrepancy should be zero (and it is in this example).  Several tests (which graph various functions) have been left in but commented out.
#
# Area of a circle of radius `r`, centered at `origin`,
# to the left of `x` and below `y`.
#
pucirc <- function(x, y, radius=1, origin=c(0,0)) {
  #
  # Helper function.
  #
  xsqrt <- function(x) {
    y <- pmax(-1, pmin(1, x))
    y * sqrt(1 - y^2) + asin(y)
  }
  #
  # Area of the unit circle to the left of `x`.
  #
  pucirc_margin <- function(x) {
    ifelse(x >= 1, pi, ifelse(x <= -1, 0, pi/2 + xsqrt(x)))
  }
  # curve(pucirc_margin(x), from=-1.5, to=1.5)
  #
  # Area of the unit circle below `y` and to the left of `x` for
  # 0 <= x,y
  #
  pucirc_ <- function(x, y) {
    u <- pmax(x, y)
    v <- pmin(x, y)
    ifelse(u >= 1, pucirc_margin(v),
           ifelse(u^2+v^2 > 1, pucirc_margin(u) + pucirc_margin(v) - pi,
                  pucirc_margin(u)/2 + pucirc_margin(v)/2 - pi/4 + u*v))
  }
  # curve(pucirc_margin(x), from=0, to=1.5, ylim=c(0, pi), lwd=2)
  # curve(pucirc_(x, 0), add=TRUE, col="Gray", lty=3, lwd=2)
  # for (y in c(1/3, 0.9, 1)) {
  #   curve(pucirc_(x, y), from=0, to=1.5, add=TRUE, col=hsv(y * 2/3), lwd=2)
  # }
  #
  # Area of the unit circle below `y` and to the left of `x` for
  # any `x` and `y`.
  #
  pucirc0 <- function(x, y) {
    ifelse(x < 0 & y < 0, pi - pucirc_margin(-x) - pucirc_margin(-y) + pucirc_(-x, -y),
           ifelse(x < 0 & y >= 0, pucirc_margin(y) - pucirc_(-x, y),
                  ifelse(x >= 0 & y < 0, pucirc_margin(x) - pucirc_(x, -y),
                         pucirc_(x, y))))
  }
  # curve(pucirc_margin(x), from=-1.5, to=1.5, ylim=c(0, pi), lwd=2)
  # curve(pucirc0(x, 0), add=TRUE, col="Gray", lty=3, lwd=2)
  # for (y in c(-1, -1/2, 1/3, 0.9, 1)) {
  #   curve(pucirc0(x, y), add=TRUE, col=hsv((y+1) * 1/3), lwd=2)
  # }
  pucirc0((x-origin[1])/radius, (y-origin[2])/radius) * radius^2
}
#
# Area of a circle overlapping a rectangle with lower left corner `ll`
# and upper right corner `ur`.
#
pucirc_rect <- function(ll, ur, radius=1, origin=c(0,0)) {
  pucirc(ur[1], ur[2], radius, origin) - 
    pucirc(ur[1], ll[2], radius, origin) -
    pucirc(ll[1], ur[2], radius, origin) +
    pucirc(ll[1], ll[2], radius, origin)
}
#------------------------------------------------------------------------------#
#
# Application: construct a rasterized circular uniform kernel.
#
library(data.table) # Fast calculation
library(ggplot2)    # Graphical display

radius <- 2.8
rmax <- floor(2*radius) / 2     # Smallest rectangle coordinate
rmin <- floor(-2*radius-1) / 2  # Larges rectangle coordinate
#
# Create a table of all lower-left cell coordinates.
#
X <- as.data.table(expand.grid(i=seq(rmin, rmax, by=1), j=seq(rmin, rmax, by=1)))
#
# Compute the overlap areas for the cells.
#
f <- Vectorize(function(i,j) pucirc_rect(c(i,j), c(i,j)+1, radius=radius))
invisible(X[, Overlap := f(i,j)])
#
# Create data to plot the circle itself.
#
theta <- seq(0, 2*pi, length.out=361)
P <- data.table(x=cos(theta)*radius, y=sin(theta)*radius)
#
# Plot everything.
#
g <- ggplot(X, aes(i,j)) + 
  geom_raster(aes(fill=Overlap, alpha=Overlap), hjust=1, vjust=1) + 
  geom_text(aes(label=round(Overlap, 3)), nudge_x=1/2, nudge_y=1/2, size=3) + 
  geom_path(data=P, aes(x, y), size=1.5, alpha=2/3, color="Orange") + 
  guides(alpha="none", fill="none") + 
  coord_fixed(ratio=1)

print(g)
#
# Check.
#
sum(X$Overlap) / (pi * radius^2) - 1 # Should be zero.

A: My suggestion is to use sampling. Create a matrix of your squares. Take samples from two uniform variables, interpret those as x and y values. Keep them only if they fall in the circle. Count those towards the square in which they belong. The estimated fraction of a square in a circle is then number of samples it caught, divided by the expected number of samples (= total samples / number of squares). Not the most efficient solution, but, very easy to implement, and it will work if you take enough samples. 
A: Below is R code. I hope it is reproducible for everyone.
Cell size is set in 4th block of code with "res"
The radius is set in the 2ndlast block/line with "width".
Missing: generalization to omit providing a reasonable estimate
require(raster)
require(rgdal) 
require(rgeos)

epsg6933 <- CRS("+init=epsg:6933")

study_extent <- extent(5836086, 14150225, -1388257, 6628905)

rasterdummy <- raster(x = study_extent, crs = epsg6933, res = c(1000, 1000))

# next comes the estimate as in the picture. could be 0s and 1s though.
focweight <- matrix(c(0,    0,    0.09, 0.25, 0.09, 0,    0, 
                      0,    0.45, 0.99, 1,    0.99, 0.45, 0,
                      0.09, 0.99, 1,    1,    1,    0.99, 0.09,
                      0.25, 1,    1,    1,    1,    1,    0.25,
                      0.09, 0.99, 1,    1,    1,    0.99, 0.09,
                      0,    0.45, 0.99, 1,    0.99, 0.45, 0,
                      0,    0,    0.09, 0.25, 0.09, 0,    0), nrow = 7)

focweight[focweight!=0] <- cellFromPolygon(object = rasterdummy, p = gBuffer(spgeom =  SpatialPoints(coords = matrix(c(study_extent@xmin + 10500, study_extent@ymax - 10500), ncol = 2), proj4string = epsg6933, bbox = bbox(rasterdummy)), quadsegs = 720, width = 2800), weights = T)[[1]][,2]

focweight # to print values to stdout

Result:

