matrix calculation reshaping and diferences between datum

i dont know how to transform the Z dimension into a 2x2 matriz at least i dont understand how to make it run

[Vasco]

heres the file

[Voss]

You may want

delta_lat_Dlx = lat_dec - lat_3_DLx.';

or

delta_lat_Dlx = lat_dec.' - lat_3_DLx;

[Vasco]

my teacher wants it this way

[Voss]

Use contour instead of contourf, and use subplot or tiledlayout.

[Vasco]

the data appears all wrong the limits of the country dont appear

i dont know whats happening

[Voss]

Where would the country border come from?

I don't know what's going on either; you'll need to step through your code and check that each variable seems reasonable (correct size, reasonable values) as it runs, and find where it goes wrong and fix it.

[Vasco]

i checked the values, the values of the other variables check, i just dont get what is happening maybe is it because the data doesnt come sorted?

the figure im making is an attempt to answer this question :"Assess the differences between ETRS89 and the local datums (D73 and DLx) regarding geodesic coordinates (Δφ, Δλ, Δh)."

Answer 1

[Cris LaPierre]

I'd use scatteredInterpolant probably.

Also, you can simplify your code a little by accessing the data in your tables more directly. See here: https://www.mathworks.com/help/matlab/matlab_prog/access-data-in-a-table.html

Here's the relevant pieces of code. Note that I am not claiming to solve your assignment for you. Just get you started.

%Inserção dos dados e criação de variáveis
file_in=readtable("RGN_ETRS89_geod.txt");

lat_dec = dms2degrees(file_in{:,["Latd","Latm","Lats"]});
lon_dec = -1*dms2degrees(file_in{:,["Lond","Lonm","Lons"]});
alt_orto = file_in.h;
alt_elips = file_in.H;

% scatterm(lat_dec,lon_dec)
k = boundary(lat_dec,lon_dec);
geoshow(lat_dec(k),lon_dec(k))
%Elipsoide de referência para ETRS89: GRS80
a=6378137;
f=1/298.257222101;
b=a*(1-f);
e=(sqrt(a^2-b^2))/a;

N=a./sqrt(1-e^2.*(sind(lat_dec)).^2);

%Cálculo de X, Y e Z
X=(N+alt_elips).*cosd(lat_dec).*cosd(lon_dec);
Y=(N+alt_elips).*cosd(lat_dec).*sind(lon_dec);
Z=(N.*(1-e^2)+alt_elips).*sind(lat_dec);

%Exercício 2
%2.1
Tx=230.994;
Ty=-102.591;
Tz=-25.199;
omegax=deg2rad(-0.633/3600);
omegay=deg2rad(0.239/3600);
omegaz=deg2rad(-0.900/3600);
D=-1.950*10^-6;
matriz=[1 -omegaz omegay; omegaz 1 -omegax; -omegay omegax 1];
for i=1:1668
    D73=[Tx;Ty;Tz]+(D+1)*matriz*[X(i);Y(i);Z(i)];
    XD73(i,1)=D73(1);
    YD73(i,1)=D73(2);
    ZD73(i,1)=D73(3);
end

%2.2
Tx=283.088;
Ty=70.693;
Tz=-117.445;
omegax=deg2rad(1.157/3600);
omegay=deg2rad(-0.059/3600);
omegaz=deg2rad(0.652/3600);
D=4.058*10^-6;
matriz=[1 -omegaz omegay; omegaz 1 -omegax; -omegay omegax 1];
for i=1:1668
    DLx=[Tx;Ty;Tz]+(D+1)*matriz*[X(i);Y(i);Z(i)];
    XDLx(i,1)=DLx(1);
    YDLx(i,1)=DLx(2);
    ZDLx(i,1)=DLx(3);
end

%Elipsóide de referência para o Datum 73 e Datum Lisboa - Hayford
a_D73= 6378388;
f_D73=1/297;
b_D73=a_D73*(1-f_D73);
e_D73=(sqrt(a_D73^2-b_D73^2))/a_D73;

%Datum Lisboa
%Determinação da latitude
%Estimativa inicial
p=sqrt(XDLx.^2+YDLx.^2);
lat_0=atand(ZDLx./p);
N_0=a_D73./sqrt(1-e_D73^2.*(sind(lat_0)).^2);

%Primeira iteração
lat_1=atand((ZDLx+e_D73^2.*N_0.*sind(lat_0))./p);
N_1=a_D73./sqrt(1-e_D73^2.*(sind(lat_1)).^2);

%Segunda iteração
lat_2=atand((ZDLx+e_D73^2.*N_1.*sind(lat_1))./p);
N_2=a_D73./sqrt(1-e_D73^2.*(sind(lat_2)).^2);

