Calculating mean value of a 2D array as a function of distance from the center in Python

Question

I'm trying to calculate the mean value of a quantity(in the form of a 2D array) as a function of its distance from the center of a 2D grid. I understand that the idea is that I identify all the array elements that are at a distance R from the center, and then add them up and divide by the number of elements. However, I'm having trouble actually identifying an algorithm to go about doing this.

I have attached a working example of the code to generate the 2d array below. The code is for calculating some quantities that are resultant from gravitational lensing, so the way the array is made is irrelevant to this problem, but I have attached the entire code so that you could create the output array for testing.

import numpy as np
import multiprocessing
import matplotlib.pyplot as plt
n = 100 # grid size

c = 3e8
G = 6.67e-11
M_sun = 1.989e30
pc = 3.086e16  # parsec
Dds = 625e6*pc          
Ds = 1726e6*pc #z=2
Dd = 1651e6*pc #z=1

FOV_arcsec = 0.0001
FOV_arcmin = FOV_arcsec/60.   


pix2rad = ((FOV_arcmin/60.)/float(n))*np.pi/180.
rad2pix = 1./pix2rad

Renorm = (4*G*M_sun/c**2)*(Dds/(Dd*Ds))
#stretch = [10, 2]




# To create a random distribution of points
def randdist(PDF, x, n):
    #Create a distribution following PDF(x). PDF and x
    #must be of the same length. n is the number of samples
    fp = np.random.rand(n,)
    CDF = np.cumsum(PDF)
    return np.interp(fp, CDF, x)


def get_alpha(args):    
    zeta_list_part, M_list_part, X, Y = args
    alpha_x = 0
    alpha_y = 0
    for key in range(len(M_list_part)):
        z_m_z_x = (X - zeta_list_part[key][0])*pix2rad
        z_m_z_y = (Y - zeta_list_part[key][1])*pix2rad
        alpha_x += M_list_part[key] * z_m_z_x / (z_m_z_x**2 + z_m_z_y**2)
        alpha_y += M_list_part[key] * z_m_z_y / (z_m_z_x**2 + z_m_z_y**2)
    return (alpha_x, alpha_y)

if __name__ == '__main__':
    # number of processes, scale accordingly
    num_processes = 1 # Number of CPUs to be used
    pool = multiprocessing.Pool(processes=num_processes)
    num = 100 # The number of points/microlenses
    r = np.linspace(-n, n, n)
    PDF = np.abs(1/r)
    PDF = PDF/np.sum(PDF)    # PDF should be normalized
    R = randdist(PDF, r, num)
    Theta = 2*np.pi*np.random.rand(num,)
    x1= [R[k]*np.cos(Theta[k])*1 for k in range(num)]
    y1 = [R[k]*np.sin(Theta[k])*1 for k in range(num)]
    # Uniform distribution
    #R = np.random.uniform(-n,n,num)
    #x1= np.random.uniform(-n,n,num)
    #y1 = np.random.uniform(-n,n,num)
    zeta_list = np.column_stack((np.array(x1), np.array(y1))) # List of coordinates for the microlenses 
    x = np.linspace(-n,n,n)
    y = np.linspace(-n,n,n)
    X, Y = np.meshgrid(x,y)
    M_list = np.array([0.1 for i in range(num)])
    # split zeta_list, M_list, X, and Y
    zeta_list_split = np.array_split(zeta_list, num_processes, axis=0)
    M_list_split = np.array_split(M_list, num_processes)
    X_list = [X for e in range(num_processes)]  
    Y_list = [Y for e in range(num_processes)]

    alpha_list = pool.map(
            get_alpha, zip(zeta_list_split, M_list_split, X_list, Y_list))
    alpha_x = 0
    alpha_y = 0
    for e in alpha_list:
        alpha_x += e[0] 
        alpha_y += e[1] 

alpha_x_y = 0
alpha_x_x = 0
alpha_y_y = 0
alpha_y_x = 0
alpha_x_y, alpha_x_x = np.gradient(alpha_x*rad2pix*Renorm,edge_order=2)
alpha_y_y, alpha_y_x = np.gradient(alpha_y*rad2pix*Renorm,edge_order=2) 
det_A = 1 - alpha_y_y - alpha_x_x + (alpha_x_x)*(alpha_y_y) - (alpha_x_y)*(alpha_y_x)
abs = np.absolute(det_A)
I = abs**(-1.)  
O = np.log10(I+1)
plt.contourf(X,Y,O,100)

The array of interest is O, and I have attached a plot of how it should look like. It can be different based on the random distribution of points.

What I'm trying to do is to plot the mean values of O as a function of radius from the center of the grid. In the end, I want to be able to plot the average O as a function of distance from center in a 2d line graph. So I suppose the first step is to define circles of radius R, based on X and Y.

def circle(x,y):
    r = np.sqrt(x**2 + y**2)
    return r 

Now I just have to figure out a way to find all the values of O, that have the same indices as equivalent values of R. Kinda confused on this part and would appreciate any help.

Answer 1

You can find the geometric coordinates of a circle with center (0,0) and radius R as such:

phi = np.linspace(0, 1, 50)
x = R*np.cos(2*np.pi*phi)
y = R*np.sin(2*np.pi*phi)

these values however will not fall on the regular pixel grid but in between.

In order to use them as sampling points you can either round the values and use them as indexes or interpolate the values from the near pixels.

Attention: The pixel indexes and the x, y are not the same. In your example (0,0) is at the picture location (50,50).