How to calculate partial area-under-the-curve?
14 次查看(过去 30 天)
显示 更早的评论
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
2021-3-15
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.
27 个评论
Jacqueline Rigatto
2021-3-15
编辑:Jacqueline Rigatto
2021-3-15
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
2021-3-15
Running the posted code:
Unrecognized function or variable 'MHWS_18'.
Waiting for it to magickally appear!
.
Star Strider
2021-3-15
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
2021-3-15
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
2021-3-15
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
2021-3-17
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
2021-3-18
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
2021-3-18
As always, my pleasure!
Not a bother at all! It’s definitely an interesting set of problems!
Jacqueline Rigatto
2021-3-18
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
2021-3-18
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
2021-3-19
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
2021-3-21
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
2021-3-21
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
2021-4-4
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
2021-4-4
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.
更多回答(2 个)
John D'Errico
2021-4-4
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.
1 个评论
Ngo Nguyen
2021-4-19
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;
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)