Finding the markersize/radius of circles, some of which overlap, needed to cover a given percentage of a polygon in Python

Question

I am trying to figure out the markersize needed in python so that the green circles in the attached figure cover a specific percentage, say 20%, of the area bounded by the blue polygon.

The green circles represent epicenters of earthquakes, with lon and lat coordinates. Markersize and radius are somewhat synonymous in this scenario, but not exactly the same because of the way markersize works.

I thought I was able to brute force this using pure math, but ended up not accounting for overlap of the circles. I then tried using a pixel method that essentially checks each pixel within the polygon to see if it is green vs any other color, but then realized that I can't figure that out because I don't know the markersize needed. I tried a monte carlo method as well, but ran into a similar issue.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
import os
from mpl_toolkits.basemap import Basemap
import sys
from bs4 import BeautifulSoup
import geopandas as gpd
import zipfile
import xml.etree.ElementTree as ET
import pickle
import chardet
import csv
import math
import random
from math import pi
from shapely.geometry import Point, Polygon

# Call in text files that include your list of earthquakes
#file_path = os.path.join("path", "to", "total_san_andreas.txt")
cat = np.loadtxt(r"C:/Users/efisher/Documents/WUS M5+ Research/total_san_andreas.txt")
#file_path = os.path.join("path", "to", "san_andreas_polygon.txt")
SA_polygon = np.loadtxt(r"C:/Users/efisher/Documents/WUS M5+ Research/san_andreas_polygon.txt")

#define the year cut-off for the before and after catalogs
def cat_date(cat, ba="after", year=2000):
    if ba == "after":
        newcat = [entry for entry in cat if float(entry[0]) >= year]
        print(f"My after catalog has {len(newcat)} entries occurring on or after the year {year}!")
    elif ba == "before":
        newcat = [entry for entry in cat if float(entry[0]) < year]
        print(f"My after catalog has {len(newcat)} entries occurring before the year {year}!")
        
    return newcat
        
#split your input earthquake file into before and after catalogs
def splitcat(cat, year):
    aftercat = cat_date(cat, year=year)
    beforecat = cat_date(cat, ba="before", year=year)
    return aftercat, beforecat

# extract columns
aftercat, beforecat = splitcat(cat, 2000)
#get coordinates of each of the earthquakes within your polygon
def getcoords(cat):
    lat = np.array([entry[1] for entry in cat])
    lon = np.array([entry[2] for entry in cat])
    in_poly = Path(SA_polygon, closed=True).contains_points(np.column_stack((lat, lon)))
    lon_in, lat_in = lon[in_poly], lat[in_poly]
    return lon_in, lat_in


blons, blats = getcoords(beforecat)
alons, alats = getcoords(aftercat)

year = cat[:, 0]
year2 = [row[0] for row in cat]
lat = cat[:, 1]
lon = cat[:, 2]
mag = cat[:, 3]

#print(type(lat))

SA_path = Path(SA_polygon)

# Find the area within your polygon
def polygon_area(coords):
    R = 6371 #earth radius in km
    lat = [coord[0]for coord in coords]
    lon = [coord[1]for coord in coords]
    area = 0
    #converts the area to square kilometers
    for pos in range(len(lat)-1):
        dlat = math.radians(lat[pos+1] - lat[pos])
        dlon = math.radians(lon[pos+1] - lon[pos])
        a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat[pos])) * math.cos(math.radians(lat[pos+1])) * math.sin(dlon/2) * math.sin(dlon/2)
        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
        area += (R**2 * (1 - math.cos(c)))/2
    return area

area = polygon_area(SA_polygon)
rounded_area = round(area) #rounds to the nearest integer value
print("Your Area of Interest has an area of", rounded_area, "km^2")

# Below is how you can find he names of relevant faults within the area of interest.
# this allows you to look up and store only the faults you wish to map that are within the pickle file
# In this case, I include the the San andreas and several of its larger/smaller branch faults.
def getfaultcoords(faultname, fcd):
    lats = []
    lons = []
    for key in fcd:
        if faultname in key:
            for coords in faultcoords[key]:
                cs = coords.split(",")
                lon = float(cs[0])
                lat = float(cs[1])
                lats.append(lat)
                lons.append(lon)
    return lats, lons

with open('faultcoords_pickle.pkl', 'rb') as f:
    faultcoords = pickle.load(f)
    
relfaults = ["San Andreas", "Garlock", "Elsinore", "San Jacinto", "Calaveras", "Hayward", "Malibu", "Big Pine", 
             "Santa Cruz", "Oak Ridge", "Imperial", "Laguna Salada", "Rodgers Creek", "San Gregorio", 
             "Mendocino", "Maacama", "Santa Ynez", "Oceanic", "Sierra Madre", "Lockhart", "Calico", 
             "Camp Rock", "Harper", "Pisgah", "Pinto", "Owens", "Wolf", "Wheeler Ridge", "Bartlett", "Whale Gulch", "Honeydew"]

