How to find the covariance matrix of a polygon? Imagine you have a polygon defined by a set of coordinates $(x_1,y_1)...(x_n,y_n)$ and its centre of mass is at $(0,0)$. You can treat the polygon as a uniform distribution with a polygonal boundary.

I'm after a method that will find the covariance matrix of a polygon. 
I suspect that the covariance matrix of a polygon is closely related to the second moment of area, but whether they are equivalent I'm not sure. The formulas found in the wikipedia article I linked seem (a guess here, it's not especially clear to me from the article) to refer to the rotational inertia around the x, y and z axes rather than the principal axes of the polygon.
(Incidentally, if anyone can point me to how to calculate the principal axes of a polygon, that would also be useful to me)
It is tempting to just perform PCA on the coordinates, but doing so runs into the issue that the coordinates are not necessarily evenly spread around the polygon, and are therefore not representative of the density of the polygon. An extreme example is the outline of North Dakota, whose polygon is defined by a large number of points following the Red river, plus only two more points defining the western edge of the state. 
 A: Edit: Didn't notice that whuber had already answered. I'll leave this up as an example of another (perhaps less elegant) approach to the problem.
The covariance matrix
Let $(X,Y)$ be a random point from the uniform distribution on a polygon $P$ with area $A$. The covariance matrix is:
$$C = \begin{bmatrix} C_{XX} & C_{XY} \\ C_{XY} & C_{YY} \end{bmatrix}$$
where $C_{XX} = E[X^2]$ is the variance of $X$, $C_{YY} = E[Y^2]$ is the variance of $Y$, and $C_{XY} = E[XY]$ is the covariance between $X$ and $Y$. This assumes zero mean, since the polygon's center of mass is located at the origin. The uniform distribution assigns constant probability density $\frac{1}{A}$ to every point in $P$, so:
$$C_{XX} = \frac{1}{A} \underset{P}{\iint} x^2 dV \quad
C_{YY} = \frac{1}{A} \underset{P}{\iint} y^2 dV \quad
C_{XY} = \frac{1}{A} \underset{P}{\iint} x y dV
\tag{1}$$
Triangulation
Instead of trying to directly integrate over a complicated region like $P$, we can simplify the problem by partitioning $P$ into $n$ triangular subregions:
$$P = T_1 \cup \cdots \cup T_n$$
In your example, one possible partitioning looks like this:

