Tugas 2 - Aplikasi Geodesi Satelit - Nabil Amirul Haq - 6016202002 [PDF]

Citation preview

TUGAS APLIKASI GEODESI SATELIT



“Mencari GAST dan Transformasi Koordinat Titik pada Scene Citra”



Dosen: Mokhamad Nur Cahyadi, S.T., M.Sc., Ph.D



Disusun Oleh: Nabil Amirul Haq



(6016202002)



PROGRAM MAGISTER TEKNIK GEOMATIKA FAKULTAS TEKNIK SIPIL, PERENCANAAN, DAN KEBUMIAN INSTITUT TEKNOLOGI SEPULUH NOPEMBER SURABAYA 2021



I.



MENGHITUNG GAST



Diketahui : 1. State Vector Time Sentinel 1A: 2021-Maret-12 11:41:18.32035 Dicari: GAST pada state vector time tersebut. Jawab: Script Matlab Perhitungan GAST clc;clear %Calculate the Julian Dates t=datetime('2021-03-12 11:41:18.32035') JD = juliandate(t) %THETAm is the mean siderial time in degrees THETAm = JD2GMST(JD); %Compute the number of centuries since J2000 T = (JD - 2451545.0)./36525; %Mean obliquity of the ecliptic (EPSILONm) % see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 3 EPSILONm = 84381448-46.8150.*T - .00059.*(T.^2) + .001813.*(T.^3); %Nutations in obliquity and longitude (degrees) % see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 4 L = 280.4665 + 36000.7698.*T; dL = 218.3165 + 481267.8813.*T; OMEGA = 125.04452 - 1934.136261.*T; %Calculate nutations using the following two equations: % see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 5 dPSI = -17.20.*sind(OMEGA) - 1.32.*sind(2.*L) - .23.*sind(2.*dL) ... + .21.*sind(2.*OMEGA); dEPSILON = 9.20.*cosd(OMEGA) + .57.*cosd(2.*L) + .10.*cosd(2.*dL) - ... .09.*cosd(2.*OMEGA); %Convert the units from arc-seconds to degrees dPSI = dPSI.*(1/3600); dEPSILON = dEPSILON.*(1/3600); %(GAST) Greenwhich apparent sidereal time expression in degrees % see http://www.cdeagle.com/ccnum/pdf/demogast.pdf equation 1 GASTdd = mod(THETAm + dPSI.*cosd(EPSILONm+dEPSILON),360) GAST = degrees2dms(GASTdd)



Dihasilkan nilai GAST dalam derajat 345.66710 atau 345040’1.7307”



II.



PROYEKSI TITIK PADA SCENE CITRA



Proyeksi titik pada scene citra Sentinel 1a akan dilakukan dengan memproyeksi koordinat geodetik pada bidang WGS84 ke koordinat UTM dan TM3. Untuk menentukan posisi titiknya, akan dilakukan menggunakan perangkat lunak ArcGIS. Citra yang digunakan pada tugas kali ini adalah citra yang sama dengan yang dipakai untuk tugas sebelumnya.



Gambar 1. Citra Sentinel 1a di ArcGIS Setelah itu, tambahkan 5 titik pada scene citra, di mana 1 titik pastinya adalah Gunung Sinabung, dan 4 lainnya posisinya bebas.



Gambar 2. Sebaran Titik pada Scene Berikut ini koordinat geodetik dari ke-5 titik tersebut. Tabel 1. Koordinat Titik Nomor Lat Long 1 3.0682246 98.44332 2 3.40038 97.70256 3 3.65379 98.70907 4 2.829309 99.54069 5 2.404576 98.09517



