I'm trying to subtract one 3d surface plot from another to see see the change in the surface over time to track the level of Coastal erosion.

I'm making 3d models of a coast to track how much erosion is occuring on the coast over a period of time. I have collected data as northing easting and elevation and created 3 large arrays make a XYZ meshgrid. I am currently trying to subtract the final meshgrid from the initial but I cant get the correct change in elevation over time. The problem is the X and Y points dont match up exactly from one date to the next because the beach eroded back on the later data collection. Does anyone know a way I can fix the problem. This is my current code. The issue is with the final component labeled DOD (Difference of Difference of elevation models) for those interesed in coastal geology.

%DEM For Jupiter 05/18/21
x = [0	0.021	0.029	2.002	2.16	3.164	4.645	7.634	11.34	19.28	27.974	34.698	41.072	44.539	48.276	48.511	52.338	56.36	58.095	59.158	59.629	60.052	60.562	60.832	154.833	154.804	154.885	156.71	161.233	167.17	174.413	184.855	193.076	201.064	209.847	218.458	225.052	232.637	235.463	236.729	237.11	200.936	200.95	201.035	203.43	204.251	206.036	213.287	220.394	227.903	235.825	245.648	254.408	264.258	272.447	279.291	282.858	285.812	287.864	290.292	185.418	185.704	185.483	187.332	190.008	191.917	193.838	198.203	204.738	207.013	210.04	210.123	214.009	220.213	225.194	230.053	234.398	240.829	249.927	257.003	264.524	270.596	273.821	274.797	238.511	238.503	238.526	239.753	241.237	241.463	242.78	248.2	256.011	259.856	262.577	265.295	267.475	271.344	275.423	279.943	283.629	289.487	295.049	301.7	306.354	310.466	313.332	316.209	319.345	322.261	326.112];
y = [672.984	672.991	673.007	673.44	673.513	673.863	674.346	675.181	676.158	678.391	681.008	683.609	685.149	686.426	687.669	687.873	689.593	690.552	691.229	691.527	691.806	692.015	692.046	692.145	297.164	297.197	297.208	297.568	299.366	301.354	304.555	309.503	313.432	317.373	321.26	324.89	328.278	332.202	333.744	334.056	334.305	154.286	154.288	154.359	154.977	155.166	155.692	157.981	160.362	163.256	165.957	169.924	173.064	175.379	177.683	178.985	179.571	179.974	180.165	180.783	211.647	211.373	211.679	212.287	212.899	213.653	214.261	215.925	218.301	219.267	220.472	220.539	221.857	224.182	225.909	227.3	228.792	230.59	233.151	235.443	237.732	239.515	240.679	240.874	0	0.029	0.075	0.117	0.327	0.552	0.726	2.056	3.809	4.879	5.726	6.48	6.978	7.786	8.295	8.907	9.647	10.475	10.783	10.743	10.693	11.095	11.377	11.639	12.007	12.814	13.615];
z = [4.31	4.285	4.297	4.064	3.98	3.62	3.35	2.066	1.741	1.743	1.778	1.82	1.478	1.171	0.801	0.678	0.66	0.555	0.64	0.487	0.607	0.709	0.649	0.352	2.849	2.832	2.81	2.833	2.63	2.459	2.301	2.133	2.094	1.806	1.999	1.816	1.633	1.535	1.502	1.541	0.573	3.435	3.457	3.293	2.945	2.894	2.562	1.95	2.195	2.1	1.67	1.607	1.628	1.602	1.49	1.353	1.194	1.129	0.648	0.083	2.929	5.266	2.871	2.677	2.24	2.06	1.792	1.924	2.218	2.218	2.04	2.023	1.776	1.814	1.757	1.69	1.806	1.772	1.664	1.646	1.686	1.647	1.571	0.508	3.376	3.356	3.376	3.388	3.164	2.564	2.314	1.874	2.031	2.197	2.173	1.019	-0.033	0.054	-0.007	-0.114	-0.021	0.103	0.171	0.202	0.34	0.48	0.277	0.373	0.329	0.232	-0.143];
figure(1)
stem3(x, y, z)
grid on
xv = linspace(min(x), max(x), 20);
yv = linspace(min(y), max(y), 20);
[X,Y] = meshgrid(xv, yv);
Z = griddata(x,y,z,X,Y);
figure(2)
surf(X, Y, Z);
grid on
set(gca, 'ZLim',[-1 15])
shading interp
ylabel('Cross Shore(m)')
xlabel('Longshore(m)')
title('DEM Jupiter Inlet South 5/18/21')
colorbar
colormap(flipud(jet))

