Matlab profil çizimini konuşuyor olacağız. Sayısal yükseklik modeli üzerinden nasıl profil çizdirebileceğinizi anlatacağım.
Matlab Profil Sonuç:
Merhaba,
Matlab ile Profil (Boykesit) Çizimi : Improfile Fonksiyonu
Matlab profil çizimlerimiz için oldukça kullanışlı araçlar sunmakta. Burada kullanacağımız ve yazımızın odak noktası olacak fonksiyonumuz improfile fonksiyonu. Şimdi uygulamamızı sırası ile kodları açıklayarak anlatalım ve daha sonra toplu bir biçimde kodları size veriyor olacağım. Öncelikle mantığını anlamak bizim için daha önemli.
Başlamadan: Eğer konuya hakimseniz daha detaylı bir uygulama için aşağıdaki konuya bakabilirsiniz.
Matlab – SYM ve ORTOFOTO ile Hacim Hesabı Uygulaması (Giriş)
1. Adım: Ortofoto ve Sayısal Yükseklik Modellerinin Okunması
% sym dosyası seçilir.
[file,path] = uigetfile('*.tif','Lütfen SYM Seçiniz');
SYM=fullfile(path,file);
% ortofoto dosyası seçilir.
[file,path] = uigetfile('*.tif','Lütfen Ortofoto Seçiniz');
ORTO=fullfile(path,file);
%sym geotif olarak okunmalı
[gtrSYM, gtrRS]=geotiffread(SYM);
[gtrORT, gtrRO]=geotiffread(ORTO);
Yukarıda ki kodu açıklayacak olursak eğer öncelikle uigetfile diyerek doyalarımızın yolunu kullanıcıya seçtirdik. Daha sonra gelen iki parametreyi fullfile fonksiyonu ile birleştirdik. Zira uigetfile fonksiyonu file değişkenine dosya ismini, path değişkenine ise dosya yolunu alıyor. Bunu dosyayolu/dosyaismi şekline getirmemiz lazım ki düzgünce kullanabilelim. Bunu da işte fullfile fonksiyonu ile yapmaktayız.
geotiffread diyerek sayısal yükseklik modelini ve ortofoto haritayı okuyoruz. Burada fonksiyon bize iki adet veri getiriyor. gtrSYM olarak adlandırdığım kısımda verinin daha çok görsel ve piksel değerlerini getirirken diğer kısımda ise jeodezik özelliklerini (koordinat bilgileri, piksel boyutları vb.) getiriyor.
2. Adım : Okunan Ortofotonun Ekranda Gösterilmesi
%sym ekrana getirilmeli imshow(ORTO); hold on
Bu kısımda harita gösterme komutlarıda kullanılabilirdi fakat ben burada imshow ile göstermeyi tercih ediyorum. Zira diğer plot vb. gibi fonksiyonların kullanımında benim işimi kolaylaştırıyor. Harita koordinatlarına nasıl geçtiğimi zaten gösteriyor olacağım.
hold on diyerek görüntüyü sürekli ekranda tutuyoruz.
3. Adım: Profil Başlangıç ve Bitiş Noktalarının Seçimi
i=1;
while i<3
[x,y,button] = ginput(1);
if button==27
break;
end
if button==3
fx=fix(x);
fy=fix(y);
hA(i,1)=gtrSYM(fy,fx);
X(i,1)=x;
Y(i,1)=y;
plot(X,Y,'--gs',...
'LineWidth',2,...
'MarkerSize',10,...
'MarkerEdgeColor','b',...
'MarkerFaceColor',[0.5 0.5 0.5]);
i=i+1;
end
end
Burada ilk etapda bir döngü kullanıyorum. ginput(2) diyerek iki adet nokta seçmesini isteyebilirdim fakat ben eğer daha farklı özellikler eklemek istiyorsam, örneğin esc tuşuna bastığında çizimden çıkması gibi, bu işi döngü ile yapmalıyım. Burada döngünün mantığı oldukça basit. Döngü i değişkeni 0’dan büyük olduğu sürece devam edecek ve ben i değişkeni sıfırdan büyük olduğu sürece nokta atabileceğim.
ginput ile biz 3 veri alabiliyoruz. x, y ve buton değeri. Burada x ve y’nin ne olduğunu zaten biliyoruz. Buton ise o ginput kullanımı esnasında mouse ile bir yere tıklanmayıp klavyeden bir tuşa basılmışsa eğer onun değerini alacaktır. Mouse ile tıkladığınızda ise buton mousenin sol tuşu için 1 sağ tuşu için 3 değerini alır.
Yukarıdaki kodu genel olarak anlatacak olursak eğer şöyle diyebiliriz. Ben esc’ye basarsam programdan çık. if button=27 olan blok bunu yapıyor. Eğer ben ancak mousenin sağ tuşuna basarsam çizim yap diyorum ki bu da if button=3 olan kısım oluyor. Mousenin sol tuşunu serbest bıraktım ki gui ekranında rahatça oraya buraya tıklayıp ekranı kaydırabilin ve yaklaştırabilin.
Günün sonunda bu kısımda elimizde x ve y fotoğraf koordinatları olmuş oluyor.
4. Adım: Improfile Fonksiyonu ile Gerekli Verileri Alma
[matlab smarttabs=”true”][cx,cy,c,xi,yi]=improfile(gtrSYM,X,Y);[/matlab]
Fonksiyonun kullanımı bu kadar arkadaşlar. Burada gtrSYM(sym pixelleri içerisindeki yükseklik değerlerini) ‘den seçtiğimiz hat boyunca olan X ve Y koordinatlarına göre yükseklik verisini alıyoruz.
5. Adım: Fotoğraf koordinatlarından Jeodezik Koordinatlara Geçiş
for i=1:length(cx)
worldY(i,1)=gtrRO.XWorldLimits(1,1)+(gtrRO.CellExtentInWorldX(1,1)*cx(i,1));
worldX(i,1)=gtrRO.YWorldLimits(1,2)-(gtrRO.CellExtentInWorldY(1,1)*cy(i,1));
worldMes(i,1)=i*0.02;
end
Burada bilmemiz gereken gtrRO.XWorldLimits gibi ifadelerin ne olduğudur. Aslında tahmin etmiş olmanız gerekiyor. Arkadaşlar bu ifadeler seçilen ortofotonun sol üst ve sağ alt köşe koordinatları. CellExtentInWorldX ise arakdaşlar pixel boyutlarıdır.
Burada yapılan şey, köşe koordinatlarına seçilen kadar pixeli, piksel boyutu ile çarp ve çıkar yada topla. Amacımız fotoğraf koordinatlarından jeodezik koordinatlara geçmek. Eğer köşe koordinatlarını biliyorsam ve piksel boyutlarımında arazide ifade ettiği değer GSD değeri biliniyorsa bu işlem oldukça basit bir hal alıyor.
İHA Fotogrametrisi Açısından Odak Uzaklığına Kısa Bir Bakış GSD ile ilgili bilgiye bu makaleden ulaşabilirsiniz. Matlab profil çizimlerimize geçelim.
6. Adım: Matlab Profil işlemlerinin Yapımı
figure(2)
subplot(2,1,1);
plot(worldMes,c);grid on, grid minor;
title("2-D Profil")
xlabel('Mesafe (m)')
ylabel('Yükseklik (m)')
subplot(2,1,2);
plot3(worldY,worldX,c);grid on, grid minor;
title("3-D Profil")
xlabel('X')
ylabel('Y')
zlabel('Yükseklik (m)');
Kodlarımız burada son buluyor. Kodların tamamına aşağıdan ulaşabilirsiniz.
% sym dosyası seçilir.
[file,path] = uigetfile('*.tif','Lütfen SYM Seçiniz');
SYM=fullfile(path,file);
% ortofoto dosyası seçilir.
[file,path] = uigetfile('*.tif','Lütfen Ortofoto Seçiniz');
ORTO=fullfile(path,file);
%sym geotif olarak okunmalı
[gtrSYM, gtrRS]=geotiffread(SYM);
[gtrORT, gtrRO]=geotiffread(ORTO);
%sym ekrana getirilmeli
imshow(ORTO);
hold on
i=1;
while i<3
[x,y,button] = ginput(1);
if button==27
break;
end
if button==3
fx=fix(x);
fy=fix(y);
hA(i,1)=gtrSYM(fy,fx);
X(i,1)=x;
Y(i,1)=y;
plot(X,Y,'--gs',...
'LineWidth',2,...
'MarkerSize',10,...
'MarkerEdgeColor','b',...
'MarkerFaceColor',[0.5 0.5 0.5]);
i=i+1;
end
end
[cx,cy,c,xi,yi]=improfile(gtrSYM,X,Y);
for i=1:length(cx)
worldY(i,1)=gtrRO.XWorldLimits(1,1)+(gtrRO.CellExtentInWorldX(1,1)*cx(i,1));
worldX(i,1)=gtrRO.YWorldLimits(1,2)-(gtrRO.CellExtentInWorldY(1,1)*cy(i,1));
worldMes(i,1)=i*0.02;
end
figure(2)
subplot(2,1,1);
plot(worldMes,c);grid on, grid minor;
title("2-D Profil")
xlabel('Mesafe (m)')
ylabel('Yükseklik (m)')
subplot(2,1,2);
plot3(worldY,worldX,c);grid on, grid minor;
title("3-D Profil")
xlabel('X')
ylabel('Y')
zlabel('Yükseklik (m)');
Sorularınızı yorum kısmından sorabilirsiniz.
Faydalı olması dileği ile herkese iyi çalışmalar diliyorum.