A. Proyeksi ke Bidang UTM Proyeksi dari koordinat Geodetik ke sistem proyeksi UTM dilakukan dengan menggunakan MATLAB, berikut ini script MATLAB yang penulis gunakan: clc;clear %parameter ellipsoid WGS84 a = 6378137 f = 298.257223563 b = a - (a*1/f) e2 = (a^2 -b^2)/a^2 e = sqrt((2*1/f)-(1/f^2)) k0 = 0.9996; n = (a-b)/(a+b) e22 = (a^2 -b^2)/b^2 ltg_dd = Nilai Lintang dalam Derajat Desimal; bjr_dd = Nilai Bujur dalam Derajat Desimal; ltg_rad = ltg_dd*pi/180; bujur_rad = bjr_dd*pi/180; N = a/((1-(e*sin(ltg_rad)^2))^0.5); zona = fix((bujur(1,1)+180)/6)+1; if bujur(1,1)>=0; B0 = (6*zona)-183; else B0 = -(6*zona)-183; end; disp ([' Koordinat terletak pada zona ' ,num2str(zona), ' dan meridian tengahnya adalah ',num2str(B0)]) p = abs(((bjr_dd*3600) - (B0*3600))*0.0001) disp('Menghitung Koefisien') sin1 = 4.84813681107637E-06 a0 = a*(1-n+(5*n*n/4)*(1-n)+(81*n^4/64)*(1-n)); b0 = (3*a*n/2)*(1-n-(7*n*n/8)*(1-n)+55*n^4/64); c0 = (15*a*n*n/16)*(1-n+(3*n*n/4)*(1-n)); d0 = (35*a*n^3/48)*(1-n+(11*n*n/16)); e0 = (315*a*n^4/51)*(1-n); m = a0*ltg_rad- b0*sin(2*ltg_rad)+c0*sin(4*ltg_rad)d0*sin(6*ltg_rad)+e0*sin(8*ltg_rad); I = k0*m II = k0*N*sin(ltg_rad)*cos(ltg_rad)*(sin1^2)*10^8/2 III = k0*N*sin(ltg_rad)*cos(ltg_rad)^3*(sin1^4)*(10^16)*(5tan(ltg_rad)^2+(9*e22)*cos(ltg_rad)^2 + 4*(e22^2)*cos(ltg_rad)^4)/24 IV = k0*N*cos(ltg_rad)*sin1*10^4 V = k0*N*cos(ltg_rad)^3*sin1^3*10^12*(1-tan(ltg_rad)^2+e22*cos(ltg_rad)^2)/6 A6 = k0*N*sin1^6*sin(ltg_rad)*cos(ltg_rad)^5*(6158*tan(ltg_rad)^2+tan(ltg_rad)^2+270*e22*cos(ltg_rad)^2330*e22*sin(ltg_rad)^2)*10^24/720 B5 = k0*N*cos(ltg_rad)^5*sin1^5*(5-18*tan(ltg_rad)^2+ tan(ltg_rad)^4 + 14*e22*cos(ltg_rad)^2 - 58*e22*sin(ltg_rad)^2)*10^20/120 disp(' Menghitung koordinat P sesuai dengan formula ') x = IV*p + V*p^3 + B5*p^5 y = I + II*p^2 + III*p^4 + A6*p^6 if lintang(1,1)>=0; Y = y; else Y = 10000000 - y; end; disp (['jadi nilai Y = ', num2str(Y) ]) if bujur(1,1) > B0; X = 500000 + x; else X = 500000 - x; end disp (['jadi nilai X = ', num2str(X) ]) disp ([' Koordinat titik ( ' num2str(X) ' , ' num2str(Y) ' ) '])



Berikut ini koordinat hasil proyeksi ke UTM: Tabel 2. Hasil Proyeksi UTM Nomor 1 2 3 4 5



X (m) 438136.0254 355848.5273 467687.0712 560098.8036 399392.2899



Y (m) 339150.3548 375945.4145 403864.3084 312740.2109 265812.8012



Berikut ini hasil plot koordinat UTM di Excel:



Gambar 3. Plot Koordinat UTM B. Proyeksi ke Bidang TM3 Proyeksi dari koordinat Geodetik ke sistem proyeksi TM3 dilakukan dengan menggunakan MATLAB, berikut ini script MATLAB yang penulis gunakan: clc;clear %parameter ellipsoid WGS84 a = 6378137 f = 298.257223563 b = a - (a*1/f) e2 = (a^2 -b^2)/a^2 e = sqrt((2*1/f)-(1/f^2)) k0 = 0.9999; n = (a-b)/(a+b) e22 = (a^2 -b^2)/b^2 ltg_dd = 2.404576; bjr_dd = 98.09517; bujur = degrees2dms(bjr_dd); lintang = degrees2dms(ltg_dd); ltg_rad = ltg_dd*pi/180; bujur_rad = bjr_dd*pi/180; N = a/((1-(e*sin(ltg_rad)^2))^0.5); zona = 47.1; B0 = 97.5;



disp ([' Koordinat terletak pada zona ' ,num2str(zona), ' dan meridian tengahnya adalah ',num2str(B0)]) p = abs(((bjr_dd*3600) - (B0*3600))*0.0001) disp('Menghitung Koefisien') sin1 = 4.84813681107637E-06 a0 = a*(1-n+(5*n*n/4)*(1-n)+(81*n^4/64)*(1-n)); b0 = (3*a*n/2)*(1-n-(7*n*n/8)*(1-n)+55*n^4/64); c0 = (15*a*n*n/16)*(1-n+(3*n*n/4)*(1-n)); d0 = (35*a*n^3/48)*(1-n+(11*n*n/16)); e0 = (315*a*n^4/51)*(1-n); m = a0*ltg_rad- b0*sin(2*ltg_rad)+c0*sin(4*ltg_rad)d0*sin(6*ltg_rad)+e0*sin(8*ltg_rad); I = k0*m II = k0*N*sin(ltg_rad)*cos(ltg_rad)*(sin1^2)*10^8/2 III = k0*N*sin(ltg_rad)*cos(ltg_rad)^3*(sin1^4)*(10^16)*(5tan(ltg_rad)^2+(9*e22)*cos(ltg_rad)^2 + 4*(e22^2)*cos(ltg_rad)^4)/24 IV = k0*N*cos(ltg_rad)*sin1*10^4 V = k0*N*cos(ltg_rad)^3*sin1^3*10^12*(1-tan(ltg_rad)^2+e22*cos(ltg_rad)^2)/6 A6 = k0*N*sin1^6*sin(ltg_rad)*cos(ltg_rad)^5*(6158*tan(ltg_rad)^2+tan(ltg_rad)^2+270*e22*cos(ltg_rad)^2330*e22*sin(ltg_rad)^2)*10^24/720 B5 = k0*N*cos(ltg_rad)^5*sin1^5*(5-18*tan(ltg_rad)^2+ tan(ltg_rad)^4 + 14*e22*cos(ltg_rad)^2 - 58*e22*sin(ltg_rad)^2)*10^20/120 disp(' Menghitung koordinat P sesuai dengan formula ') x = IV*p + V*p^3 + B5*p^5 y = I + II*p^2 + III*p^4 + A6*p^6 Y = 1500000 + y; disp (['jadi nilai Y = ', num2str(Y) ]) X = 200000 + x; disp (['jadi nilai X = ', num2str(X) ]) disp ([' Koordinat titik ( ' num2str(X) ' , ' num2str(Y) ' ) '])



Berikut ini koordinat hasil proyeksi ke TM3: Tabel 2. Hasil Proyeksi TM3 Nomor 1 2 3 4 5



X (m) 304865.9183 222510.2132 334338.3185 426939.6587 266195.0251



Y (m) 1839282.2586 1875963.7554 1904070.6253 1812864.1485 1765873.6599



Berikut ini hasil plotting koordinat TM3 di Excel:



Gambar 4. Plot Koordinat TM3