MATLAB: Using ‘fit’ to add a trend line to series of scatter points

fitlatlonplotployfitpolyvaltrend linexy

Hi there,
I am trying to add a trendline to a series of scatter points that represent geographic locations and have an estimate of the goodness of fit.
Here is a simplified version of my data- please can someone help me pin point why/ where it is going wrong. This is the output I get and as you can see this cannot be right from looking at the values on my y-axis
f =
Linear model Poly1:
f(x) = p1*x + p2
Coefficients (with 95% confidence bounds):
p1 = 0.1116 (-0.04946, 0.2726)
p2 = 37.37 (32.47, 42.27)
clear all
close all
lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000]
lon= [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000]
h=[30.3 30.45 40.72 40.8];
lat_lon_proportions(h) %function taken from Jonathen Sullivan - available from https://uk.mathworks.com/matlabcentral/fileexchange/32462-correctly-proportion-a-lat-lon-plot
% plot(lat,lon,'k.')
x=lon
y=lat
hold on
f=fit(y,x,'poly1')
plot(f,y,x)

Best Answer

That is the correct first-order least-squares fit. But I think in your plot the '.' marker on the far left side is coincident with the y axis. The result is it looks like the fit is off. Try using 'o' markers instead, as shown below.
Also, as a general practice to help with debugging these kinds of things, I recommend trying not to change variable names unless there's a really compelling reason to do so. Setting x=lon and y=lat just adds confusion. It's particularly confusing in this case because the syntax plot(f,y,x) plots y on the x axis and x on the y axis, but automatically labels the x data as y, and labels the y data as x. I recommend this rewrite:
lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000];
lon = [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000];
f=fit(lat,lon,'poly1')
plot(f,lat,lon,'o')
xlabel latitude
ylabel longitude
Or, if you want to keep longitudes on the x axis, as is common, this:
lat = [30.4056270000000;30.3991500000000;30.4000740000000;30.4046990000000;30.4102510000000;30.4093280000000;30.3889660000000;30.3982300000000;30.4074760000000;30.4195050000000;30.4139500000000;30.4120980000000;30.4083980000000;30.3936010000000;30.3926810000000];
lon = [40.7602680000000;40.7609660000000;40.7616710000000;40.7623780000000;40.7623830000000;40.7602710000000;40.7644750000000;40.7574480000000;40.7616770000000;40.7609820000000;40.7637930000000;40.7651980000000;40.7637880000000;40.7595540000000;40.7560360000000];
f=fit(lon,lat,'poly1')
plot(f,lon,lat,'o')
xlabel longitude
ylabel latitude