1

I'm looking to calculate the area within a contour of an unknown function. Values of my function are held within a Numpy array, which when plotted looks like this: enter image description here

I would like to calculate the cross-sectional area at different heights, for example 50% of the peak height.

Also, sometimes there is noise around the feature and therefore other peaks, meaning that a contour at a specific height can lead to more than one area being created. I'm only interested in the area within the contour that surrounds the "peak of interest", if looking from above.

The "peak of interest" is the one shown with the red dot at it's apex in the diagram above. I have an algorithm that identifies the peak of interest, so I always know which particular peak I'm interested in.

Can anyone help me with this? Many thanks

I found this question, but I can't get the given example to run. I'm having trouble with the line vs = contour.get_paths()[0].vertices, the list index is out of range.

I would also still need to adapt the code to make sure that the area encloses the peak location

This is the section of code that tries and find the areas levels = [0.5, 0.75] 

calc_levels = (radial[x_peak,y_peak])*np.transpose(levels)
cs = plt.contour(Z,X_,radial, levels = calc_levels,colors='k')
plt.clabel(cs, inline=2, fontsize=10)
for loop in range(len(levels)):
    vs = None
    contour_ = None
    contour_ = cs.collections[loop]
    vs = contour_.get_paths()[0].vertices
    ax.scatter(vs[:,0],vs[:,1],calc_levels[loop],c='red',s=20)
    # Compute area enclosed by vertices
    a = area(vs)
    features['Radial_Area_' + str(levels[loop]*100)] = area(vs)
    print("r = {}: a = {}".format(levels[loop],a))

where

def area(vs):
    a = 0
    x0,y0 = vs[0]
    for [x1,y1] in vs[1:]:
        dx = x1-x0
        dy = y1-y0
        a += 0.5*(y0*dx - x0*dy)
        x0 = x1
        y0 = y1
    return a

as taken for the example in the linked post in this question. However, it code always crashes on the line contour_ = cs.collections[loop], apparently the list index is out of range. I've tried every modification of code I can think of to get this working, but I'm now out of ideas.

Interestingly, if I run the code line by line then it doesn't crash and does give me output, even though it is an improbably small value

Interestingly, I eventually got this code working on the simple case, where only the peak of interest is enclosed. I have no idea what I changed from the above code to get it working. I would still like some help identifying which area to calculate in the case of multiple peaks being present. Please let me know if you are able to help me with this

Community
  • 1
  • 1
jlt199
  • 2,349
  • 6
  • 23
  • 43
  • A naive approach would be to count the number of elements of the array which are greater than the value at which you're taking the cross section. – Denziloe Mar 13 '17 at 20:34
  • How do you ensure it's in the contour surrounding the peak and not in another contour due to noise? – jlt199 Mar 13 '17 at 22:12
  • I don't understand the question in that case. What's "the" peak? If there are multiple cross sections then there are multiple peaks. – Denziloe Mar 13 '17 at 22:47
  • "the peak" is the one shown in the image, identified with the red dot. Now suppose the dataset is noisy and has several other large peaks, I'm still only interested in the contour that surrounds the peak of interest (the one with the red dot) - I've got a whole other algorithm to identify that particular peak, so let's assume that it works and the correct peak is always identified correctly. Now I want the area of the contour which is at 50% peak height and surrounds the peak of interest. – jlt199 Mar 13 '17 at 23:05

0 Answers0