%DEM For Jupiter 07/19/21
x2 = [0	0.077	1.706	4.06	4.616	8.058	15.114	23.843	30.606	38.668	43.17	45.993	48.502	56.899	63.84	68.794	73.953	79.599	83.535	90.726	95.493	99.077	100.832	154.756	154.81	157.89	159.827	161.288	162.737	163.887	164.847	169.442	178.022	187.066	196.365	201.609	209.007	213.425	218.9	221.104	223.74	226.474	230.323	233.46	200.984	201.098	203.613	205.924	206.682	208.653	210.062	221.259	227.132	230.58	241.219	251.34	260.693	267.388	271.533	273.74	276.167	280.078	283.24	285.064	238.711	238.745	241.125	241.602	243.658	251.758	258.072	263.226	264.287	267.333	268.922	274.877	279.748	281.677	285.916	294.186	300.944	303.29	304.192	306.415	309.747	312.644	315.125];
y2 = [672.961	672.955	673.522	674.527	674.681	676.032	677.945	680.266	682.328	684.825	686.141	686.974	687.576	689.307	690.414	691.091	691.742	692.597	693.105	693.778	694.177	694.668	694.84	297.157	297.161	297.234	297.596	297.944	298.303	298.646	298.928	300.379	303.184	306.121	309.145	310.847	313.307	314.838	317.588	317.563	318.222	318.823	319.909	320.732	154.172	154.117	154.911	155.518	155.714	156.147	156.505	158.742	160.11	160.711	162.378	163.633	165.162	166.63	167.773	168.203	168.881	169.811	170.512	170.808	0.01	0	0.258	0.24	0.663	1.83	2.88	3.635	3.868	4.542	4.701	5.505	6.221	6.516	7.01	8.19	8.617	8.822	8.983	9.214	9.496	9.721	9.872];
z2 = [4.308	4.262	4.117	3.404	3.288	2.046	1.673	1.728	1.743	1.63	1.376	1	0.706	0.86	1.073	1.09	0.755	0.341	0.185	0.119	-0.237	-0.714	-0.95	2.81	2.756	2.77	2.697	2.786	2.715	2.536	2.446	2.482	2.372	2.215	2.077	1.956	1.748	1.691	1.754	0.325	0.098	-0.215	-0.746	-1.122	3.214	3.267	3.063	2.699	2.356	2.195	2.025	2.041	2.054	1.77	1.56	1.55	1.555	1.436	1.357	0.136	-0.161	-0.602	-0.974	-1.201	3.345	3.296	3.138	2.517	2.124	1.794	2.022	2.139	1.784	0.253	-0.126	-0.139	-0.095	0.145	0.333	0.631	0.747	0.755	0.164	-0.165	-0.655	-1.064	-1.401];
figure(3)
stem3(x2, y2, z2)
grid on
xv2 = linspace(min(x2), max(x2), 20);
yv2 = linspace(min(y2), max(y2), 20);
[X2,Y2] = meshgrid(xv2, yv2);
Z2 = griddata(x2,y2,z2,X2,Y2);
figure(4)
surf(X2, Y2, Z2);
grid on
set(gca, 'ZLim',[-1 15])
shading interp
ylabel('Cross Shore(m)')
xlabel('Longshore(m)')
title('DEM South Inlet Park 7/19/21')
colorbar
colormap(flipud(jet))

%DOD
Z4 = Z2-Z;
figure(5)
surf(X, Y, Z4);
grid on
set(gca, 'ZLim',[-1 15])
shading interp
ylabel('Cross Shore(m)')
xlabel('Longshore(m)')
title('Jupiter DOD Last-First')
colorbar
colormap(flipud(jet))
caxis([-2.5 1])

Answer 1

[Roshan Mathew Tom]

I would suggest you to use the interpolation function to interpolate data onto a fixed set of reference points. The code would go as follows:

[x1, y1, z1] = DataFromDay1;
[x2, y2, z2] = DataFromDay2;

z2New = interp2(x2,y2,z2,x1,y1);

Diff = z2New - z1;

surf(x1, y1, Diff);

[Austin Scheinkman]

Sorry if this is a dumb question but im pretty new to matlab. How exactly do I put my data in DataFromDay1 and DataFromDay2 in your code. Like what format?

[Roshan Mathew Tom]

In your code, that would be

DataFromDay1 = [x;y;z]';
DataFromDay2 = [x2;y2;z2]';