rfs = {}
for fault in relfaults:
    lats, lons = getfaultcoords(fault, faultcoords)
    rfs[fault] = {"lats": lats, "lons": lons}


#create basemap object and plot
fig, ax = plt.subplots()

#Give the lat/lon coords based on corners of the map you are trying to create
#can set resolution to h or l (high/low)
#EPSG code is for a mercator projection taken from the "European Petroleum Survey Group"
#can change the image quality by changing the pixels. The higher the pixels, the longer it will take to process
m = Basemap(projection='merc', llcrnrlat=30, urcrnrlat=42, llcrnrlon=-126, urcrnrlon=-112, resolution='h', ax = ax, epsg=3857)
m.arcgisimage(service='World_Imagery', xpixels=1000, verbose=False)

#draw coastlines, countriy lines, states, parallels, etc.
m.drawcoastlines()
m.drawcountries()
m.drawstates()
m.drawparallels(np.arange(30, 42, 3), labels=[True, False, False, False])
m.drawmeridians(np.arange(-126, -112, 3.5), labels=[False, False, False, True])

# plot data
#converts lats/lons in SA_polygon to x,y coordinates. Plots as blue line (b-)
plons, plats = m(SA_polygon[:, 1], SA_polygon[:, 0])
m.plot(plons, plats, 'b-', markersize=1)

#converts lats and lons in the before/after catalogs to x,y coordinates
blons, blats = m(blons, blats)
alons, alats = m(alons, alats)

for fault in rfs:
    lons = rfs[fault]['lons']
    lats = rfs[fault]['lats']
    mlons, mlats = m(lons, lats)
    m.plot(mlons[:-1], mlats[:-1], 'o', color='orange', markersize=0.2)

# find the markersize needed for the "before" earthquakes to cover a specified percentage of the map/polygon area
percentage = int(input("What percentage of the polygon/area of interest would you like covered by before-cat earthquakes? "))
def get_markersize(rounded_area, percentage):
    blons, blats = getcoords(beforecat)
    polygon_area = rounded_area
    beforecat_area = len(blons) / len(beforecat) * polygon_area
    markersize = (((rounded_area) * percentage)/100)/len(beforecat)
    return markersize

markersize = get_markersize(rounded_area, percentage)
radius = (markersize/pi)**(0.5)
rounded_radius = round(radius, 2)
print("The markersize needed to cover ", percentage, " percent of the map/polygon is ", rounded_radius)



# find the markersize needed for the "before" earthquakes to cover a specified percentage of the map/polygon area
def monte_carlo(blats, blons, percentage):
    n = 1000000 # number of random points to generate
    count = 0
    for i in range(n):
        # Generate random points within the SA_polygon
        lon = random.uniform(min(blons), max(blons))
        lat = random.uniform(min(blats), max(blats))
        point = Point(lon, lat)
        if point.within(SA_path):
            count += 1
    fraction_inside = count / n
    radius = math.sqrt(fraction_inside * rounded_area / pi / percentage)
    rounded_radius = round(radius, 2) #round to two decimal places
    print(f"The radius of the circles needed to cover {percentage}% of the area is approximately {rounded_radius} km")
    return rounded_radius

rounded_radius = monte_carlo(blats, blons, percentage)
#print(f"The radius of the circles needed to cover {percentage}% of the area is approximately {rounded_radius} km")
#points = monte_carlo(point)




#plots coordinates of before and after earthquake catalogs as green circle (go) and red + (r+) respectively
m.plot(blons, blats, 'go', markersize=rounded_radius)
m.plot(alons, alats, 'r+', markersize=5)
# ax.set_aspect('equal')

# Add a title, axes, and legend to the figure.
# Labelpad moves the labels. Postive values move down on x-axis and left on y-axis, and vice versa.
ax.set_title("M5+ Seismicity along the San Andreas Fault and it's Branch Faults")
ax.set_xlabel("Longitude", labelpad=10)
ax.set_ylabel("Latitude", labelpad=40, rotation=90, va='center', ha='right')

# Create custom legend entries with different shapes and colors
green = plt.Line2D([], [], color='green', marker='o', linestyle='None', label='Earthquakes Before the Year 2000')
red= plt.Line2D([], [], color='red', marker='+', linestyle='None', label='Earthquakes Year 2000-Present')
orange = plt.Line2D([], [], color='orange', marker='_', linestyle='None', label='Relevant Faults')
blue = plt.Line2D([], [], color='blue', marker='_', linestyle='None', label='San  Andreas Polygon')

# Add a legend to the plot with custom entries
plt.legend(handles=[green, red, orange, blue],
           labels=['Earthquakes Before the Year 2000', 'Earthquakes Year 2000-Present', 'Relevant Faults', 'San Andreas Polygon'],
           bbox_to_anchor=(1.05, 0.5), loc='center left')
#ax.legend(loc='center right', title="Legend", "line1", "line2", "line3", "line4")
plt.show()