Спойлер
Код: Выделить всё
clear
clc
[Z_UP, R_UP] = geotiffread('srtm_49_01.tif'); %считываем первый растр
[Z_DOWN, R_DOWN] = geotiffread('srtm_49_02.tif'); %считываем второй растр
GRID_SIZE=size(Z_UP,1);
map=[Z_UP;Z_DOWN]; % склеиваем их в одну матрицу
mapref=R_UP; % сохраняем матрицу привязки растра
if GRID_SIZE==6000
mapref.RasterSize=[12000 6000]; % меняем размер, потому что теперь у нас не два растра 6000x6000, а один размером 12000x6000
mapref.LatitudeLimits=[50 60]; %корректируем значения границ по широте
mapref.LongitudeLimits=[60 65]; %корректируем значения границ по широте
else
mapref.RasterSize=[12002 6001]; % меняем размер, потому что теперь у нас не два растра 6001x6001, а один размером 12002x6001
mapref.LatitudeLimits=[50 60]; %корректируем значения границ по широте
mapref.LongitudeLimits=[60 65]; %корректируем значения границ по широте
end
%Точка начала профиля
lat1=54.99889;
lon1=60.37778;
%Точка окончания профиля
lat2=55.12726;
lon2=60.57104;
count_points=100;
lats=linspace(lat1,lat2,count_points);
longs=linspace(lon1,lon2,count_points);
val = ltln2val(double(map),mapref,lats,longs,'bilinear');
dist=distance(lat1,lon1,lat2,lon2,6378137,'degree');
dist=linspace(0,dist,count_points);
X=dist;
Y=val;
lambda=0.75;
D=sqrt((X(1)-X(end))^2+(Y(1)-Y(end))^2);
r=linspace(0,D,count_points);
for ii=1:length(r)
y(ii)=sqrt((1/5)*(lambda*r(ii)*(D-r(ii)))*(1/D));
x(ii)=r(ii);
end
% FY=[y -y];
% FX=[x x];
% FY=FY+Y(1)+20;
FX=x;
FY=y;
FY=FY+Y(1)+20;
x_center=X(1);
y_center=Y(1)+20;
b1=glmfit([X(1) X(end)],[Y(1) Y(end)]);
theta=atan(b1(2));
R = [cos(theta) -sin(theta);
sin(theta) cos(theta)];
v=[FX;FY];
center = repmat([x_center; y_center], 1, length(FX));
vo = R*(v - center) + center;
x_rotated = vo(1,:);
y_rotated = vo(2,:);
% Вырезаем из SRTM кусок, чтобы построить 3D поверхность
[X,Y] = meshgrid(longs,lats);
Z = ltln2val(double(map),mapref,Y,X);
dist=distance(lat1,lon1,lat1,lon2,6378137,'degree');
Y=linspace(0,dist,count_points);
dist=distance(lat1,lon1,lat2,lon1,6378137,'degree');
X=linspace(0,dist,count_points);
[X,Y] = meshgrid(X,Y);
figure('render','opengl');
surf(X,Y,Z);
hold on
grid on
for ii=1:count_points
X_C(ii)=X(1,ii);
Y_C(ii)=Y(ii,end);
end
plot3(X_C,Y_C,val+1,'-r','LineWidth',2);
plot3(X_C,Y_C,linspace(y_rotated(1),y_rotated(end),count_points),'--k','LineWidth',2);
CONV=0.0174532925;
hip=sqrt(X_C(end)^2+Y_C(end)^2);
angel_z=asin(Y_C(end)/hip)/CONV;
V=[];
FY=FY-Y(1)-20;
L=min(FY);
FY=FY-L;
v=[FX;zeros(1,length(FX));FY];
%Формируем 3D поверхность по элипсу Френеля
% Она будет выглядеть сплющенной с боков, но это норм...
% У нас высота идёт в диапазоне 340-380 метров, а по длине идут 100 метров
% в итоге, чтобы хоть что-то было понятно, приходится рисовать так.
for ii=1:359
R=rotx(ii);
v=R*v;
V=[V v];
end
b1=glmfit([X_C(1,1) X_C(end,end)],[y_rotated(1) y_rotated(end)]);
%:Доворачиваем по осям X и Y
R = rotz(angel_z);
V=R*V;
R=roty(-atan(b1(2))/CONV);
V=R*V;
V(3,:)=V(3,:)+Y(1)+20+L;
plot3(V(1,:),V(2,:),V(3,:),'-b','LineWidth',2);