Calculating the fraction of the area of multiple squares overlapped by a circle

Question

This is a geometrical question based on a programming problem I have. Basically, I have a MySQL database full of latitude and longitude points, spaced out to be 1km from each other, corresponding to a population of people who live within the square kilometer around each point. I then want to know the relative fraction of each of those grids taken up by a circle of arbitrary size that overlaps them, so I can figure out how many people roughly live within a given circle.

Here is a practical example of one form of the problem (distances not to scale):

I am interested in knowing the population of people who live within a radius of point X. My database figures out that its entries for points A and B are close enough to point X to be relevant. Point A in this example is something like 40.7458, -74.0375, and point B is something like 40.7458, -74.0292. Each of those green lines from A and B to its grid edge represents 0.5 km, so that the gray circle around A and B each represent 1 km^2 respectively.

Point X is at around 40.744, -74.032, and has a radius (in purple) of 0.05 km.

Now I can easily calculate the red lines shown using geographic trig functions. So I know that the line AX is about .504 km, and the distance line BX is about .309 km, for whatever that gets me.

So my question is thus: what's a solid way for calculating the fraction of grid A and the fraction of grid B taken up by the purple circle inscribed around X?

Ultimately I will be taking the population totals and multiplying them by this fraction. So in this case, the 1 km^2 grid around corresponds to 9561 people, and the grid around B is 10763 people. So if I knew (just hypothetically) that the radius around X covered 1% of the area of A and 3% of the area of B, I could make a reasonable back-of-the-envelope estimate of the total population covered by that circle by multiplying the A and B populations by their respective fractions and just summing them.

I've only done it with two squares above, but depending on the size of the radius (which can be arbitrary), there may be a whole host of possible squares, like so, making it a more general problem:

In some cases, where it is easy to figure out that the square grid in question is 100% encompassed by the radius, it is in principle pretty easy (e.g. if the distance between AX was smaller than the radius around X, I know I don't have to do any further math).

Now, it's easy enough to figure out which points are within the range of the circle. But I'm a little stuck on figuring out what fractions of their corresponding areas are.

Thank you for your help.

Answer 1

I ended up coming up with what worked out to be a pretty good approximate solution, I think. Here is how it looks in PHP:

//$p is an array of latitude, longitude, value, and distance from the centerpoint 
//$cx,$cy are the lat/lon of the center point, $cr is the radius of the circle
//$pdist is the distance from each node to its edge (in this case, .5 km, since it is a 1km x 1km grid)

function sum_circle($p, $cx, $cy, $cr, $pdist) {
    $total = 0; //initialize the total
    $hyp = sqrt(($pdist*$pdist)+($pdist*$pdist)); //hypotenuse of distance

    for($i=0; $i<count($p); $i++) { //cycle over all points
        $px = $p[$i][0];    //x value of point
        $py = $p[$i][1];    //y value of point
        $pv = $p[$i][2];    //associated value of point (e.g. population)
        $dist = $p[$i][3];  //calculated distance of point coordinate to centerpoint

        //first, the easy case — items that are well outside the maximum distance
        if($dist>$cr+$hyp) { //if the distance is greater than circle radius plus the hypoteneuse
            $per = 0; //then use 0% of its associated value
        } else if($dist+$hyp<=$cr) { //other easy case - completely inside circle (distance + hypotenuse <= radius)
            $per = 1; //then use 100% of its associated value
        } else { //the edge cases
            $mx = ($cx-$px); $my = ($cy-$py); //calculate the angle of the difference
            $theta = abs(rad2deg(atan2($my,$mx)));
            $theta = abs((($theta + 89) % 90 + 90) % 90 - 89); //reduce it to a positive degree between 0 and 90
            $tf = abs(1-($theta/45)); //this basically makes it so that if the angle is close to 45, it returns 0, 
                                      //if it is close to 0 or 90, it returns 1
            $hyp_adjust = ($hyp*(1-$tf)+($pdist*$tf)); //now we create a mixed value that is weighted by whether the
                                                        //hypotenuse or the distance between cells should be used                                                   

            $per = ($cr-$dist+$hyp_adjust)/100; //lastly, we use the above numbers to estimate what percentage of 
                                                //the square associated with the centerpoint is covered
            if($per>1) $per = 1; //normalize for over 100% or under 0%
            if($per<0) $per = 0; 
        }
        $total+=$per*$pv;   //add the value multiplied by the percentage to the total
    }  
    return $total;
}

This seems to work and is pretty fast (even though it does use some trig on the edge cases). The basic logic is that when calculating edge cases, the two extreme possibilities is that the circle radius is either exactly perpendicular to the grid, or exactly at 45 degree angles from it. So it figures out roughly where between those extremes it falls and then uses that to figure out roughly what percentage of the grid square is covered. It gives plausible results in my testing.

For the size of the squares and circles I am using, this seems to be adequate?

I wrote a little application in Processing.js to try and help me work this out. Without explaining all of it, you can see how the algorithm is "thinking" by looking at this screenshot:

Basically, if the circle is yellow it means it has already figured out it is 100% in, and if it is red it is already quickly screened as 100% out. The other cases are the edge cases. The number (ranging from 0 to 1) under the dot is the (rounded) percentage of coverage calculated using the above method, while the number under that is the calculated theta value used in the above code.

For my purposes I think this approximation is workable.