How to calculate partial area-under-the-curve?

Hi, I'm trying to find the area under the curve (ABC) for a part of a graph. I use the "trapz" function, but this function calculates the ABC for an entire area below the selected part of the graph. Any tips on how I can calculate only part of it (not the whole part up to the x axis)? Please see the figure. Interested are to be calculate is in red.

Thanks, Jacqueline

Answer 1

[Star Strider]

Use trapz twice, once to integrate (A, B, C), and again to integrate (A, C).

Then subtract the (A, C) result from the (A, B, C) result.

It would definitely help to have the data, preferably as a .mat file.

[Jacqueline Rigatto]

x_32= [0;
    0.0500;
    5.0000;
    5.0500;
   10.1500;
   10.2000;
   12.9200;
   15.1600;
   17.0600;
   22.3700;
   28.8800;
   43.3600];

z_MHWS_32=[4.1526;
    4.0126;
    3.9426;
    3.8246;
    4.8826;
    4.8096;
    4.1466;
    2.7946;
    2.5326;
    1.8326;
    1.0296;
   -0.4854];


NM_18=0.82;
MHWS_18=1.5-NM_18;
MHWS_vector_18=repmat(MHWS_18,79,1);

figure(1)
plot (x_32,z_MHWS_32,'Green')
hold on
plot(MHWS_vector_18,'Blue')
hold off

x=x_32;  
z=z_MHWS_32; 

xi=5.05; 
xf=12.92; 

xRange = [xi,xf];
% find logical indices for range of interest
idl = x >= xRange(1) & x <= xRange(2);
% provide just values of interest to trapz
selected_area=trapz(x(idl),z(idl))

This is my code above. I used "trapz" between A and C to calculate the top part, but it calculated up to the x axis.

[Star Strider]

Running the posted code:

Unrecognized function or variable 'MHWS_18'.

Waiting for it to magickally appear!

.

[Jacqueline Rigatto]

Edited code. Thank you for letting me know

[Star Strider]

My pleasure!

Try this:

xRange = [xi,xf];
% find logical indices for range of interest
idl = x >= xRange(1) & x <= xRange(2);
% provide just values of interest to trapz
selected_area = trapz(x(idl),z(idl))
idlidx = find(idl);                                                                     % Numeric Indices
selected_baseline = interp1(xRange,z_MHWS_32([idlidx(1),idlidx(end)]),x_32(idl));       % Interpolate
baseline_area = trapz(x(idl),selected_baseline)                                         % X-Axis To Baseline Area
selected_area_above_baseline = selected_area - baseline_area                            % Desired Result

producing:

selected_area = 
   34.6261

baseline_area = 
   31.3667

selected_area_above_baseline = 
    3.2594

The additional code usess the ends of the ‘xRange’ variable to interpolate the intervening ‘x_32’ values, then integrates them and does the subtraction to find the area between the baseline and the desired area. The result appears (to me) to be reasonable.

[Jacqueline Rigatto]

Thank you for your help, Star Strider!

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Jacqueline Rigatto]

Star Strider, consider the same graphic and the same code, but in a different location (figure in the attachment). I calculated the red area but it was negative using:

xi=15.16; 
xf=32;

producing:

selected_area =

   25.9669


baseline_area =

   28.4773


selected_area_above_baseline =

   -2.5104
   

Is this negative value correct or do I have to use the module?

[Star Strider]

Using the ‘new’ ‘xi’ and ‘xf’:

selected_baseline = interp1(xRange,MHWS_vector_18([idlidx(1),idlidx(end)]),x_32(idl));  % Interpolate
                                       ↑ ← CHANGE VECTOR

producing:

selected_area = 
   25.9669

baseline_area = 
    9.3296

selected_area_above_baseline = 
   16.6373

Here, ‘baseline_area’ is between the x-axis and the blue horizontal line, not with respect to different points on the green curve, as in the original problem.

[Jacqueline Rigatto]