There are various ways to produce a triangulation (see here). For example, you could compute the Delaunay triangulation of the vertices, then discard edges that fall outside $P$ (since it may be nonconvex as in the example).
Integrals over $P$ can then be split into sums of integrals over the triangles:
$$C_{XX} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} x^2 dV \quad
C_{YY} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} y^2 dV \quad
C_{XY} = \frac{1}{A} \sum_{i=1}^n \underset{T_i}{\iint} x y dV
\tag{2}$$
A triangle has nice, simple boundaries so these integrals are easier to evaluate.
Integrating over triangles
There are various ways to integrate over triangles. In this case, I used a trick that involves mapping a triangle to the unit square. Transforming to barycentric coordintes might be a better option.
Here are solutions for the integrals above, for an arbitrary triangle $T$ defined by vertices $(x_1,y_1), (x_2,y_2), (x_3,y_3)$. Let:
$$v_x = \left[ \begin{smallmatrix} x_1 \\ x_2 \\ x_3 \end{smallmatrix} \right]
\quad
v_y = \left[ \begin{smallmatrix} y_1 \\ y_2 \\ y_3 \end{smallmatrix} \right]
\quad
\vec{1} = \left[ \begin{smallmatrix} 1 \\ 1 \\ 1 \end{smallmatrix} \right]
\quad
L = \left[ \begin{smallmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{smallmatrix} \right]$$
Then:
$$\underset{T}{\iint} x^2 dV =
\frac{A}{6} \text{Tr}(v_x v_x^T L) \quad
\underset{T}{\iint} y^2 dV =
\frac{A}{6} \text{Tr}(v_y v_y^T L) \quad
\underset{T}{\iint} x y dV =
\frac{A}{12} (\vec{1}^T v_x v_y^T \vec{1} + v_x^T v_y)
\tag{3}$$
Putting everything together
Let $v_x^i$ and $v_y^i$ contain the x/y coordinates of the vertices for each triangle $T_i$, as above. Plug $(3)$ into $(2)$ for each triangle, noting that the area terms cancel out. This gives the solution:
$$C_{XX} =
\frac{1}{6} \sum_{i=1}^n \text{Tr} \big( v_x^i (v_x^i)^T L \big) \quad
C_{YY} =
\frac{1}{6} \sum_{i=1}^n \text{Tr} \big( v_y^i (v_y^i)^T L \big) \quad
C_{XY} =
\frac{1}{12} \sum_{i=1}^n \big( \vec{1}^T v_x^i (v_y^i)^T \vec{1} + (v_x^i)^T v_y^i \big)
\tag{4}$$
Principal axes
The principal axes are given by the eigenvectors of the covariance matrix $C$, just as in PCA. Unlike PCA, we have an analytic expression for $C$, rather than having to estimate it from sampled data points. Note that the vertices themselves are not a representative sample from the uniform distribution on $P$, so one can't simply take the sample covariance matrix of the vertices. But, $C$ *is* a relatively simple function of the vertices, as seen in $(4)$.
A: Let's do some analysis first.
Suppose within the polygon $\mathcal{P}$ its probability density is proportional function $p(x,y).$  Then the constant of proportionality is the inverse of the integral of $p$ over the polygon,
$$\mu_{0,0}(\mathcal{P})=\iint_{\mathcal P} p(x,y) \mathrm{d}x\,\mathrm{d}y.$$
The barycenter of the polygon is the point of average coordinates, computed as their first moments.  The first one is
$$\mu_{1,0}(\mathcal{P})=\frac{1}{\mu_{0,0}(\mathcal{P})} \iint_{\mathcal P} x\,p(x,y)\mathrm{d}x\,\mathrm{d}y.$$
The inertial tensor can be represented as the symmetric array of second moments computed after translating the polygon to put its barycenter at the origin: that is, the matrix of central second moments
$$\mu^\prime_{k,l}(\mathcal{P}) = \frac{1}{\mu_{0,0}(\mathcal{P})} \iint_{\mathcal P} \left(x - \mu_{1,0}(\mathcal{P})\right)^k\,\left(y - \mu_{0,1}(\mathcal{P})\right)^l\,p(x,y)\mathrm{d}x\,\mathrm{d}y$$
where $(k,l)$ range from $(2,0)$ to $(1,1)$ to $(0,2).$  The tensor itself--aka covariance matrix--is
$$I(\mathcal{P}) = \pmatrix{\mu^\prime_{2,0}(\mathcal{P})  & \mu^\prime_{1,1}(\mathcal{P}) \\ \mu^\prime_{1,1}(\mathcal{P})  &  \mu^\prime_{0,2}(\mathcal{P})}.$$
A PCA of $I(\mathcal{P})$ yields the principal axes of $\mathcal{P}:$ these are the unit eigenvectors scaled by their eigenvalues.

Next, let's work out how to do the calculations.  Because the polygon is presented as a sequence of vertices describing its oriented boundary $\partial\mathcal P,$ it is natural to invoke

Green's Theorem: $$\iint_{\mathcal{P}} \mathrm{d}\omega = \oint_{\partial\mathcal{P}}\omega$$ where $\omega = M(x,y)\mathrm{d}x + N(x,y)\mathrm{d}y$ is a one-form defined in a neighborhood of $\mathcal{P}$ and $$\mathrm{d}\omega = \left(\frac{\partial}{\partial x}N(x,y) - \frac{\partial}{\partial y}M(x,y)\right)\mathrm{d}x\,\mathrm{d}y.$$

For instance, with $\mathrm{d}\omega = x^k y^l \mathrm{d}x\mathrm{d}y$ and constant (i.e., uniform) density $p,$ we may (by inspection) select one of the many solutions, such as $$\omega(x,y) = \frac{-1}{l+1}x^k y^{l+1}\mathrm{d}x.$$
The point of this is that the contour integral follows the line segments determined by the sequence of vertices.  Any line segment from vertex $\mathbf{u}$ to vertex $\mathbf{v}$ can be parameterized by a real variable $t$ in the form
$$t \to \mathbf{u} + t\mathbf{w}$$
where $\mathbf{w} \propto \mathbf{v}-\mathbf{u}$ is the unit normal direction from $\mathbf{u}$ to $\mathbf{v}.$ The values of $t$ therefore range from $0$ to $|\mathbf{v}-\mathbf{u}|.$  Under this parameterization $x$ and $y$ are linear functions of $t$ and $\mathrm{d}x$ and $\mathrm{d}y$ are linear functions of $\mathrm{d}t.$  Thus the integrand of the contour integral over each edge becomes a polynomial function of $t,$ which is easily evaluated for small $k$ and $l.$

Implementing this analysis is as straightforward as coding its components.  At the lowest level we will need a function to integrate a polynomial one-form over a line segment.  Higher level functions will aggregate these to compute the raw and central moments to obtain the barycenter and inertial tensor, and finally we can operate on that tensor to find the principal axes (which are its scaled eigenvectors).  The R code below performs this work.  It makes no pretensions of efficiency: it is intended only to illustrate the practical application of the foregoing analysis.  Each function is straightforward and the naming conventions parallel those of the analysis.
Included in the code is a procedure to generate valid closed, simply connected, non-self-intersecting polygons (by randomly deforming points along a circle and including the starting vertex as its final point in order to create a closed loop).  Following this are a few statements to plot the polygon, display its vertices, adjoin the barycenter, and plot the principal axes in red (largest) and blue (smallest), creating a polygon-centric positively-oriented coordinate system.

#
# Integrate a monomial one-form x^k*y^l*dx along the line segment given as an 
# origin, unit direction vector, and distance.
#
lintegrate <- function(k, l, origin, normal, distance) {
  # Binomial theorem expansion of (u + tw)^k
  expand <- function(k, u, w) {
    i <- seq_len(k+1)-1
    u^i * w^rev(i) * choose(k,i)
  }
  # Construction of the product of two polynomials times a constant.
  omega <- normal[1] * convolve(rev(expand(k, origin[1], normal[1])), 
                                expand(l, origin[2], normal[2]),
                                type="open")
  # Integrate the resulting polynomial from 0 to `distance`.
  sum(omega * distance^seq_along(omega) / seq_along(omega))
}
#
# Integrate monomials along a piecewise linear path given as a sequence of
# (x,y) vertices.
#
cintegrate <- function(xy, k, l) {
  n <- dim(xy)[1]-1 # Number of edges
  sum(sapply(1:n, function(i) {
    dv <- xy[i+1,] - xy[i,]               # The direction vector
    lambda <- sum(dv * dv)
    if (isTRUE(all.equal(lambda, 0.0))) {
      0.0
    } else {
      lambda <- sqrt(lambda)              # Length of the direction vector
      -lintegrate(k, l+1, xy[i,], dv/lambda, lambda) / (l+1)
    }
  }))
}
#
# Compute moments of inertia.
#
inertia <- function(xy) {
  mass <- cintegrate(xy, 0, 0)
  barycenter = c(cintegrate(xy, 1, 0), cintegrate(xy, 0, 1)) / mass
  uv <- t(t(xy) - barycenter)   # Recenter the polygon to obtain central moments
  i <- matrix(0.0, 2, 2)
  i[1,1] <- cintegrate(uv, 2, 0)
  i[1,2] <- i[2,1] <- cintegrate(uv, 1, 1)
  i[2,2] <- cintegrate(uv, 0, 2)
  list(Mass=mass,
       Barycenter=barycenter,
       Inertia=i / mass)
}
#
# Find principal axes of an inertial tensor.
#
principal.axes <- function(i.xy) {
  obj <- eigen(i.xy)
  t(t(obj$vectors) * obj$values)
}
#
# Construct a polygon.
#
circle <- t(sapply(seq(0, 2*pi, length.out=11), function(a) c(cos(a), sin(a))))
set.seed(17)
radii <- (1 + rgamma(dim(circle)[1]-1, 3, 3))
radii <- c(radii, radii[1])  # Closes the loop
xy <- circle * radii
#
# Compute principal axes.
#
i.xy <- inertia(xy)
axes <- principal.axes(i.xy$Inertia)
sign <- sign(det(axes))
#
# Plot barycenter and principal axes.
#
plot(xy, bty="n", xaxt="n", yaxt="n", asp=1, xlab="x", ylab="y",
     main="A random polygon\nand its principal axes", cex.main=0.75)
polygon(xy, col="#e0e0e080")
arrows(rep(i.xy$Barycenter[1], 2), 
       rep(i.xy$Barycenter[2], 2),
       -axes[1,] + i.xy$Barycenter[1],     # The -signs make the first axis .. 
       -axes[2,]*sign + i.xy$Barycenter[2],# .. point to the right or down.
       length=0.1, angle=15, col=c("#e02020", "#4040c0"), lwd=2)
points(matrix(i.xy$Barycenter, 1, 2), pch=21, bg="#404040")

