How can i calculate the area under two curves that intersect?

I have two curves that plot by x-y coordinates

How can I calculate this area?

[Dyuman Joshi]

Do you want to area for the 2 highlighted areas or all intersecting areas?

[Sanufee]

Yess

[Sanufee]

That's a two diffent cross section of the beach, So I want to know the volume that transiton

Answer 1

[Star Strider]

Integrate them using trapz then do what you want with the results (keep them as they are, add them, add their absolute values, etc.) —

x = linspace(0, 40, 120);
y1 = sin(2*pi*x/30);
y2 = cos(2*pi*x/30);

L = numel(x);
idx = find(diff(sign(y1-y2)));                                              % indices Of Intersections
% q = x(idx)
for k = 1:numel(idx)-1                                                      % Integrate Between Indices
    idxrng = max(1,idx(k)) : min(L,idx(k+1));                               % Index Range
    A(k) = trapz(x(idxrng), y1(idxrng)) - trapz(x(idxrng), y2(idxrng));
end

A

A = 1×2

13.4886 -13.4862

for k = 1:numel(idx)-1                                                      % Interpolate Then Integrate
    idxrng1 = max(1,idx(k)-1) : min(L,idx(k)+1);                            % Index Range (First Inmtersection)
    idxrng2 = max(1,idx(k+1)-1) : min(L,idx(k+1)+1);                        % Index Range (Second Inmtersection)
    xi(1) = interp1(y1(idxrng1)-y2(idxrng1), x(idxrng1), 0);                % Find First 'x' Intersection
    xi(2) = interp1(y1(idxrng2)-y2(idxrng2), x(idxrng2), 0);                % Find Second 'x' Intersection
    xv = x((x>=xi(1)) & (x<=xi(2)));                                        % 'x' Vector Between Them
    y1v = interp1(x, y1, xv);                                               % Corresponding 'y1' Vector
    y2v = interp1(x, y2, xv);                                               % Corresponding 'y2' Vector
    A(k) = trapz(xv, y1v) - trapz(xv, y2v);                                 % Integrate
end

A

A = 1×2

13.4771 -13.4956

figure
plot(x, y1)
hold on
plot(x, y2)
hold off
grid

Here are two options. I suggest using the second one (interpolate first), since that seems to match your data more closely.

EDIT — (21 Feb 2024 02:35)

Guessing the data —

x1 = [2 4 7 9 10 12 14 17 19 22 24 27 29 32 34 37];
y1 = [-0.25 -0.4 -0.45 -0.6  -1.0  -1.2  -1.4  -1.45 -1.45 -1.5  -1.51  -1.6  -1.9  -2.3 -2.7 -2.9];
x2 = [2 4 7 9 10 12 13 17 19 22 23 27 29 34 36];
y2 = [-0.26 -0.4 -0.50 -0.7  -0.8  -1.1  -1.5  -1.7 -1.7  -1.7 -1.7  -1.8  -2.1 -2.6 -2.8];
% x1y1 = [x1; y1]
% x2y2 = [x2; y2]

xv = linspace(min([x1 x2]), max([x2 x2]), 100);                             % Common 'x' Vector
y1v = interp1(x1, y1, xv, 'pchip');                                         % Interpolate Using 'pchip'
y2v = interp1(x2, y2, xv, 'pchip');                                         % Interpolate Using 'pchip'


x = xv;
y1 = y1v;
y2 = y2v;

L = numel(xv);
idx = find(diff(sign(y1-y2)));                                              % indices Of Intersections

