- You will see updates in your activity feed.
- You may receive emails, depending on your notification preferences.

388 views (last 30 days)

Show older comments

Jacqueline Rigatto
on 15 Mar 2021

Answered: Ngo Nguyen
on 19 Apr 2021 at 5:36

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

Star Strider
on 15 Mar 2021

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
on 15 Mar 2021

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
on 15 Mar 2021

Running the posted code:

Unrecognized function or variable 'MHWS_18'.

Waiting for it to magickally appear!

.

Star Strider
on 15 Mar 2021

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.

Star Strider
on 15 Mar 2021

My pleasure!

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

.

Jacqueline Rigatto
on 15 Mar 2021

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
on 15 Mar 2021

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
on 17 Mar 2021

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
on 18 Mar 2021

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

.

Star Strider
on 18 Mar 2021

As always, my pleasure!

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

Jacqueline Rigatto
on 18 Mar 2021

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
on 18 Mar 2021

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.

Star Strider
on 19 Mar 2021

Jacqueline Rigatto
on 21 Mar 2021

Star Strider
on 21 Mar 2021

Jacqueline Rigatto
on 4 Apr 2021

Star Strider
on 4 Apr 2021

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.

John D'Errico
on 4 Apr 2021

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.

Ngo Nguyen
on 19 Apr 2021 at 5:36

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

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;

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)