How to calculate area between curve and horizontal line?

14 views (last 30 days)
Apologies if this has already been answered in another form elsewhere on the forums, however I couldn't find any answers that helped me directly with my code. I currently have code that maps the profile and area of a channel. I wish to have a horizontal y value across my channel width to simulate the fluctuating water level, the example line is shown in the image as the green line. So for example I would want to calculate and shade the area underneath the greenline only.
With the function would I be able to loop the same function for a range of values for the horizontal values, to provide an output I could then plot? A vector form would be ideal if possible.
I've already tried to use the area and fill commands unsuccesfully, or maybe I was not using the correct format for them.
I've also attached the relevant code below, any help would be much appreciated.
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = -[0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0]; % y values in image were produced by rand function
subplot(2, 1, 1);
area(x,-y)
grid on;
hold on;
for k = 1 : length(x)
xline(x(k), 'Color', 'k');
end
yt = yticks;
for k = 1 : length(yt)
yline(yt(k), 'Color', 'k');
end
plot (x, -y, 'r.-', 'LineWidth', 3, 'MarkerSize', 30)
hold on
yline(-0.5,'g',11);
xlabel('Chainage (m)', 'FontSize', 11);
ylabel('Depth (m)', 'FontSize', 11);
An = (trapz(x,y)); % Overall area.
% Compute areas within each pair of x locations.
midpoints = (y(1:end-1) + y(2:end))/2;
deltax = diff(x);
Ani = deltax .* midpoints;
subplot(2, 1, 2);
bar(x(1:end-1)-x(1)+ 0.5, Ani);
title('Areas of Slices', 'FontSize', 14);
xlabel('x', 'FontSize', 11);
ylabel('Area (m^2)', 'FontSize', 11);
grid on;

Accepted Answer

esat gulhan
esat gulhan on 21 Aug 2020
clc;clear;
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = [0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0];
h0= -0.5; %the axis in the graph
h=ones(size(x))*h0;
Int=pchip(x,y);
xx=linspace(x(1),x(end),500);
yy=y-h0;
yyy=pchip(x,yy,xx);
gt=yyy<0;
theArea =-trapz(xx(gt), (yyy(gt)))
hold off
area(xx, min(ppval(Int,xx), h(1)), h(1), 'EdgeColor', 'none', 'FaceColor', 'b'),grid,
hold on
plot(xx,ppval(Int,xx),'k',xx,h0,'Linewidth',2)
set(gca, 'XLim', [x(1) x(end)], 'YLim', [min(y) 0]);
title(strcat('BlueArea=', num2str(theArea)))
xlabel('Chainage (m)', 'FontSize', 11);
ylabel('Depth (m)', 'FontSize', 11);
You can change h0 which is horizontal line. when h0=-1
when h0=-0.5
when h0=0
I RESHAPED YOUR X AND Y DATAS 500 INTERVALS WITH PCHIP INTERPOLATION, IF NOT THE DATA IS NOT ENOUGH FOR GOOD SOLUTION.
  4 Comments
Hamzah Mahmood
Hamzah Mahmood on 22 Aug 2020
I'll look into that, cheers.
Another thing has come up, how would I be able to then keep these values as an array, im trying this however it isn't working.
clc;
clear;
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = [0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0];
for stage= [-2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0] %the axis in the graph
stage_dim= length(h)
h=ones(size(x))*stage;
Int=pchip(x,y);
xx=linspace(x(1),x(end),500);
yy=y-stage;
yyy=pchip(x,yy,xx);
gt=yyy<0;
Area =zeros(stage_dim);
for i = 1:stage_dim
Area(i) =-trapz(xx(gt), (yyy(gt)))
end
hold off
area(xx, min(ppval(Int,xx), h(1)), h(1), 'EdgeColor', 'none', 'FaceColor', 'b'),grid,
hold on
plot(xx,ppval(Int,xx),'k',xx,stage,'Linewidth',2)
set(gca, 'XLim', [x(1) x(end)], 'YLim', [min(y) 0]);
title(strcat('Area=', num2str(Area)))
xlabel('Chainage (m)', 'FontSize', 11);
ylabel('Depth (m)', 'FontSize', 11);
end
esat gulhan
esat gulhan on 23 Aug 2020
Sorry, i dont understand your question. If you ask this question in new stage, other users can see and answer...

Sign in to comment.

More Answers (2)

esat gulhan
esat gulhan on 21 Aug 2020
Edited: esat gulhan on 21 Aug 2020
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12]
y = -[0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0]
Int=pchip(x,y)
plot(x,ppval(Int,x))
MIItemp = (fnint(Int))
Area = fnval(MIItemp,12)%Area between 0 12
If you want to find area between 0-5 , Area = fnval(MIItemp,5)
There are other ways but this is the easiest way .
Note:If your matlab is latetest version you can write makima(x,y),instead of phcip(x,y). If your lines are strict use makima
  1 Comment
Hamzah Mahmood
Hamzah Mahmood on 21 Aug 2020
Hi Esat,
I think you may have misinterpered which axes I'm looking to cacluate the area from. The -0.5 value is a line that shifts up and down the y-axis rather than along the x-axis like so (the shade blue area).
Would your code allow me to calculate the area n the verticle directon if moodified slightly( I tried for a box with y = ones(1,13) however this didnt work). I'm relatively new to Matlab, so your code has helped introduce me to some new commands I hadn't come across before too.
Thanks.

Sign in to comment.


esat gulhan
esat gulhan on 22 Aug 2020
Did you see my second answer. I think it helps.
If you add pause to the code it will seem like an animation.
clc;clear;
x = [ 0 1 2 3 4 5 6 7 8 9 10 11 12];
y = [0 -0.5 -0.8 -0.8 -1 -1.1 -1.2 -1.2 -1.4 -1.2 -1.1 -1 0];
for h0= 0:-0.04:-2; %the axis in the graph
h=ones(size(x))*h0;Int=pchip(x,y);xx=linspace(x(1),x(end),500);
yy=y-h0;yyy=pchip(x,yy,xx);gt=yyy<0;theArea =-trapz(xx(gt), (yyy(gt)))
hold off
area(xx, min(ppval(Int,xx), h(1)), h(1), 'EdgeColor', 'none', 'FaceColor', 'b'),grid,
hold on
plot(xx,ppval(Int,xx),'k',xx,h0,'Linewidth',2)
set(gca, 'XLim', [x(1) x(end)], 'YLim', [min(y) 0]);
title(strcat('BlueArea=', num2str(theArea)))
xlabel('Chainage (m)', 'FontSize', 11);ylabel('Depth (m)', 'FontSize', 11);pause(0.00001)
end

Categories

Find more on Graph and Network Algorithms in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!