for k = 1:numel(idx)-1                                                      % Interpolate Then Integrate
    idxrng1 = max(1,idx(k)-1) : min(L,idx(k)+1);                            % Index Range (First Inmtersection)
    idxrng2 = max(1,idx(k+1)-1) : min(L,idx(k+1)+1);                        % Index Range (Second Inmtersection)
    xi(1) = interp1(y1(idxrng1)-y2(idxrng1), x(idxrng1), 0);                % Find First 'x' Intersection
    xi(2) = interp1(y1(idxrng2)-y2(idxrng2), x(idxrng2), 0);                % Find Second 'x' Intersection
    xv = x((x>=xi(1)) & (x<=xi(2)));                                        % 'x' Vector Between Them
    y1v = interp1(x, y1, xv);                                               % Corresponding 'y1' Vector
    y2v = interp1(x, y2, xv);                                               % Corresponding 'y2' Vector
    A(k) = trapz(xv, y1v) - trapz(xv, y2v);                                 % Integrate
end

figure
plot(x, y1, '-b', 'MarkerFaceColor','b')
hold on
plot(x, y2, '-r', 'MarkerFaceColor','r')
hold off
grid
text(10.5, -1, sprintf('\\leftarrow Area = %.3f', A(3)))
text(26, -1.7, sprintf('\\leftarrow Area = %.3f', A(4)))
text(20, -0.5, sprintf('Total Area = %.3f', sum(abs(A([3 4])))))

.

[Sanufee]

[Sanufee]

This happended to me 😢

[Star Strider]

It would help to have your data.

My variables are all row vectors, so my guess is that yours are column vectors.

In that instance, one option would be:

xv = linspace(min([x1; x2]), max([x2; x2]), 100); % Common 'x' Vector

and another:

xv = linspace(min([x1 x2],[],'all'), max([x2 x2],[],'all'), 100); % Common 'x' Vector

At least one of those should work.

If you are still having problems with your data and my code, please provide your data as a text, spreadhseet, or .mat file.

Converting my initial vectors to column vectors works with the first option, and no other changes in the rest of my code —

x1 = [2 4 7 9 10 12 14 17 19 22 24 27 29 32 34 37].';
y1 = [-0.25 -0.4 -0.45 -0.6  -1.0  -1.2  -1.4  -1.45 -1.45 -1.5  -1.51  -1.6  -1.9  -2.3 -2.7 -2.9].';
x2 = [2 4 7 9 10 12 13 17 19 22 23 27 29 34 36].';
y2 = [-0.26 -0.4 -0.50 -0.7  -0.8  -1.1  -1.5  -1.7 -1.7  -1.7 -1.7  -1.8  -2.1 -2.6 -2.8].';
% x1y1 = [x1; y1]
% x2y2 = [x2; y2]

xv = linspace(min([x1; x2]), max([x2; x2]), 100);                           % Common 'x' Vector
y1v = interp1(x1, y1, xv, 'pchip');                                         % Interpolate Using 'pchip'
y2v = interp1(x2, y2, xv, 'pchip');                                         % Interpolate Using 'pchip'


x = xv;
y1 = y1v;
y2 = y2v;

L = numel(xv);
idx = find(diff(sign(y1-y2)));                                              % indices Of Intersections

for k = 1:numel(idx)-1                                                      % Interpolate Then Integrate
    idxrng1 = max(1,idx(k)-1) : min(L,idx(k)+1);                            % Index Range (First Inmtersection)
    idxrng2 = max(1,idx(k+1)-1) : min(L,idx(k+1)+1);                        % Index Range (Second Inmtersection)
    xi(1) = interp1(y1(idxrng1)-y2(idxrng1), x(idxrng1), 0);                % Find First 'x' Intersection
    xi(2) = interp1(y1(idxrng2)-y2(idxrng2), x(idxrng2), 0);                % Find Second 'x' Intersection
    xv = x((x>=xi(1)) & (x<=xi(2)));                                        % 'x' Vector Between Them
    y1v = interp1(x, y1, xv);                                               % Corresponding 'y1' Vector
    y2v = interp1(x, y2, xv);                                               % Corresponding 'y2' Vector
    A(k) = trapz(xv, y1v) - trapz(xv, y2v);                                 % Integrate
end