Thanks again for your help! helped a lot, Star Strider!

[Star Strider]

As always, my pleasure!

[Jacqueline Rigatto]

Last case (attached): suppose that the two areas outside the rectangle have already been calculated; to know the area of ​​the red rectangle I calculate the whole area and then do two subtractions? The total area I use xi = 0.0500 and xf = 32 using the routine below and the result will come out in "selected_area_above_baseline" or "selected_area"?

xRange = [xi,xf];
% find logical indices for range of interest
idl = x >= xRange(1) & x <= xRange(2);
% provide just values of interest to trapz
selected_area = trapz(x(idl),z(idl))
idlidx = find(idl);                                                                     % Numeric Indices
selected_baseline = interp1(xRange,MHWS_vector_18([idlidx(1),idlidx(end)]),x_32(idl));  % Interpolate
baseline_area = trapz(x(idl),selected_baseline)                                         % X-Axis To Baseline Area
selected_area_above_baseline = selected_area - baseline_area  

[Star Strider]

The red rectangle is significantly different from the other areas, not bordered by both, and not related to the first area, so it would not be accurate to calculate it from those. It would be more closely related to the second area calculation..

Perhaps:

idx = x_32 <= 15.16;
nidx = find(idx);
RedRectangleArea = trapz(x_32(idx), ones(size(x_32(idx)))*z_MHWS_32(nidx(end))) - trapz(x_32(idx), MHWS_vector_18(idx))

producing:

RedRectangleArea =
   32.0573

.

[Jacqueline Rigatto]

Thank you so much for all your help and sorry for bothering you

[Star Strider]

As always, my pleasure!

Not a bother at all! It’s definitely an interesting set of problems!

[Jacqueline Rigatto]

Thank you very much from the heart!

Below is a figure and a code of how to calculate each area (1 to 3). My results:

Area 1: 0.9655;

Area 2: 10.5941;

Area 3: 16.6373;

Total area (sum of areas 1, 2 and 3): 33.4274

If you add the areas from 1 to 3 in the "hand", give 28.1969, but if you do it by the code, give 33.4274. It was not to give the same value? Is there a way to find the total area using the trapz?

Area 3 and the total area is the same code.

clc; clear all; close all


x_32= [0;
    0.0500;
    5.0000;
    5.0500;
   10.1500;
   10.2000;
   12.9200;
   15.1600;
   17.0600;
   22.3700;
   28.8800;
   43.3600];

z_MHWS_32=[4.1526;
    4.0126;
    3.9426;
    3.8246;
    4.8826;
    4.8096;
    4.1466;
    2.7946;
    2.5326;
    1.8326;
    1.0296;
   -0.4854];

NM_18=0.82;
MHWS_18=1.5-NM_18;
MHWS_vector_18=repmat(MHWS_18,79,1);

figure(1)
plot (x_32,z_MHWS_32,'Green')
hold on
plot(MHWS_vector_18,'Blue')
hold off

x=x_32;  
z=z_MHWS_32; 

%%%%%%%%%%%%%%%% 1
% xi=10.15; 
% xf=15.16; 

%%%%%%%%%%%%%%%% 2
% xi=10.15; 
% xf=15.16; 

%%%%%%%%%%%%%%%% 3
% xi=15.16; 
% xf=32;

%%%%%%%%%%%%%%%% all areas (1, 2 and 3)
xi=10.15;
xf=32;
%%%%%%%%%%%%%%%%%%%%%% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% xRange = [xi,xf];
% % find logical indices for range of interest
% idl = x >= xRange(1) & x <= xRange(2);
% % provide just values of interest to trapz
% selected_area = trapz(x(idl),z(idl))
% idlidx = find(idl);                                                                     % Numeric Indices
% selected_baseline = interp1(xRange,z_MHWS_32([idlidx(1),idlidx(end)]),x_32(idl));       % Interpolate
% baseline_area = trapz(x(idl),selected_baseline)                                         % X-Axis To Baseline Area
% selected_area_above_baseline = selected_area - baseline_area       


