How can I count number of points between 2 contours in Python

Question

I am using matplotlib.pyplot to interpolate my data and create contours. Following this answer/example (about how to calculate area within a contour), I am able to get the vertices of a contour line. Is there a way to use that information, i.e., the vertices of a line, to count how many points fall between two given contours? These points will be different from the data used for deriving the contours.

Answer 1

I am not sure if I understand what points do you want to check, but, if you have the line vertices (two points) and want to check if the third point falls in between the two, you can take a simple (not efficient) approach and calculate the area of the triangle formed by the three. If the area is 0 then the point falls on the same line. Also, you can calculate the distance between the points and see if the point is in between on the line or outside (on the extended) line.

Hope this helps!

Answer 2

Usually, you do not want to reverse engineer your plot to obtain some data. Instead you can interpolate the array that is later used for plotting the contours and find out which of the points lie in regions of certain values.

The following would find all points between the levels of -0.8 and -0.4, print them and show them in red on the plot.

import numpy as np; np.random.seed(1)
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf

X, Y = np.meshgrid(np.arange(-3.0, 3.0, 0.1), np.arange(-2.4, 1.0, 0.1))
Z1 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = mlab.bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = 10.0 * (Z2 - Z1)

points = np.random.randn(15,2)/1.2
levels = [-1.2, -0.8,-0.4,-0.2]

# interpolate points
f = Rbf(X.flatten(), Y.flatten(), Z.flatten()) 
zi = f(points[:,0], points[:,1])
# add interpolated points to array with columns x,y,z
points3d = np.zeros((points.shape[0],3))
points3d[:,:2] = points
points3d[:,2] = zi
# masking condition for points between levels
filt = (zi>levels[1]) & (zi <levels[2])
# print points between the second and third level
print(points3d[filt,:])

### plotting
fig, ax = plt.subplots()

CS = ax.contour(X, Y, Z, levels=levels)
ax.clabel(CS, inline=1, fontsize=10)
#plot points between the second and third level in red:
ax.scatter(points[:,0], points[:,1], c=filt.astype(float), cmap="bwr" )

plt.show()

Answer 3

Ampere's law of magnetic fields can help here - although it can be computationally expensive. This law says that the path integral of a magnetic field along a closed loop is proportional to the current inside the loop.

Suppose you have a contour C and a point P (x0,y0). Imagine an infinite wire located at P perpendicular to the page (current going into the page) carrying some current. Using Ampere's law we can prove that the magnetic field produced by the wire at a point P (x,y) is inversely proportional to the distance from (x0,y0) to (x,y) and tangential to the circle centered at point P passing through point (x,y). Therefore if the wire is located outside the contour the path integral is zero.

import numpy as np
import pylab as plt

# generating a mesh and values on it
delta = 0.1
x = np.arange(-3.1*2, 3.1*2, delta)
y = np.arange(-3.1*2, 3.1*2, delta)
X, Y = np.meshgrid(x, y)
Z = np.sqrt(X**2 + Y**2)

# generating the contours with some levels
levels = [1.0]
plt.figure(figsize=(10,10))
cs = plt.contour(X,Y,Z,levels=levels)

# finding vertices on a particular level
contours = cs.collections[0]
vertices_level = contours.get_paths()[0].vertices # In this example the 
shape of vertices_level is (161,2)
# converting points into two lists; one per dimension. This step can be optimized
lX, lY = list(zip(*vertices_level))
          
# computing Ampere's Law rhs 
def AmpereLaw(x0,y0,lX,lY):
    S = 0
    for ii in range(len(lX)-1):
        dx = lX[ii+1] - lX[ii]
        dy = lY[ii+1] - lY[ii]
        ds = (1/((lX[ii]-x0)**2+(lY[ii]-y0)**2))*(-(lY[ii]-y0)*dx+(lX[ii]-x0)*dy)
        if -1000 < ds < 1000: #to avoid very lare numbers when denominator is small
            S = S + ds
    return(S)
    
# we know point (0,0) is inside the contour
AmpereLaw(0,0,lX,lY)
# result: -6.271376740062852

# we know point (0,0) is inside the contour
AmpereLaw(-2,0,lX,lY)
# result:  0.00013279920934375876

You can use this result to find points inside one contour but outside of the other.