figure
plot(x, y1, '-b', 'MarkerFaceColor','b')
hold on
plot(x, y2, '-r', 'MarkerFaceColor','r')
hold off
grid
text(10.5, -1, sprintf('\\leftarrow Area = %.3f', A(3)))
text(26, -1.7, sprintf('\\leftarrow Area = %.3f', A(4)))
text(20, -0.5, sprintf('Total Area = %.3f', sum(abs(A([3 4])))))

.

[Sanufee]

[Sanufee]

And then I try to put the data by manual. This happended, the area calculated only 1 show

[Sanufee]

[Sanufee]

x1 = [0.8961 3.3569 5.8782 8.3384 10.8515 13.3112 15.8254 18.3049 20.8188 23.3092 25.771 28.2613 29.7442 32.2171 33.1798 35.9226].';

y1 = [-0.3838 -0.5309 -0.762 -0.9384 -1.1361 -1.2114 -1.4222 -1.2912 -1.3113 -1.5192 -1.8071 -2.0278 -2.1672 -2.4575 -2.5841 -2.9109].';

x2 = [2.0252 4.4448 6.8848 9.4151 12.0184 14.3916 17.01 19.3367 21.8709 24.3984 26.8193 29.3898 31.8715 34.2408 35.8614].';

y2 = [-0.3586 -0.5645 -0.75 -1.0408 -1.0724 -1.1438 -1.164 -1.5574 -1.8426 -1.7756 -1.7499 -1.9891 -2.2325 -2.3882 -2.5566].';

[Sanufee]

Can you text for this data pls/

[Sanufee]

This is my spreadsheet hope this will able to do 🥺

[Sanufee]

I am a beginner for matlab my professor tell me to try to solve it

[Star Strider]

Thank you!

If you post your data, I can tweak my code to produce the areas you want. The ‘A’ vector will have areas for each region my code identifies (between crossings). (In the last image, there should be about 7.)

Are there specific criteria as to which areas you particularly want to identify?

[Sanufee]

Could you help to give me a +value on green area (Its mean in Post-monsoon the sand has been increasing) and -value on orange area (Its mean in Post-monsoon the sand has been decreasing.) Thankks.

[Star Strider]

I needed to tweak my code a bit to include the areas at both ends (defined by only one crossing). I could not understand what information you want from the spreadsheet files, so I used the supplied vectors.

Try this —

% Pre = readtable('Premonsoon.xlsx')
% Post = readtable('Postmonsoon.xlsx')

x1 = [0.8961 3.3569 5.8782 8.3384 10.8515 13.3112 15.8254 18.3049 20.8188 23.3092 25.771 28.2613 29.7442 32.2171 33.1798 35.9226].';
y1 = [-0.3838 -0.5309 -0.762 -0.9384 -1.1361 -1.2114 -1.4222 -1.2912 -1.3113 -1.5192 -1.8071 -2.0278 -2.1672 -2.4575 -2.5841 -2.9109].';
x2 = [2.0252 4.4448 6.8848 9.4151 12.0184 14.3916 17.01 19.3367 21.8709 24.3984 26.8193 29.3898 31.8715 34.2408 35.8614].';
y2 = [-0.3586 -0.5645 -0.75 -1.0408 -1.0724 -1.1438 -1.164 -1.5574 -1.8426 -1.7756 -1.7499 -1.9891 -2.2325 -2.3882 -2.5566].';

xv = linspace(min([x1; x2]), max([x2; x2]), 100);                           % Common 'x' Vector
y1v = interp1(x1, y1, xv, 'pchip');                                         % Interpolate Using 'pchip'
y2v = interp1(x2, y2, xv, 'pchip');                                         % Interpolate Using 'pchip'
x1a = x1;
y1a = y1;
x2a = x2;
y2a = y2;

x = xv;
y1 = y1v;
y2 = y2v;

