Serba-serbi tentang data ADCP, bagaimana membaca dan mengolahnya

Saya sedang punya proyek untuk menulis hal-hal penting terkait data ADCP dan bagaimana cara membaca dan mengolahnya di dalam Blog OSEANOGRAFI ini. Rencananya tulisan ini akan dibagi menjadi 2 bagian utama, yaitu:

Bagian 1: tentang jenis-jenis ADCP, sistem akuisi, dan karakteristik datanya.

Bagian 2: tentang cara membaca dan mengolah data ADCP serta perangkat lunak yang dapat digunakan untuk keperluan tersebut.

Sepertinya 2 bagian itu sudah cukup banyak. Atau mungkin ada hal lain yang perlu ditambahkan? silahkan kasih saya masukan ya kalau berkenan.

Proyek penulisan ini mungkin perlu waktu yang agak lama, tapi saya berjanji untuk secara disiplin dan konsisten mem-posting-nya di blog ini. Motivasinya sederhana saja. Motivasi pertama, saya ingin menghidupkan kembali blog yang sudah berdiri sejak tahuan 2000-an ini dan mengisinya dengan informasi-informasi tentang oseanografi yang mudah-mudahan dapat bermanfaat. Motivasi kedua, saya ingin meneruskan untuk belajar lebih detail lagi tentang ADCP yang sudah saya awali di tahun 2007 ketika mengikuti pelayaran di Samudera Pasifik dengan Kapal Riset Mirai dari JAMSTEC, Jepang. Motivasi ketiga, saya ingin membuat panduan teknis yang semoga bermanfaat untuk para calon Oseanografer Indonesia.

Akhir kata, mumpung masih hangat, saya ingin mengucapkan “Selamat Hari Nusantara, 13 Desember 2017!” Semoga semangat Deklarasi Djuanda selalu bergelora di dada kita semua. :)

Script MATLAB untuk Prediksi Elevasi Pasang Surut

Bongkar-bongkar harddisk dan laptop lama, ternyata saya menemukan sebuah script yang dibuat tahun 2011. Script ini dibuat untuk memprediksi elevasi pasang surut berdasarkan data amplitudo dan fasa dari 9 komponen pasang surut utama yang banyak digunakan di Indonesia. Daripada script ini hanya disimpan dan jamuran, berikut saya tuliskan saja di blog OSEANOGRAFI yang sudah jablay ini, siapa tahu bermanfaat buat mereka yang tidak sengaja berkunjung atau menemukan blog ini melalui mesin pencari.

Oh iya, di dalam script ini saya memanfaatkan code-code yang dibangun oleh Egbert dan Erofeeva (2002) dalam Oregon Tide Prediction Software (OTPS).

Enjoy!

% program prediksi elevasi pasang surut
% dengan input amplitudo (meter) dan fasa (derajat)
% untuk 9 komponen utama pasang surut
% m2,s2,k1,o1,n2,p1,k2
% agus setiawan - maret 2011

clear all

pi=3.14159265358979;
rad=pi/180.0;
ncmx=29;
nc=9; % jumlah komponen pasut untuk prediksi
ind=[1,2,3,4,5,6,7,21,22]; % nomor indeksnya dalam constid

% constid
% 'm2 ','s2 ','k1 ','o1 ',
% 'n2 ','p1 ','k2 ','q1 ',
% '2n2 ','mu2 ','nu2 ','l2 ',
% 't2 ','j1 ','m1 ','oo1 ',
% 'rho1','mf ','mm ','ssa ',
% 'm4 ','ms4 ','mn4 ','m6 ',
% 'm8 ','mk3 ','s6 ','2sm2',
% '2mk3'

% example data di stasiun tg. priok

% MSL
s0 = 0.0; % in meter

% amplitudo dan fasa
ampi(1) = 0.074; phasi(1) = 174.10; % m2
ampi(2) = 0.025; phasi(2) = 117.10; % s2
ampi(3) = 0.262; phasi(3) = 32.70; % k1
ampi(4) = 0.125; phasi(4) = 16.80; % o1
ampi(5) = 0.020; phasi(5) = 109.00; % n2
ampi(6) = 0.072; phasi(6) = 21.60; % p1
ampi(7) = 0.009; phasi(7) = 56.00; % k2
ampi(8) = 0.00; phasi(8) = 0.00; % m4
ampi(9) = 0.00; phasi(9) = 0.00; %ms4