%%%%%%%%%%%%%%%%%%%%%% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% xRange = [xi,xf];
% idx = x >= xRange(1) & x <= xRange(2);
% nidx = find(idx);
% RedRectangleArea = trapz(x_32(idx), ones(size(x_32(idx)))*z_MHWS_32(nidx(end))) - trapz(x_32(idx), MHWS_vector_18(idx))


%%%%%%%%%%%%%%%%%%%%%% 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% xRange = [xi,xf];
% % find logical indices for range of interest
% idl = x >= xRange(1) & x <= xRange(2);
% % provide just values of interest to trapz
% selected_area = trapz(x(idl),z(idl))
% idlidx = find(idl);                                                                     % Numeric Indices
% selected_baseline = interp1(xRange,MHWS_vector_18([idlidx(1),idlidx(end)]),x_32(idl));  % Interpolate
% baseline_area = trapz(x(idl),selected_baseline)                                         % X-Axis To Baseline Area
% selected_area_above_baseline = selected_area - baseline_area      

%%%%%%%%%%%%%%%%%%%%%% all areas (1, 2 and 3) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xRange = [xi,xf];
% find logical indices for range of interest
idl = x >= xRange(1) & x <= xRange(2);
% provide just values of interest to trapz
selected_area = trapz(x(idl),z(idl))
idlidx = find(idl);                                                                     % Numeric Indices
selected_baseline = interp1(xRange,MHWS_vector_18([idlidx(1),idlidx(end)]),x_32(idl));  % Interpolate
baseline_area = trapz(x(idl),selected_baseline)                                         % X-Axis To Baseline Area
selected_area_above_baseline = selected_area - baseline_area      

[Star Strider]

These assignments after figure(1):

x=x_32;  
z=z_MHWS_32; 

[gmax,gmix] = max(z_MHWS_32);
idxrng = gmix:numel(z_MHWS_32);                                 % Vector Indices From ‘maxg’ To End
xv = x_32(idxrng);                                              % ‘x’ Vector
bv = MHWS_vector_18(idxrng);                                    % Blue Line Vector
gv = z_MHWS_32(idxrng);                                         % Green Line Vector
bgi = interp1(bv-gv, xv, 0);                                    % ‘x’-Coordinate Of The Intersection Of Blue & Green Lines

gintrpv = linspace(xv(1), bgi);                                 % Interpolation Vector Spanning Range (10.15, 32.22)
gintrp = interp1(xv, gv, gintrpv);                              % Interpolated Green Line
bintrp = interp1(xv, bv, gintrpv);                              % Interpolated Blue Line
GreenArea = trapz(gintrpv, gintrp);                             % From x-Axis To Green Line
BelowBlueLine = trapz(gintrpv, bintrp);                         % From x-Axis To Blue Line
AreaGreenGreaterThanBlue = GreenArea - BelowBlueLine

produces:

AreaGreenGreaterThanBlue =
    34.0175

So the 33.4274 value is closer. This latest computation uses a significantly greater number of points ( because the intersection of the blue and green lines is not exactly at an index of ‘x_32’ ) so interpolation is necessary to calculate the areas.

[Jacqueline Rigatto]

In my case, is it better to use trapz or interpolation?

[Star Strider]

I used both here. I would interpolate first to get increased resolution of the data (use the default 'linear' interpolation method, since the data do not deviate much from being linear in the various segments), then do whatever other processing would be necessary, and then use trapz to get the areas. The areas using the interpolated data would be more precise than the areas using the original data, so the values would differ slightly. The areas using the interpolated data will be more accurate, however the actual values may not differ much from those using the original data.

[Jacqueline Rigatto]

Thanks for all the help, Star Strider

[Star Strider]

As always, my pleasure!

[Jacqueline Rigatto]