L = numel(xv);
idx = find(diff(sign(y1-y2)));                                              % indices Of Intersections
idx = [1 idx L];                                                            % Include Both Ends

for k = 1:numel(idx)-1                                                      % Interpolate Then Integrate
    idxrng1 = max(1,idx(k)-1) : min(L,idx(k)+1);                            % Index Range (First Inmtersection)
    idxrng2 = max(1,idx(k+1)-1) : min(L,idx(k+1)+1);                        % Index Range (Second Inmtersection)
    xi(1) = interp1(y1(idxrng1)-y2(idxrng1), x(idxrng1), 0);                % Find First 'x' Intersection
    xi(2) = interp1(y1(idxrng2)-y2(idxrng2), x(idxrng2), 0);                % Find Second 'x' Intersection
    if isnan(xi(1))                                                         % The 'NaN' Values Result Of There Are No Intersections (At Each End), So Test For Them & Make Appropriate Substitutions
        xi(1) = x(1);
    elseif isnan(xi(2))
        xi(2) = x(end);
    end
    xv = x((x>=xi(1)) & (x<=xi(2)));                                        % 'x' Vector Between Them
    y1v = interp1(x, y1, xv);                                               % Corresponding 'y1' Vector
    y2v = interp1(x, y2, xv);                                               % Corresponding 'y2' Vector
    A(k) = trapz(xv, y1v) - trapz(xv, y2v);                                 % Integrate
    xc(k) = mean(xv);                                                       % Used For 'text' Loations
    yc(k) = mean([y1v y2v]);                                                % Used For 'text' Loations
end

A                                                                           % Segment Areas

A = 1×5

-0.5386 0.0191 -1.1436 2.0686 -2.0872

figure
hp(1) = plot(x1a, y1a, '.b', 'MarkerFaceColor','b', 'DisplayName','Premonsoon');
hold on
plot(x, y1,  '-b')
hp(2) = plot(x2a, y2a, '.r', 'MarkerFaceColor','r', 'DisplayName','Postmonsoon');
plot(x, y2, '-r')
hold off
grid
axis([0 45  -3.5 0])
legend(hp,'Location','best')
for k = 1:numel(A)
    text(xc(k), yc(k), sprintf('\\leftarrow Area = %.3f', A(k)))
end
text(7, -2.25, sprintf('Increasing  = %.3f\nDecreasing = %.3f',abs(sum(A([1:3 5]))),abs(A(4))))

I did not use patch to colour the individual areas, although I can if that is important.

[Sam Chak]

Hi @Sanufee,

If you plan to present the patchy areas in a technical report, it would be helpful to include a detailed explanation of the underlying numerical method used to estimate these areas. This will enhance the clarity and rigor of your report.

[Sanufee]

@Star Strider Omg! That's work. Thank you so much. 🥹🥹🥹

If If it's not disturbing. Can you help me to code for ploting the surface (Beach surface) with 3 Coordinates (x,y,z) data

here's example

I don't know much about interpolate func in MatLAB.

If you can solve that please let me see 🥲

[Sanufee]

@Sam Chak Thanks for the advice sir! 🫡

[Star Strider]

@Sanufee — I may be able to. I need a file with the necessary information. (It looks like a surf plot.)

[Sanufee]

This surface plot by (E,N) coordinate in my file