% convert x and y components
for i=1:nc,
xx(i)=ampi(i)*cos(phasi(i)*rad);
yy(i)=ampi(i)*sin(phasi(i)*rad);
end

omega_d=[1.405189e-04,1.454441e-04,7.292117e-05,6.759774e-05,...
1.378797e-04,7.252295e-05,1.458423e-04,6.495854e-05,...
1.352405e-04,1.355937e-04,1.382329e-04,1.431581e-04,...
1.452450e-04,7.556036e-05,7.028195e-05,7.824458e-05,...
6.531174e-05,0.053234e-04,0.026392e-04,0.003982e-04,...
2.810377e-04,2.859630e-04,2.783984e-04,4.215566e-04,...
5.620755e-04,2.134402e-04,4.363323e-04,1.503693e-04,...
2.081166e-04];

phase_mkB=[1.731557546,0.000000000,0.173003674,1.558553872,...
6.050721243,6.110181633,3.487600001,5.877717569,...
4.086699633,3.463115091,5.427136701,0.553986502,...
0.052841931,2.137025284,2.436575100,1.929046130,...
5.254133027,1.756042456,1.964021610,3.487600001,...
3.463115091,1.731557546,1.499093481,5.194672637,...
6.926230184,1.904561220,0.000000000,4.551627762,...
3.809122439];

index=[30,35,19,12,27,17,37,10,25,26,28,33,34,...
23,14,24,11,5,3,2,45,46,44,50,0,42,51,40,0];

% awal tanggal peramalan
iyyy = 2010; % year
mm = 3; % month
id = 1; % day
hh = 0; % hour
mi = 0; % minute
ss = 0; % second
tgl0=datenum(iyyy,mm,id,hh,mi,ss);

% convert date to modified julian day
dpm(1) = 31; % jan
dpm(2) = 28; % feb
dpm(3) = 31; % mar
dpm(4) = 30; % apr
dpm(5) = 31; % may
dpm(6) = 30; % jun
dpm(7) = 31; % jul
dpm(8) = 31; % aug
dpm(9) = 30; % sep
dpm(10)= 31; % oct
dpm(11)= 30; % nov
dpm(12)= 31; % dec

mjd = 0;
% no earlier dates
if(iyyy<1858),iyyy=1858; end if(iyyy==1858 & mm>11),mm=11;end
if(iyyy==1858 & mm==11 & id>17),id=17;end

days=0;
for i=1:mm-1,
   days=days+dpm(i);
   if(i==2 & fix(iyyy/4)*4==iyyy),days=days+1;end
end
days=days+id-321;

% leap day correction
for k=1900:100:iyyy,
   if(iyyy==k & mm>2),days=days-1;end
end

for k=2000:400:iyyy,
   if(iyyy==k & mm>2),days=days+1;end
end

% each 4th year is leap year
nleap=fix((iyyy-1-1860)*0.25);
if(iyyy>1860),nleap=nleap+1;end

% except
for k=1900:100:iyyy-1,
   if(iyyy>k),nleap=nleap-1;end
   if(iyyy==k & mm>2),days=days-1;end
end

% but each in the row 2000:400:... is leap year again
for k=2000:400:iyyy-1,
   if(iyyy>k),nleap=nleap+1;end
   if(iyyy==k & mm>2),days=days+1;end
end
mjd=365*(iyyy-1858)+nleap+days;
time_mjd=mjd+hh/24.0+mi/(24.0*60.0)+ss/(24.0*60.0*60.0);