Star Strider, what would areas 1, 2 and 3 look like if you do interpolation? The interpolation you did was for the 3 areas together?

[Star Strider]

The areas would likely not be significantly different from the original calculations, since all the calculations were with respect to specific points. Actually, several fo them (for example 15 Mar 2021 at 18:33, and 15 Mar 2021 at 23:37 ) used interpolation, the first to create a ‘lower border’ vector the same size as the upper border of that area, and that last for the same reason I used interpolation for the entire area because the end-point was not a specific index, so the code needed to create one.

[Jacqueline Rigatto]

Thanks again for your help, Star Strider

[Star Strider]

As always, my pleasure!

[Jacqueline Rigatto]

Hi Star Strider, can you take one more question? On 18 Mar 2021 at 20:00, how do i calculate area 2 with initial x and final x? Will be the same on 18 Mar 2021 at 0:25?

[Star Strider]

I have forgotten where we were on this.

The easiest way would be:

Area_2 = (15.16 - 10.15)*(2.7946-0.68)

producing:

Area_2 = 
    10.5941

If you prefer a different approach, I will have to remember what we were doing with this, and go back and reconstruct the trapz calls.

[Jacqueline Rigatto]

Thanks for all the help, Star Strider

Answer 2

[John D'Errico]

After all of those comments, to be honest, sorry, but you were both working far too hard on a moderately simple problem.

Simplest is to just create a polygon, then use polyarea of that object. Alternatively, create a polyshape object. Again, compute the area of said polygonal region.

help polyshape

  POLYSHAPE Create polyshape object
 
  A polyshape is a 2-D polygon that can consist of one or more solid regions,
  and can have holes. polyshape objects are created from the x- and 
  y-coordinates that describe the polygon's boundaries.   
 
  Boundaries that make up a polyshape should have no self-intersections or 
  intersections with other boundaries, and all boundaries should be properly 
  nested. Otherwise, the polyshape is ill-defined, and can lead to inaccurate
  or unexpected results.
 
  PG = POLYSHAPE(X, Y) creates a polyshape object from 2-D vertices given
  by two vectors X and Y. X and Y must have the same size.
 
  PG = POLYSHAPE(P) creates a polyshape from 2-D points contained in the
  2-column array P. 
 
  PG = POLYSHAPE({X1, X2, ..., Xn}, {Y1, Y2, ..., Yn}) creates a polyshape 
  with n boundaries from two cell arrays. Each cell array contains n vectors, 
  where each Xi contains the x-coordinates and Yi contains the y-coordinates
  of the i-th boundary.
 
  PG = POLYSHAPE(..., 'SolidBoundaryOrientation', DIR) specifies the 
  direction convention for determining solid versus hole boundaries. DIR 
  can be one of the following:
   'auto' (default) - Automatically determine direction convention
   'cw' - Clockwise vertex direction defines solid boundaries
   'ccw' - Counterclockwise vertex direction defines solid boundaries
 
  This name-value pair is typically only specified when creating a polyshape 
  from data that was produced by other software that uses a particular 
  convention.  
 
  PG = POLYSHAPE(..., 'Simplify', tf) specifies how ill-defined polyshape 
  boundaries are handled. tf can be one of the following:
   true (default) - Automatically alter boundary vertices to create a 
                    well-defined polygon
   false - Do not alter boundary vertices even though the polyshape is 
           ill-defined. This may lead to inaccurate or unexpected results.
 
  PG = POLYSHAPE(..., 'KeepCollinearPoints', tf) specifies how to treat 
  consecutive vertices lying along a straight line. tf can be one of the 
  following:
    true  - Keep all collinear points as vertices of PG.
    false - Remove collinear points so that PG contains the fewest number
            of necessary vertices.
  The value of the 'KeepCollinearPoints' parameter is carried over when the
  ADDBOUNDARY, INTERSECT, SIMPLIFY, SUBTRACT, UNION, and XOR functions are
  used with input PG.
 
  PG = POLYSHAPE() creates a polyshape object with no vertices.
 
  Polyshape properties:
    Vertices      - 2-D polyshape vertices
    NumRegions    - Number of regions in the polyshape
    NumHoles      - Number of holes in the polyshape
 
  Methods for modifying a polyshape:
    addboundary     - Add boundaries to a polyshape
    rmboundary	    - Remove boundaries in a polyshape
    rmholes         - Remove all the holes in a polyshape
    rmslivers       - Clean up degeneracies in a polyshape
    simplify        - Fix degeneracies and intersections in a polyshape 
 
  Methods for querying a polyshape:
    ishole          - Determine if a boundary is a hole 
    issimplified    - Determine if a polyshape is simplified
    numsides        - Find the total number of sides in a polyshape
    numboundaries   - Find the total number of boundaries in a polyshape 
 
  Methods for geometric information:
    area            - Find the area of a polyshape 
    boundary        - Get the x- and y-coordinates of a boundary 
    boundingbox	    - Find the bounding box of a polyshape 
    centroid	    - Find the centroid of a polyshape 
    convhull	    - Find the convex hull of a polyshape 
    holes           - Convert all holes in a polyshape into an array of polyshapes 
    isinterior	    - Determine if a point is inside a polyshape 
    nearestvertex   - Find the vertex of a polyshape nearest to a point
    overlaps        - Determine if two polyshapes overlap
    perimeter       - Get the perimeter of a polyshape 
    regions         - Put all regions in a polyshape into an array of polyshapes
    triangulation   - Construct a 2-D triangulation from polyshape 
 
  Methods for transformation:
    rotate          - Rotate a polyshape by angle with respect to a center
    scale           - Scale a polyshape by a factor 
    translate       - Translate a polyshape by a vector 
 
   Methods for Boolean operations:
    intersect	    - Find the intersection of two polyshapes or between polyshape and line 
    subtract        - Find the difference of two polyshapes 
    union           - Find the union of two polyshapes 
    xor             - Find the exclusive or of two polyshapes 
 
   Other methods: 
    polybuffer      - Create buffer zone around polyshape
    isequal         - Determine if two polyshapes are identical
    plot            - Plot polyshape and return a Polygon object handle 
    sortboundaries  - Sort boundaries in polyshape
    sortregions     - Sort regions in polyshape
    turningdist     - Find the turning distance of two polyshapes
 
  Example: Find the area and centroid of a polyshape and plot it
   %create a polyshape with 4 vertices
   quadShp = polyshape([0 0 1 3], [0 3 3 0]);
   %compute area and centroid
   a = area(quadShp);
   [cx, cy] = centroid(quadShp);
   figure; plot(quadShp);
   hold on
   %plot centroid point as '*'
   plot(cx, cy, '*');
 
  See also area, centroid, nsidedpoly

    Documentation for polyshape
       doc polyshape

help polyshape/area

  AREA Find the area of a polyshape
 
  A = AREA(pshape) returns the area of a polyshape. The area is the sum
  of the area of each solid region of pshape. 
 
  A = AREA(pshape, I) returns the area of the I-th boundary of pshape. This 
  syntax is only supported when pshape is a scalar polyshape object.
 
  See also perimeter, centroid, sortboundaries, polyshape

    Documentation for polyshape/area
       doc polyshape/area

The result will be a piecewise linear approximation to the area. That is the same result you will get from trapz anyway, since trapz is simply a tool that performs integration using linear segments to connect the points.

[Jacqueline Rigatto]

Thank you for the help, John D'Errico

Answer 3

[Ngo Nguyen]

You can use determinants to calculate the area of any polygons.

https://math.stackexchange.com/questions/264736/find-area-of-pentagon-using-determinants

For example:

function S = areaGate(P)

S = 0;

for i = 1:length(P)-1

S = S + det(P(i:i+1,:));

end

S = abs(S)/2;