and that have Elevation (It's should be z axes value in 3D surface) of each point then i want combine it to 3D surface sir.

[Star Strider]

It is best to use the scatteredInterpolant function to calcualte the surface from the data. (I did it two ways here with the same scatteredInterpolant function, first with the original unsorted vectors to do the interpolation to calculate the surface, and second with sorted vectors. The interpolation method choices are 'linear', 'nearest', and 'natural'. I experimented with them, and like 'natural' best in this instance.) Choose the appropriate colormap and view to depict it as you want. I also used the surfc function to draw the surfaces, since the contour plot lets me see details that I would otherwise not be able to. Use it or surf to drtaw the surfaces.

Post = readtable('Postmonsoon.xlsx')

Post = 529×4 table
    Point      E          N       Elevation
    _____    ______    _______    _________

     NaN     2.0642    -0.0163     -0.1885 
     NaN          0          0           0 
       1     4.4896    -0.0456     -0.3672 
       2     4.4866    -0.0455     -0.3671 
       3     6.9964    -0.0443     -0.4542 
       4     9.4778    -0.0158     -0.7015 
       5     12.005    -0.0133     -0.8593 
       6     14.603    -0.0125     -1.0918 
       7     16.983    -0.0117     -1.4529 
       8       19.3     -0.011     -1.6625 
       9     21.899    -0.0101     -1.7106 
      10      23.91    -0.0145     -1.6612 
      11      26.93    -0.0143     -1.8578 
      12     29.437     -0.014      -2.123 
      13     34.331    -0.0138      -2.564 
      14       36.5     -0.049     -2.8317 

E = Post.E(3:end);
N = Post.N(3:end);
Elev = Post.Elevation(3:end);



F = scatteredInterpolant(E,N,Elev,'natural')

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

F = 
  scatteredInterpolant with properties:

                 Points: [526×2 double]
                 Values: [526×1 double]
                 Method: 'natural'
    ExtrapolationMethod: 'linear'

[Em,Nm] = ndgrid(E, N);
Elevm = F(Em,Nm);

figure
surfc(Em, Nm, Elevm, 'EdgeColor','none')
colormap(turbo)
xlabel('E')
ylabel('N')
zlabel('Elevation')
title('Unsorted Vectors')
view(135, 30)

[Em,Nm] = ndgrid(sort(E), sort(N));
Elevm = F(Em,Nm);

figure
surfc(Em, Nm, Elevm, 'EdgeColor','interp')
colormap(turbo)
xlabel('E')
ylabel('N')
zlabel('Elevation')
title('Sorted Vectors')
view(135, 30)

.

[Sanufee]

That's workkkk! It's look so good than my expect. Can you tell me how to add the Z=-2 plane in to it. I have tried it but it split 2 window of figure. It's not ploting on the same figure😢

Post = readtable('Postmonsoon.xlsx')

E = Post.E(3:end);
N = Post.N(3:end);
Elev = Post.Elevation(3:end);



F = scatteredInterpolant(E,N,Elev,'natural')

[Em,Nm] = ndgrid(E, N);
Elevm = F(Em,Nm);

%figure
%surfc(Em, Nm, Elevm, 'EdgeColor','none')
%colormap(turbo)
%xlabel('E')
%ylabel('N')
%zlabel('Elevation')
%title('Unsorted Vectors')
%view(135, 30)



[Em,Nm] = ndgrid(sort(E), sort(N));
Elevm = F(Em,Nm);

figure
surfc(Em, Nm, Elevm, 'EdgeColor','interp')
colormap(turbo)
xlabel('E')
ylabel('N')
zlabel('Elevation')
title('Sorted Vectors')
view(135, 30)

% Adding the 0 0 0 plane (I asked chatgpt)
[X,Y] = meshgrid(min(Em(:)):max(Em(:)), min(Nm(:)):max(Nm(:)));
Z = zeros(size(X));
surf(X, Y, Z, 'FaceAlpha', 0.3, 'FaceColor', 'r', 'EdgeColor', 'none');

[Sanufee]

and could you help me to expand the E axes?

[Sanufee]

Ohh, sir I got it! I already known how to combine it. I use hold on and hold off. but I still don't know how to expand E axes.

[Star Strider]

The Z = -2 plane is straightforward. You can either use surf or patch. The'FaceAlpha' and 'EdgeAlpha' name-value pairs control the transparency of the surface. I set them both to 0.5 here, increasing that value makes them less transparent, and decreasing it makes them more transparent. The allowed range is from 0 to 1.

I am not certain what you want with respect to expanding the E axis. One option is to set the limits to be larger with

xlim([-10 50])

since the ‘x’ axis is the ‘E’ axis, however it would not be a good idea to expand the surface to match it. That would require extrapolation, and that is never appropriate because it is impossible to know what the surface (or any such function) actually does outside the region that originally defines it. Interpolating is appropriate, however extrapolating never is appropriate.

Post = readtable('Postmonsoon.xlsx')

Post = 529x4 table
    Point      E          N       Elevation
    _____    ______    _______    _________

     NaN     2.0642    -0.0163     -0.1885 
     NaN          0          0           0 
       1     4.4896    -0.0456     -0.3672 
       2     4.4866    -0.0455     -0.3671 
       3     6.9964    -0.0443     -0.4542 
       4     9.4778    -0.0158     -0.7015 
       5     12.005    -0.0133     -0.8593 
       6     14.603    -0.0125     -1.0918 
       7     16.983    -0.0117     -1.4529 
       8       19.3     -0.011     -1.6625 
       9     21.899    -0.0101     -1.7106 
      10      23.91    -0.0145     -1.6612 
      11      26.93    -0.0143     -1.8578 
      12     29.437     -0.014      -2.123 
      13     34.331    -0.0138      -2.564 
      14       36.5     -0.049     -2.8317 

E = Post.E(3:end);
N = Post.N(3:end);
Elev = Post.Elevation(3:end);

F = scatteredInterpolant(E,N,Elev,'natural');

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

[Em,Nm] = ndgrid(sort(E), sort(N));
Elevm = F(Em,Nm);

figure
surfc(Em, Nm, Elevm, 'EdgeColor','interp')
hold on
surf(Em, Nm, -2*ones(size(Em)), 'FaceAlpha',0.5, 'EdgeAlpha',0.5, 'EdgeColor','none')
hold off
colormap(turbo)
xlabel('E')
ylabel('N')
zlabel('Elevation')
title('Sorted Vectors')
view(135, 30)
xlim([-10 50])

.

[Sanufee]

Now, I want to plot the 2 surface in same fig for comparison. So i cant change different color of surface (Colormap). Can you help please?

I tried it by change in fiqure manualy but It's not color that i want. I want to change to different colormap to see the different clearly.

Then that's better if I can show the legend also.

Thankss

Here's my code

Pre = readtable('Premonsoon.xlsx')
Post = readtable('Postmonsoon.xlsx')
E1 = Pre.E(3:end);
N1 = Pre.N(3:end);
Elev1 = Pre.Elev(3:end);

E2 = Post.E(3:end);
N2 = Post.N(3:end);
Elev2 = Post.Elevation(3:end);


F1 = scatteredInterpolant(E1,N1,Elev1,'natural')

[Em1,Nm1] = ndgrid(sort(E1), sort(N1));
Elevm1 = F1(Em1,Nm1);

F2 = scatteredInterpolant(E2,N2,Elev2,'natural')

[Em2,Nm2] = ndgrid(sort(E2), sort(N2));
Elevm2 = F2(Em2,Nm2);

% Plot
hold on
% Pre-monsoon surface
surfc(Em1, Nm1, Elevm1, 'EdgeColor','interp')
colormap("turbo") % Change colormap for pre-monsoon surface
title('Pre-monsoon')
view(140, 30)
hold off

hold on
% Post-monsoon surface
surfc(Em2, Nm2, Elevm2, 'EdgeColor','interp')
colormap("parula") % Change colormap for post-monsoon surface
title('Post-monsoon')
view(140, 30)
hold off

% Legend
h = colorbar;
ylabel(h,'Elevation')


hold on
[X,Y] = meshgrid(min(Em1(:)):max(Em1(:)), min(Nm1(:)):max(Nm1(:)));
Z = -2 * ones(size(X)); % Set Z values to -2
surf(X, Y, Z, 'FaceAlpha', 0.3, 'FaceColor', '#4DBEEE', 'EdgeColor', 'none');
hold off

[Star Strider]

That is apparently not possible. The colormap has to apply to an axis, not a specific surf plot.

One option is to set each surf plot to a different colour, howefvr the best approach may be to interpolate the larger values (‘Premonsoon’) to the size of the smaller values (‘Postmonsoon’) so they are both the same size (this requires a third scatteredInterpolant call) and then subtract them, as I did at the end (using a slightly different approach to plotting it) —

Pre = readtable('Premonsoon.xlsx')

Pre = 606×4 table
    No      E         N        Elev  
    __    ______    ______    _______

     1    2.0028    0.0003     -0.247
     2    4.4752    0.1113    -0.3718
     3    6.9702    0.1733    -0.4969
     4     9.471    0.2356    -0.6752
     5    11.946    0.2971    -0.9908
     6     14.44     0.359    -1.1287
     7    16.893    0.4018    -1.3546
     8    19.425    0.4804    -1.4187
     9    21.903    0.5552    -1.4402
    10    24.411    0.6455    -1.5291
    11    26.882     0.673    -1.7164
    12    29.372    0.7274    -1.9758
    13    31.862    0.7887    -2.2727
    14    34.376    0.8114     -2.646
    15    36.864     0.836    -2.9736
    16    1.7236    9.9816    -0.2711

Post = readtable('Postmonsoon.xlsx')

Post = 531×6 table
    Var1     Var2      Var3       Var4              Var5            Var6
    ____    ______    _______    _______    ____________________    ____

    NaN          0          0          0    {'10 Feb, Saturday'}    NaN 
      1     4.4896    -0.0456    -0.3672    {0×0 char          }    NaN 
      2     4.4866    -0.0455    -0.3671    {0×0 char          }    NaN 
      3     6.9964    -0.0443    -0.4542    {0×0 char          }    NaN 
      4     9.4778    -0.0158    -0.7015    {0×0 char          }    NaN 
      5     12.005    -0.0133    -0.8593    {0×0 char          }    NaN 
      6     14.603    -0.0125    -1.0918    {0×0 char          }    NaN 
      7     16.983    -0.0117    -1.4529    {0×0 char          }    NaN 
      8       19.3     -0.011    -1.6625    {0×0 char          }    NaN 
      9     21.899    -0.0101    -1.7106    {0×0 char          }    NaN 
     10      23.91    -0.0145    -1.6612    {0×0 char          }    NaN 
     11      26.93    -0.0143    -1.8578    {0×0 char          }    NaN 
     12     29.437     -0.014     -2.123    {0×0 char          }    NaN 
     13     34.331    -0.0138     -2.564    {0×0 char          }    NaN 
     14       36.5     -0.049    -2.8317    {0×0 char          }    NaN 
     15     2.1108     9.9106    -0.2643    {0×0 char          }    NaN 

E1 = Pre.E(3:end);
N1 = Pre.N(3:end);
Elev1 = Pre.Elev(3:end);

% E2 = Post.E(3:end);
% N2 = Post.N(3:end);
% Elev2 = Post.Elevation(3:end);
E2 = Post{3:end,2};                                     % No Appropriate Variable Names In 'Post'
N2 = Post{3:end,3};
Elev2 = Post{3:end,4};


F1 = scatteredInterpolant(E1,N1,Elev1,'natural');

[Em1,Nm1] = ndgrid(sort(E1), sort(N1));
Elevm1 = F1(Em1,Nm1);

F2 = scatteredInterpolant(E2,N2,Elev2,'natural');

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

[Em2,Nm2] = ndgrid(sort(E2), sort(N2));
Elevm2 = F2(Em2,Nm2);

% Plot
hold on
% Pre-monsoon surface
surfc(Em1, Nm1, Elevm1, 'EdgeColor','interp', 'FaceColor','r')
% colormap(turbo) % Change colormap for pre-monsoon surface
title('Pre-monsoon')
view(140, 30)
hold off

hold on
% Post-monsoon surface
surfc(Em2, Nm2, Elevm2, 'EdgeColor','interp', 'FaceColor','g');
% colormap(parula) % Change colormap for post-monsoon surface
title('Post-monsoon')
view(140, 30)
hold off

% Legend
h = colorbar;
ylabel(h,'Elevation')


hold on
[X,Y] = meshgrid(min(Em1(:)):max(Em1(:)), min(Nm1(:)):max(Nm1(:)));
Z = -2 * ones(size(X)); % Set Z values to -2
surf(X, Y, Z, 'FaceAlpha', 0.3, 'FaceColor', '#4DBEEE', 'EdgeColor', 'none');
hold off

size(Elevm1)

ans = 1×2

604 604

size(Elevm2)

ans = 1×2

529 529

F3 = scatteredInterpolant(Em1(:),Nm1(:),Elevm1(:));

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

Elevm1i = F3(Em2,Nm2);                                      % 'Elev1m' Interpolated To The dimensions of 'Elev2m'

figure
hsl = surfl(Em2,Nm2,Elevm2-Elevm1i, 'light')

hsl = 
  1×2 graphics array:

    Surface    Light  

hsl(1).EdgeColor = 'none';
% surfc(Em2,Nm2,Elevm2-Elevm1i, 'EdgeColor','interp')
colormap(turbo)
Ax = gca;
Ax.YDir = 'reverse';
Ax.XDir = 'reverse';
title('Post-Pre Difference: Elevm2-Elevm1 (Uses Interpolated Elevm1)')
zlim([-2 1])

figure
contourf(Em2,Nm2,Elevm2-Elevm1i, (-2:0.25:1), 'ShowText','on')
Ax = gca;
colormap(turbo)
colorbar
Ax.YDir = 'reverse';
Ax.XDir = 'reverse';
title('Post-Pre Difference: Elevm2-Elevm1 (Uses Interpolated Elevm1)')

.

[Sanufee]

It's so useful for me! Thanks so much. maybe I could give my graduate certificate to you. 🥹

[Star Strider]

As always, my pleasure!

(I already have a few of my own! I am happy to be able to help you.)

Answer 2

[John D'Errico]

1. Subtract the two curves.

2. Take the absolute value of the difference.

3. Compute the integral.

Note that you may need the curves in a functional form if you want to do this accurately, so you might need an interpolation tool to do that. This is often in the form of a spline.

X = linspace(0,10,11);
Y1 = rand(1,11);
Y2 = rand(1,11);
f1 = spline(X,Y1);
f2 = spline(X,Y2);
fnplt(f1,'r')
hold on
fnplt(f2,'b')
plot(X,Y1,'ko',X,Y2,'kx')
hold off

AreaBetween = integral(@(x) abs(fnval(f1,x) - fnval(f2,x)),0,10)

AreaBetween = 3.2668

Note the improtance of using integral here, and of the use of functions in the form of splines to interpolate. The problem is, if you do not, then when the curves cross, you will get an error, and that error could be significant. For example, we could just use trapz directly on the data.

trapz(X,abs(Y1 - Y2))

ans = 3.5494

You should see there is a very significant difference.

plot(X,abs(Y1 - Y2),'o-')
hold on
fplot(@(x) abs(fnval(f1,x) - fnval(f2,x)),[0,10])
hold off

Where those curves actually crossed, the difference function should go all the way to zero. That is what the red curve does. In fact, the red curve is what you actually want to integrate.

But if the curves cross between two of the data points, the result will be the blue line. And that will totally mess up the trapezoidal integration, since it will see only the blue curve. You can see here, the difference is approximately a 10% error in the area estimate. The error from trapz will usually be an overestimate of the area.

[Sanufee]

Thanks!! I will try it