ntime=31*24*60; % lama prediksi nday*24hr*60 (in minute)
for k=1:1:ntime,
   tgl(k)=tgl0+k/1440;

   % astrol
   T=((k-1)/(24*60)+time_mjd)-51544.4993;
   circle=360.0;

   SHPN(1)=218.3164+13.17639648*T;
   SHPN(2)=280.4661+0.98564736*T;
   SHPN(3)=83.3535+0.11140353*T;
   SHPN(4)=125.0445-0.05295377*T;

   SHPN(1)=mod(SHPN(1),circle);
   SHPN(2)=mod(SHPN(2),circle);
   SHPN(3)=mod(SHPN(3),circle);
   SHPN(4)=mod(SHPN(4),circle);

   if(SHPN(1)<0.0),SHPN(1)=SHPN(1)+circle;end
   if(SHPN(2)<0.0),SHPN(2)=SHPN(2)+circle;end
   if(SHPN(3)<0.0),SHPN(3)=SHPN(3)+circle;end
   if(SHPN(4)<0.0),SHPN(4)=SHPN(4)+circle;end s=SHPN(1); 
   h=SHPN(2); 
   p=SHPN(3); 
   omega=SHPN(4); 
   pp=282.94; 

   % nodal 
   hour=(((k-1)/(24.0*60)+time_mjd)-fix(((k-1)/(24.0*60)+time_mjd)))*(24.0*60); 
   t1=15.0*hour; t2=30.0*hour; 
   arg(1)=h-pp; % Sa 
   arg(2)=2.0*h; % Ssa 
   arg(3)=s-p; % Mm 
   arg(4)=2.0*s-2.0*h; % MSf 
   arg(5)=2.0*s; % Mf 
   arg(6)=3.0*s-p; % Mt 
   arg(7)=t1-5.0*s+3.0*h+p-90.0; % alpha1 
   arg(8)=t1-4.0*s+h+2.0*p-90.0; % 2Q1 
   arg(9)=t1-4.0*s+3.0*h-90.0; % sigma1 
   arg(10)=t1-3.0*s+h+p-90.0; % q1 
   arg(11)=t1-3.0*s+3.0*h-p-90.0; % rho1 
   arg(12)=t1-2.0*s+h-90.0; % o1 
   arg(13)=t1-2.0*s+3.0*h+90.0; % tau1 
   arg(14)=t1-s+h+90.0; % M1 
   arg(15)=t1-s+3.0*h-p+90.0; % chi1 
   arg(16)=t1-2.0*h+pp-90.0; % pi1 
   arg(17)=t1-h-90.0; % p1 
   arg(18)=t1+90.0; % s1 
   arg(19)=t1+h+90.0; % k1 
   arg(20)=t1+2.0*h-pp+90.0; % psi1 
   arg(21)=t1+3.0*h+90.0; % phi1 
   arg(22)=t1+s-h+p+90.0; % theta1 
   arg(23)=t1+s+h-p+90.0; % J1 
   arg(24)=t1+2.0*s+h+90.0; % OO1 
   arg(25)=t2-4.0*s+2.0*h+2.0*p; % 2N2 
   arg(26)=t2-4.0*s+4.0*h; % mu2 
   arg(27)=t2-3.0*s+2.0*h+p; % n2 
   arg(28)=t2-3.0*s+4.0*h-p; % nu2 
   arg(29)=t2-2.0*s+h+pp; % M2a 
   arg(30)=t2-2.0*s+2.0*h; % M2 
   arg(31)=t2-2.0*s+3.0*h-pp; % M2b 
   arg(32)=t2-s+p+180.0; % lambda2 
   arg(33)=t2-s+2.0*h-p+180.0; % L2 
   arg(34)=t2-h+pp; % t2 
   arg(35)=t2; % 
   S2 arg(36)=t2+h-pp+180.0; % R2 
   arg(37)=t2+2.0*h; % K2 
   arg(38)=t2+s+2.0*h-pp; % eta2 
   arg(39)=t2-5.0*s+4.0*h+p; % MNS2 
   arg(40)=t2+2.0*s-2.0*h; % 2SM2 
   arg(41)=1.5*arg(30); % M3 
   arg(42)=arg(19)+arg(30); % MK3 
   arg(43)=3.0*t1; % S3 
   arg(44)=arg(27)+arg(30); % MN4 
   arg(45)=2.0*arg(30); % M4 
   arg(46)=arg(30)+arg(35); % MS4 
   arg(47)=arg(30)+arg(37); % MK4 
   arg(48)=4.0*t1; % S4 
   arg(49)=5.0*t1; % S5 
   arg(50)=3.0*arg(30); % M6 
   arg(51)=3.0*t2; % S6 
   arg(52)=7.0*t1; % S7 
   arg(53)=4.0*t2; 

   sinn=sin(omega*rad); 
   cosn=cos(omega*rad); 
   sin2n=sin(2.0*omega*rad); 
   cos2n=cos(2.0*omega*rad); 
   sin3n=sin(3.0*omega*rad); 

   f(1)=1.0; % Sa 
   f(2)=2.0; % Ssa 
   f(3)=1.0-0.130*cosn; % Mm 
   f(4)=1.0; % MSf 
   f(5)=1.043+0.414*cosn; % Mf 
   f(6)=sqrt((1.0+.203*cosn+.040*cos2n)^2+... 
        (.203*sinn+.040*sin2n)^2); % Mt 
   f(7)=1.0; % alpha1 
   f(8)=sqrt((1.+.188*cosn)^2+(.188*sinn)^2); % 2Q1 
   f(9)=f(8); % sigma1 
   f(10)=f(8); % q1 
   f(11)=f(8); % rho1 
   f(12)=sqrt((1.0+0.189*cosn-0.0058*cos2n)^2+...
         (0.189*sinn-0.0058*sin2n)^2); % O1 
   f(13)=1.0; % tau1 
   tmp1=1.36*cos(p*rad)+.267*cos((p-omega)*rad); % Ray's 
   tmp2= 0.64*sin(p*rad)+.135*sin((p-omega)*rad); 
   f(14)=sqrt(tmp1^2 + tmp2^2); % M1 
   f(15)=sqrt((1.+.221*cosn)^2+(.221*sinn)^2); % chi1 
   f(16)=1.0; % pi1 
   f(17)=1.0; % P1 
   f(18)=1.0; % S1 
   f(19)=sqrt((1.+.1158*cosn-.0029*cos2n)^2+...
         (.1554*sinn-.0029*sin2n)^2); % K1 
   f(20)=1.0; % psi1 
   f(21)=1.0; % phi1 
   f(22)=1.0; % theta1 
   f(23)=sqrt((1.+.169*cosn)^2+(.227*sinn)^2); % J1 
   f(24)=sqrt((1.0+0.640*cosn+0.134*cos2n)^2+...
         (0.640*sinn+0.134*sin2n)^2); % OO1 
   f(25)=sqrt((1.-.03731*cosn+.00052*cos2n)^2+...
         (.03731*sinn-.00052*sin2n)^2); % 2N2 
   f(26)=f(25); % mu2 
   f(27)=f(25); % N2 
   f(28)=f(25); % nu2 
   f(29)=1.0; % M2a 
   f(30)=f(25); % M2 
   f(31)=1.0; % M2b 
   f(32)=1.0; % lambda2 
   temp1=1.-0.25*cos(2.0*p*rad)...
        -0.11*cos((2.0*p-omega)*rad)-0.04*cosn; 
   temp2=0.25*sin(2.0*p)+0.11*sin((2.0*p-omega)*rad)...
        +0.04*sinn; 
   f(33)=sqrt(temp1^2+temp2^2); % L2 
   f(34)=1.0; % t2 
   f(35)=1.0; % S2 
   f(36)=1.0; % R2 
   f(37)=sqrt((1.+.2852*cosn+.0324*cos2n)^2+...
        (.3108*sinn+.0324*sin2n)^2); % K2 
   f(38)=sqrt((1.+.436*cosn)^2+(.436*sinn)^2); % eta2 
   f(39)=f(30)^2; % MNS2 
   f(40)=f(30); % 2SM2 
   f(41)=1.0; % M3 
   f(42)=f(19)*f(30); % MK3 
   f(43)=1.0; % S3 
   f(44)=f(30)^2; % MN4 
   f(45)=f(44); % M4 
   f(46)=f(44); % MS4 
   f(47)=f(30)*f(37); % MK4 
   f(48)=1.0; % S4 
   f(49)=1.0; % S5 
   f(50)=f(30)^3; % M6 
   f(51)=1.0; % S6 
   f(52)=1.0; % S7 
   f(53)=1.0; % S8 
   
   u(1)=0.0; % Sa 
   u(2)=0.0; % Ssa 
   u(3)=0.0; % Mm 
   u(4)=0.0; % MSf 
   u(5)=-23.7*sinn+2.7*sin2n-0.4*sin3n; % Mf 
   u(6)=atan(-(.203*sinn+.040*sin2n)/...
       (1.0+.203*cosn+.040*cos2n))/rad; % Mt 
   u(7)=0.0; % alpha1 
   u(8)=atan(.189*sinn/(1.+.189*cosn))/rad; % 2Q1 
   u(9)=u(8); % sigma1 
   u(10)=u(8); % q1 
   u(11)=u(8); % rho1 
   u(12)=10.8*sinn-1.3*sin2n+0.2*sin3n; % O1 
   u(13)=0.0; % tau1 
   u(14)=atan2(tmp2,tmp1)/rad; % M1 
   u(15)=atan(-.221*sinn/(1.+.221*cosn))/rad; % chi1 
   u(16)=0.0; % pi1 
   u(17)=0.0; % P1 
   u(18)=0.0; % S1 
   u(19)=atan((-.1554*sinn+.0029*sin2n)/...
        (1.+.1158*cosn-.0029*cos2n))/rad; % K1 
   u(20)=0.0; % psi1 
   u(21)=0.0; % phi1 
   u(22)=0.0; % theta1 
   u(23)=atan(-.227*sinn/(1.+.169*cosn))/rad; % J1 
   u(24)=atan(-(.640*sinn+.134*sin2n)/...
        (1.+.640*cosn+.134*cos2n))/rad; % OO1 
   u(25)=atan((-.03731*sinn+.00052*sin2n)/...
        (1.-.03731*cosn+.00052*cos2n))/rad; % 2N2 
   u(26)=u(25); % mu2 
   u(27)=u(25); % N2 
   u(28)=u(25); % nu2 
   u(29)=0.0; % M2a 
   u(30)=u(25); % M2 
   u(31)=0.0; % M2b 
   u(32)=0.0; % lambda2 
   u(33)=atan(-temp2/temp1)/rad; % L2 
   u(34)=0.0; % t2 
   u(35)=0.0; % S2 
   u(36)=0.0; % R2 
   u(37)=atan(-(.3108*sinn+.0324*sin2n)/...
        (1.+.2852*cosn+.0324*cos2n))/rad; % K2 
   u(38)=atan(-.436*sinn/(1.+.436*cosn))/rad; % eta2 
   u(39)=u(30)*2.0; % MNS2 
   u(40)=u(30); % 2SM2 
   u(41)=1.50*u(30); % M3 
   u(42)=u(30)+u(19); % MK3 
   u(43)=0.0; % S3 
   u(44)=u(30)*2.0; % MN4 
   u(45)=u(44); % M4 
   u(46)=u(30); % MS4 
   u(47)=u(30)+u(37); % MK4 
   u(48)=0.0; % S4 
   u(49)=0.0; % S5 
   u(50)=u(30)*3.0; % M6 
   u(51)=0.0; % S6 
   u(52)=0.0; % S7 
   u(53)=0.0; % S8 

   pu=zeros(1,ncmx); 
   pf=ones(1,ncmx); 

   for i=1:ncmx, 
      if(index(i)>0),
         pu(i)=u(index(i))*rad;
         pf(i)=f(index(i));
      end
   end

   % end of nodal

   time=(((k-1)/(24.0*60)+time_mjd)-48622.0)*86400;

   % make_a

   for j=1:ncmx,
      omegi(j)=omega_d(j);
      phase(j)=phase_mkB(j);
   end

   for j=1:nc,
      i=ind(j);
      if(i~=0),
         xc(j)=pf(i)*cos(omegi(i)*time+phase(i)+pu(i));
         yc(j)=pf(i)*sin(omegi(i)*time+phase(i)+pu(i));
      end
   end

   sum=0.0;
   for i=1:nc,
      sum=sum+xx(i)*xc(i)+yy(i)*yc(i);
   end
   zpred(k)=sum;
end

% gambarkan hasil prediksi
dateFormat='dd/mm/yyyy';
figure('Color','w');
plot(tgl,zpred);
grid on
set(gca, 'YLim', [-0.5, 0.5]);
set(gca, 'XLim', [tgl0 tgl0+15]);
title('Tidal Elevation Prediction');
ylabel('Elevation (m)');
datetick('x',dateFormat)