All of the data collected by the United States during its decennial census is freely available from the Census Bureau. The data that we used throughout this paper comes from the 2010 data set. Our redistricting approach is driven primarily by the population density function, \(\rho\), of a particular state.

The Census Bureau provides population data down to the resolution of a single census tract. The Bureau also provides the geographic shapes of each census tract as they were during the collection of the data. If we let \(T\) be a census tract, \(A(T)\) be the area of the census tract (as projected upon the GCS North American 1983 [arcGIS]), and \(p(T)\) be the population of the tract, we define the population density function as \[\rho(\phi,\theta) = p(T)/A(T),\qquad (\phi,\theta)\in T\] The density function, \(\rho\), is a piecewise constant function defined within the boundaries of a single state. For definiteness, we assume that \(\phi\in [0,\pi]\) is latitude (with \(0\) at the North Pole) and \(\theta \in (-\pi,\pi]\) is longitude (with \(0\) at the Prime Meridian and \(\theta\) increasing to the east). Though not a sphere, we will assume a spherical approximation of the earth with radius \(R=3958.755\) miles [Moritz].

To facilitate numerical algorithms, we further discretize \(\rho\) using a uniform grid in the latitude/longitude domain as in the figure below. For simplicity, we assume that the state that we are redistricting is bounded by a spherical patch (some rectangle in the \((\phi,\theta)\)-plane) that is completely contained in both the northern hemisphere and in the western hemisphere. Hence, the latitudes spanning the state are contained in the interval \([\phi_{\min},\phi_{\max}] \subset [0,\pi/2]\) and the longitudes are contained in the interval \([\theta_{\min},\theta_{\max}] \subset (-\pi, 0]\). I.e. \(\phi_{\min}\) is the northern most latitude of the state and \(\theta_{\min}\) is the western most longitude of the state.

We discretize the spherical patch as \begin{eqnarray*} \phi_i & = & \phi_{\min} + \frac{\phi_{\max}-\phi_{\min}}{M},~i=0,1,\ldots M \\ \theta_j & = & \theta_{\min} + \frac{\theta_{\max}-\theta_{\min}}{N},~j=0,1,\ldots N. \end{eqnarray*}

**A uniform latitude/longitude grid on the surface of a sphere**

We then generate a discrete density function on the patch \(\rho_{ij} = \rho(x_i,y_j).\)

In order to reconstitute the population of a particular region, it is necessary to know the area of of each grid square in the discretization. Recall that each grid square of latitude and longitude represents a small patch on the surface of the earth. Furthermore, though the grid squares are uniform in the latitude and longitude domain, the actual surface area represented by each square depends upon its latitude again see this figure and this figure. We approximate the area of the patch with upper-left hand corner at \((x_i, y_j)\) by assuming that the earth is a sphere and using a spherical surface integral.

**A flat projection of a uniform latitude longitude grid**

For later reference, the colored dots indicate a Moore-type neighborhood set. The red dot is the center of the neighborhood. The yellow dots are the neighbors of the center dot. Every dot is the center of an associated Moore neighborhood.

First we convert \((x_i,y_j)\) to spherical coordinates. \begin{eqnarray*} \phi_i & = & (90^\circ - x_i)\cdot \frac{\pi}{180^\circ} \\ \theta_j & = & y_j \cdot \frac{\pi}{180^\circ}. \end{eqnarray*} Then, \begin{eqnarray*} A(x_i,y_j) & = & \int_{\theta_j}^{\theta_{j+1}} \int_{\phi_i}^{\phi_{i+1}} R^2 \sin\phi \, d\phi d\theta\\ & = & R^2(\theta_{j+1} - \theta_j) (\cos(\phi_{i+1}) - \cos(\phi_i)). \end{eqnarray*} It follows immediately that the population of a particular grid square is well-approximated by \(\rho_{ij} A(x_i,y_j)\).