%Terceira iteração
lat_3_DLx=atand((ZDLx+e_D73^2.*N_2.*sind(lat_2))./p);

%DLX figura e D73 figura
delta_lat_Dlx = lat_dec - lat_3_DLx;

% Create a function for predicting delta_lat_Dlx based on lat & lon data
F = scatteredInterpolant(lon_dec,lat_dec,delta_lat_Dlx);

% Create a 2d grid of Lat and Lon values
[LON,LAT] = meshgrid(-10:0.1:-6,36:0.05:43);
% Predict gridded values for delta_lat_Dlx
Z = F(LON,LAT);

% set values outside portugal to nan (so they do not appear in contour plot
in = inpolygon(LAT,LON,lat_dec(k),lon_dec(k));
Z(~in) = nan;

% Create a contour plot
hold on
contourm(LAT,LON,Z,10,'r','ShowText','on')
hold off

[Vasco]

First of all thank you for your suport.

I tried doing it for the lat, lon and h

the figure appears but the limits of the country dont Im gonna send the full code again if you could help me fix the problem, it would be live saving.

close all;
clear all;
clc;
%Exercício 1
%Inserção dos dados e criação de variáveis
file_in=readtable("RGN_ETRS89_geod.txt");


lat_dec = dms2degrees(file_in{:,["Latd","Latm","Lats"]});
lon_dec = -1*dms2degrees(file_in{:,["Lond","Lonm","Lons"]});
alt_orto = file_in.h;
alt_elips = file_in.H;

%Transformação dos dados em decimais
k = boundary(lat_dec,lon_dec);
geoshow(lat_dec(k),lon_dec(k));
%Elipsoide de referência para ETRS89: GRS80
a=6378137;
f=1/298.257222101;
b=a*(1-f);
e=(sqrt(a^2-b^2))/a;

%Cálculo de N (grande normal)
N=a./sqrt(1-e^2.*(sind(lat_dec)).^2);

%Cálculo de X, Y e Z
X=(N+alt_elips).*cosd(lat_dec).*cosd(lon_dec);
Y=(N+alt_elips).*cosd(lat_dec).*sind(lon_dec);
Z=(N.*(1-e^2)+alt_elips).*sind(lat_dec);
%% 
%Exercício 2
%2.1
Tx=230.994;
Ty=-102.591;
Tz=-25.199;
omegax=deg2rad(-0.633/3600);
omegay=deg2rad(0.239/3600);
omegaz=deg2rad(-0.900/3600);
D=-1.950*10^-6;
matriz=[1 -omegaz omegay; omegaz 1 -omegax; -omegay omegax 1];
for i=1:1668
    D73=[Tx;Ty;Tz]+(D+1)*matriz*[X(i);Y(i);Z(i)];
    XD73(i,1)=D73(1);
    YD73(i,1)=D73(2);
    ZD73(i,1)=D73(3);
end
%%
%2.2
Tx=283.088;
Ty=70.693;
Tz=-117.445;
omegax=deg2rad(1.157/3600);
omegay=deg2rad(-0.059/3600);
omegaz=deg2rad(0.652/3600);
D=4.058*10^-6;
matriz=[1 -omegaz omegay; omegaz 1 -omegax; -omegay omegax 1];
for i=1:1668
    DLx=[Tx;Ty;Tz]+(D+1)*matriz*[X(i);Y(i);Z(i)];
    XDLx(i,1)=DLx(1);
    YDLx(i,1)=DLx(2);
    ZDLx(i,1)=DLx(3);
end
%%
%Elipsóide de referência para o Datum 73 e Datum Lisboa - Hayford
a_D73= 6378388;
f_D73=1/297;
b_D73=a_D73*(1-f_D73);
e_D73=(sqrt(a_D73^2-b_D73^2))/a_D73;

%Cálculo da longitude
lon_D73=atand(YD73(:,1)./XD73(:,1));

%Determinação da latitude
%Estimativa inicial
p=sqrt((XD73(:,1)).^2+(YD73(:,1)).^2);
lat_0_73=atand(ZD73(:,1)./p);
N_0_73=a_D73./sqrt(1-e_D73^2.*(sind(lat_0_73)).^2);

%Primeira iteração
lat_1_73=atand((ZD73(:,1)+e_D73^2.*N_0_73.*sind(lat_0_73))./p);
N_1_73=a_D73./sqrt(1-e_D73^2.*(sind(lat_1_73)).^2);

%Segunda iteração
lat_2_73=atand((ZD73(:,1)+e_D73^2.*N_1_73.*sind(lat_1_73))./p);
N_2_73=a_D73./sqrt(1-e_D73^2.*(sind(lat_2_73)).^2);

%Terceira iteração
lat_3_D73=atand((ZD73(:,1)+e_D73^2.*N_2_73.*sind(lat_2_73))./p);
N_3_73=a_D73./sqrt(1-e_D73^2.*(sind(lat_3_D73)).^2);
%Os valores já não se alteram da segunda para a terceira iteração, fazendo
%desta a que fornece o valor definitivo da latitude.

%Cálculo de h
h_D73=(p./cosd(lat_3_D73))-N_3_73;
%%
%Datum Lisboa
%Cálculo da longitude
lon_DLx=atand(YDLx./XDLx);

%Determinação da latitude
%Estimativa inicial
p=sqrt(XDLx.^2+YDLx.^2);
lat_0_lx=atand(ZDLx./p);
N_0_lx=a_D73./sqrt(1-e_D73^2.*(sind(lat_0_lx)).^2);

%Primeira iteração
lat_1_lx=atand((ZDLx+e_D73^2.*N_0_lx.*sind(lat_0_lx))./p);
N_1_lx=a_D73./sqrt(1-e_D73^2.*(sind(lat_1_lx)).^2);

%Segunda iteração
lat_2_lx=atand((ZDLx+e_D73^2.*N_1_lx.*sind(lat_1_lx))./p);
N_2_lx=a_D73./sqrt(1-e_D73^2.*(sind(lat_2_lx)).^2);

%Terceira iteração
lat_3_DLx=atand((ZDLx+e_D73^2.*N_2_lx.*sind(lat_2_lx))./p);
N_3_lx=a_D73./sqrt(1-e_D73^2.*(sind(lat_3_DLx)).^2);
%Os valores já não se alteram da segunda para a terceira iteração, fazendo
%desta a que fornece o valor definitivo da latitude.

%Cálculo de h
h_DLx=(p./cosd(lat_3_DLx))-N_3_lx;
%%

%----------------
%----- etrs89 and datum lisboa
%DLX figura e D73 figura

delta_lat_Dlx = (lat_dec - lat_3_DLx).*3600;
delta_lon_Dlx = (lon_dec - lon_D73).*3600;
delta_h_Dlx = (alt_elips-h_D73);
%%
figure(1);
% Create a function for predicting delta_lat_Dlx based on lat & lon data
F = scatteredInterpolant(lon_dec,lat_dec,delta_lat_Dlx);
% Create a 2d grid of Lat and Lon values
[LON,LAT] = meshgrid(-10:0.1:-6,36:0.05:43);
% Predict gridded values for delta_lat_Dlx
Z = F(LON,LAT);
% set values outside portugal to nan (so they do not appear in contour plot
in = inpolygon(LAT,LON,lat_dec(k),lon_dec(k));
Z(~in) = nan;
% Create a contour plot
hold on
contourm(LAT,LON,Z,10,'r','ShowText','on')
hold off

%%
figure(2);
% Create a function for predicting delta_lat_Dlx based on lat & lon data
G = scatteredInterpolant(lon_dec,lat_dec,delta_lon_Dlx);
% Create a 2d grid of Lat and Lon values
[LON_1,LAT_1] = meshgrid(-10:0.1:-6,36:0.05:43);
% Predict gridded values for delta_lat_Dlx
V = G(LON_1,LAT_1);
% set values outside portugal to nan (so they do not appear in contour plot
in_1 = inpolygon(LAT_1,LON_1,lat_dec(k),lon_dec(k));
V(~in_1) = nan;
% Create a contour plot
hold on
contourm(LAT_1,LON_1,V,10,'r','ShowText','on')
hold off

%%
figure(3);
% Create a function for predicting delta_lat_Dlx based on lat & lon data
L = scatteredInterpolant(lon_dec,lat_dec,delta_h_Dlx);
% Create a 2d grid of Lat and Lon values
[LON_2,LAT_2] = meshgrid(-10:0.1:-6,36:0.05:43);
% Predict gridded values for delta_lat_Dlx
O = L(LON_2,LAT_2);
% set values outside portugal to nan (so they do not appear in contour plot
in = inpolygon(LAT_2,LON_2,lat_dec(k),lon_dec(k));
O(~in_2) = nan;

Unrecognized function or variable 'in_2'.

% Create a contour plot
hold on
contourm(LAT_2,LON_2,O,10,'r','ShowText','on')
hold off

[Cris LaPierre]

You need to redraw the border in each figure. You currently only draw it for the first one.

Also, the error message indicates that you do not have a variable 'in_2'

[Vasco]

how do i redraw the border?

Ive seem to forget to put the variable

[Cris LaPierre]

The figure command creates a new, empty figure. Identfy the code that draws the border, and repeat it after each figure command

[Vasco]

thank you very much its fixed

i try to do it with subplot but theres some type of error because it loads the lines but it doesnt load the values nor the borders.

[Cris LaPierre]

I wouldn't expect subplot to fix the issue. Please try what I suggested.