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)

 

Menghitung salinitas dari data CTD

Beberapa minggu yang lalu seorang kawan, teman kuliah dulu yang sekarang sudah bergelar MEng. Oseanografi bertanya: “punyakah saya program untuk menghitung harga salinitas dari hasil pengamatan CTD? (jika ada dalam vb-nya excel, itu akan lebih baik). Saya yang memang tidak pernah bersentuhan dengan data CTD dengan singkat menjawab “tidak”, tapi saya ingat bahwa di websitenya Pak Tomzcak ada fasilitas itu dalam bentuk kalkulator online. Tapi ya tidak mangkus juga kalau data CTD yang banyak itu harus dimasukkan satu persatu ke kalkulator onlinenya Pak Tomzcak itu, bisa gempor dong tangan ini…

Setelah berlalu beberapa pekan, dan kebetulan punya waktu luang, saya mencoba mencari-cari di internet tentang cara menghitung salinitas dari data CTD, dan saya menemukan 2 alamat website ini: “Exercise 5: Using EXCEL to process CTD data” dan sebuah source code dalam Fortran untuk menghitung salinitas dan beberapa parameter fisis oseanografi lainnya dari CSIRO.

Karena saya terbiasa menggunakan Matlab untuk urusan hitung-menghitung dalam bidang oseanografi, akhirnya saya mencoba untuk mengubah kedua sumber berharga itu ke dalam fungsi-fungsi Matlab, dan terbentuklah dengan sukses function sal78, sal, dan conduc. Function sal78 adalah fungsi untuk mengubah rasio konduktivitas ke salinitas (jika M=0) dan mengubah salinitas ke rasio konduktivitas (jika M=1) berdasarkan pada UNESCO Report #37, 1981, juga ada dalam Practical Salinity Scale 1978 oleh E.L. Lewis, IEEE Ocean Eng. Jan, 1980. Prosedurnya mudah dan banyak tersedia rumusnya di dunia maya ini.

Untuk mengubah harga konduktivitas, temperatur dan tekanan dari data CTD, saya mencoba membuat fungsi sal, dimana dalam fungsi tersebut dipanggil sal78. Sedangkan untuk mengubah harga salinitas ke konduktivitas, saya buat fungsi conduc, dimana di dalamnya juga dipanggil fungsi sal78. Saya sudah coba bandingkan hasilnya dengan kalkulator onlinenya Pak Tomzcak, alhamdulillah harganya sama. Berikut adalah isi dari fungsi-fungsi tersebut,semoga bermanfaat:

1.Function sal78.m

function sal=sal78(c,t,p,m)
% notes from me:
% --------------
% some parts are taken from
% http://marine.csiro.au/datacentre/process/formats/sal78.f
% and the other parts are taken from:
% http://faculty.washington.edu/blewis/ocn499/EXER05.htm
%
% i just try to combine it for my needs and purposes
%
% this function is called by sal to get the salinity value from
% conductivity --> s=sal78(g,t,p,0)
%
% this function is called by conduc to get the conductivity value from
% salinity --> c=sal78(s,t,p,1)
%
% the sal and conduc codes are taken from:
% http://marine.csiro.au/datacentre/process/formats/sal78.f
%
% ***agus setiawan, sep.2006***
%
% this header and explanation is from:
% http://faculty.washington.edu/blewis/ocn499/EXER05.htm
%
% sal78 converts conductivity to salinity
% the conductivity ratio (c)=1.0000000 for salinity=35
% PSS-78. temperature=15.0 deg. celcius, and atmospheric
% pressure.
%
% references: also located in UNESCO Report NO. 37 1981
% Practical Salinity Scale 1978: E.L. Lewis IEEE Ocean Eng.
% Jan. 1980
%
% units:
% pressure p decibars
% temperature t deg. celcius (IPTS-68)
% conductivity c ratio (m=0)
% conductivity c mmho/sec (m=1)
% salinity sal78 (PSS-78) (m=0)
% checkvalues:
% sal78=1.888091:c=40.0000,t=40degC,p=10000dcbrs:m=1
% sal78=40.00000:c=1.888091,t=40degC,p=10000dcbrs:m=0
%
% sal78 ratio: returns zero for conductivity ratio: <0.0005
% sal78: returns zero for salinity: <0.02
%
% Practical Salinity Scale 1978 definition with temperature
% correction
%
% convert conductivity to salinity
%
% this part of code is from:
% http://faculty.washington.edu/blewis/ocn499/EXER05.htm
r=c;
dt=t-15.0;
rt35=(((1.0031E-9*t-6.9698e-7)*t+1.104259e-4)*t +2.00564e-2)*t+0.6766097;
c=((3.989e-15*p-6.370e-10)*p+2.070e-5)*p;
b=(4.464e-4*t+3.426e-2)*t+1.0;
a=-3.107e-3*t+0.4215;
if m==0,
rt=r/(rt35*(1.0+c/(b+a*r)));
rt=sqrt(abs(rt));
sal=((((2.7081*rt-7.0261)*rt+14.0941)*rt+25.3851)...
 *rt -0.1692)*rt+0.0080...
  +(dt/(1.0+0.0162*dt))*(((((-0.0144*rt+0.0636)*rt-0.0375)...
 *rt-0.0066)*rt-0.0056)*rt+0.0005);
else
% this part of code is from:
% http://marine.csiro.au/datacentre/process/formats/sal78.f
rt=sqrt(r/35.0);
for n=1:10,
 si=((((2.7081*rt-7.0261)*rt+14.0941)*rt+25.3851)*rt...
   -0.1692)*rt+0.0080...
  +(dt/(1.0+0.0162*dt))*(((((-0.0144*rt...
   +0.0636)*rt-0.0375)*rt-0.0066)*rt-0.0056)*rt+0.0005);
  dsal=((((13.5405*rt-28.1044)*rt+42.2823)*rt+50.7702)*rt...
   -0.1692)+(dt/(1.0+0.0162*rt))*((((-0.0720*rt+0.2544)*rt...
   -0.1125)*rt-0.0132)*rt-0.0056);
   rt=rt+(r-si)/dsal;
end
rtt=rt35*rt*rt;
c=rtt*(c+b);
b=b-rtt*a;
rr=sqrt(abs(b*b+4.0*a*c))-b;
sal=0.5*rr/a;
end  % if m==0

2.Function sal.m

function s=sal(c,t,p)
%
% function s=sal(c,t,p)
% this function derives salinity from a value of the
% in situ conductivity (as determined, for example, by a CTD),
% temperature and pressure.
%
% units:
% conductivity (c) : mmho/cm
% temperature (t)  : degC
% pressure (p)     : dbar --> 1 bar=10 dbar=100 kPa
%
%
% c35150 is the conductivity of standard seawater at 35 PSU,
% 15 degC, and atmospheric pressure
c35150=42.914;

if(abs(c)<0.01), s=-1, return, end
g=(c*(1.-6.5E-06*(t-2.8)+1.5E-8*(p-3000.)))/c35150;
s=sal78(g,t,p,0);

3.Function conduc.m

function c=conduc(s,t,p)
%
% function c=conduc(s,t,p)
% this function derives conductivity from a value of the
% salinity, temperature and pressure.
%
% units:
% salinity (s)     : -
% temperature (t)  : degC
% pressure (p)     : dbar  --> 1 bar = 10 dbar = 100 kPa
%
%
% c35150 is the conductivity of standard seawater at 35 PSU,
% 15 degC, and atmospheric pressure
c35150=42.914;
g=sal78(s,t,p,1);
c=g*c35150/(1.0-6.5e-6*(t-2.8)+1.5e-8*(p-3000.));

Contoh prog ekstrak netcdf

Pada tulisan terdahulu saya sudah memberikan cara bagaimana menginstal netcdf toolbox dan juga contoh untuk menggunakannya. Sekarang saya ingin memberikan sedikit uraian yang lebih detail tentang bagaimana mengekstrak data netcdf dari NCEP untuk wilayah Indonesia dan me-resample-nya.

Hal-hal yang perlu diingat jika kita mengekstrak data dari NCEP adalah tentang format gridnya. Ada data NCEP yang menggunakan format fixed grid, dan ada juga yang menggunakan format gaussian grid. Pada format fixed grid (contohnya data slp), jumlah selnya adalah 73 x 144 dengan resolusi 2.5 x 2.5 derajat. Sedangkan pada format gaussian grid (contohnya data angin), jumlah selnya adalah 94 x 192 (resolusi T62). Spasi sel arah zonal (bujur) adalah seragam sebesar 1.875 derajat, sedangkan spasi sel arah meridional (lintang) adalah tak seragam seperti di bawah ini:

88.5420 86.6532 84.7532 82.8508 80.9474
79.0435 77.1393 75.2351 73.3307 71.4262
69.5217 67.6171 65.7125 63.8079 61.9033
59.9986 58.0940 56.1893 54.2846 52.3799
50.4752 48.5705 46.6658 44.7611 42.8564
40.9517 39.0470 37.1422 35.2375 33.3328
31.4281 29.5234 27.6186 25.7139 23.8092
21.9044 19.9997 18.0950 16.1902 14.2855
12.3808 10.4760  8.5713  6.6666  4.7618
2.8571  0.9524

untuk BBU, dan dengan angka yang sama tapi dengan tanda minus untuk BBS.

Jika kita menggunakan ncdump -h pada file netcdf, hal ini bisa dilihat pada bagian awal header. Jika data berformat fixed grid, maka di bagian awal header akan tertulis:

dimensions:
 lon = 144 ;
 lat = 73 ;

Sedangkan jika data berformat gaussian, maka akan tertulis:

dimensions:
 lon = 192 ;
 lat = 94 ;

Berikut adalah contoh programnya:

% contoh ekstraksi data netcdf dari ncep utk wilayah indonesia
% dan bagaimana cara meresample sesuai resolusi yang kita inginkan
%
% program oleh agus setiawan - 2005
% silahkan digunakan sesuka hati asal jangan diaku milik sendiri

clear all

% baca file ncep
nc=netcdf('vwnd.10m.gauss.1997.nc','nowrite');
ygrid=nc{'lat'}(:);
xgrid=nc{'lon'}(:);
scalef=nc{'vwnd'}.scale_factor(:);
addoff=nc{'vwnd'}.add_offset(:);
v_wind=nc{'vwnd',1}(1,:,:);  % contoh hanya diambil jam pertama saja
v_wind=v_wind.*scalef+addoff;

nc=netcdf('uwnd.10m.gauss.1997.nc','nowrite');
u_wind=nc{'uwnd',1}(1,:,:);  % contoh hanya diambil jam pertama saja
u_wind=u_wind.*scalef+addoff;

% perhatikan baik-baik!
% ---------------------
% data ncep terdiri dari 2 format grid
% fixed (2.5 x 2.5 derajat dengan jumlah sel 73 x 144)
% dan gaussian (T62 dengan jumlah sel 94 x 192)
%
% rubah di sini utk menentukan mau ngolah format yang mana
%
% set sel='fixed' atau sel='gauss'
sel='gauss';

if(sel=='fixed'),

dx=2.5;  % selang sel dalam derajat

% definisi untuk daerah indonesia
lat_awal  =   7.5;   % sel baris ke-34
lat_akhir = -12.5;   % sel baris ke-42
lon_awal  =  95.0;   % sel kolom ke-39
lon_akhir = 142.5;   % sel kolom ke-58

%definisikan sel dalam bentuk lintang-bujur untuk keperluan resampling
[lat,lon]=meshgrid(lat_awal:-dx:lat_akhir,lon_awal:dx:lon_akhir);
lat=lat'; lon=lon';

% indeks untuk ekstraksi
% nilai ini bisa dicari dengan cara:
% for i=1:144, lon(i)=(i-1)*360/144;
% for i=1:73, lat(i)=(i-1)*180/73;
baris_awal = 34;
baris_akhir= 42;
kolom_awal = 39;
kolom_akhir= 58;

% ekstrak daerah indonesia
u_wind_ind=u_wind(baris_awal:baris_akhir,kolom_awal:kolom_akhir);
v_wind_ind=v_wind(baris_awal:baris_akhir,kolom_awal:kolom_akhir);

elseif(sel=='gauss'),

lat_temp=[8.5713081359863   6.6665730476379  4.7618379592896  2.8571028709412  0.9523676037788...
         -0.9523676037788  -2.8571028709412 -4.7618379592896 -6.6665730476379 -8.5713081359863...
        -10.4760417938232 -12.3807764053345];
lat_temp=lat_temp';
lon_temp=93.75:1.875:142.5;
lat1=ones(size(lat_temp,1),size(lon_temp,2));
for i=1:size(lat1,2), lat(:,i)=lat_temp(:)+lat1(:,i); end
lon1=ones(size(lat_temp,1),size(lon_temp,2));
for i=1:size(lon1,2), lon(:,i)=lon_temp(i)+lon1(:,i); end

% indeks untuk ekstraksi
% nilai ini bisa dicari dengan cara:
% for i=1:192, lon(i)=(i-1)*360/192;
% [xlat,dlat,sinc]=gauss2lats(nlat);   % perlu fungsi eksternal gauss2lats
baris_awal = 43;
baris_akhir= 54;
kolom_awal = 51;
kolom_akhir= 77;

% ekstrak daerah indonesia
u_wind_ind=u_wind(baris_awal:baris_akhir,kolom_awal:kolom_akhir);
v_wind_ind=v_wind(baris_awal:baris_akhir,kolom_awal:kolom_akhir);

end

% definisikan parameter resampling
dx_r=0.5;
lat_r_awal  =   6.0;  % batas utara
lat_r_akhir = -11.0;  % batas selatan
lon_r_awal  =  95.0;  % batas barat
lon_r_akhir = 141.0;  % batas timur

[lat_r,lon_r]=meshgrid(lat_r_awal:-dx_r:lat_r_akhir,lon_r_awal:dx_r:lon_r_akhir);
lat_r=lat_r'; lon_r=lon_r';

% resampling dan interpolasing
% ----------------------------
% catatan:
% tersedia 4 metode interpolasi:
% 'nearest' - nearest neighbor interpolation
% 'linear'  - bilinear interpolation
% 'cubic'   - bicubic interpolation
% 'spline'  - spline interpolation
u_wind_ind_r=interp2(lon,lat,u_wind_ind,lon_r,lat_r,'spline');
v_wind_ind_r=interp2(lon,lat,v_wind_ind,lon_r,lat_r,'spline');

% cek : plot vektor angin
quiver(lon_r,lat_r,u_wind_ind_r,v_wind_ind_r)

Selamat mencoba dan mengembangkannya.

Baca unformatted dgn Matlab

Misalkan kita memiliki program dalam fortran untuk menyimpan file dalam format unformatted sebagai berikut:

integer n,m,i,j
parameter (n=5,m=10)
real a(n,m),c,d

c… proses perhitungan anda…
c… (terserah seperti apa)…

c… misalnya hasil disimpan dalam format seperti di bawah ini…
open(11,file=’tes.bin’,form=’unformatted’,status=’unknown’)
write(11) c,d
write(11) a
close(11)

stop
end

Dan kita ingin mengolahnya lebih lanjut dengan menggunakan Matlab (misalnya untuk menampilkan gambarnya). Maka kita dapat membaca hasil keluaran dari fortran tersebut dengan perintah sebagai berikut:

fid=fopen(‘tes.bin’,’rb’);
fseek(fid,4,’bof’);
c=fread(fid,4,’float’)
d=fread(fid,4,’float’);
fseek(fid,8,’cof’);
a=fread(fid,[5,10],’float’);

Script Matlab background

Bagaimana menjalankan program yang anda buat dengan menggunakan matlab secara background (biasanya dipakai untuk program dengan waktu running yang cukup lama).

pada sistem operasi linux:

buat sebuah script (misalnya dengan nama file ‘matbg’) sbb:

#!/bin/csh -f

# Clear the DISPLAY.
unsetenv DISPLAY

# Call MATLAB with the appropriate input and output,
# make it immune to hangups and quits using ”nohup”,
# and run it in the background.
nohup matlab $2 &

cara menggunakan:

dari direktori tempat program itu berada jalankan:

matbg

pada sistem operasi windows:

buat sebuah batch file (misalnya dengan nama file ‘matbg.bat’) sbb:

matlab -nosplash -nodesktop -minimize -r nama_program -logfile log_file

(jangan lupa ganti ‘nama_program’ dengan nama program matlab yang mau anda jalankan dan ‘log_file’ dengan nama logfile yang anda kehendaki).

jalankan program dengan men-double click ‘matbg.bat’ atau dapat juga ketikan ‘matbg’ melalui dos prompt.

Baca file ASCII dgn Matlab

Pada tulisan kali ini saya akan mencoba membahas tentang membaca file dalam format ASCII dalam berbagai bentuk yang umumnya ada dalam bidang oseanografi.

Jika jumlah kolom pada setiap baris dari file ASCII yang hendak dibaca seragam (jadi membentuk sebuah matriks m baris x n kolom), maka kita dapat memanggilnya dengan cara berikut:

misalnya kita memiliki data kecepatan dan arah arus laut ‘arus.dat’ :

1.0 180.0
1.5 220.0
2.0 270.0

cara pertama:

arus=load(‘arus.dat’);
kec = arus(:,1);
arah = arus(:,2);

cara kedua:

[kec arah] = textread(‘arus.dat’,’%f %f’)

Namun kadangkala juga kita mempunyai file data dengan kategori di atas namun pada baris pertama terdapat header seperti (misalnya nama file data adalah ‘arus.dat’):

kec arah
1.0 180.0
1.5 220.0
2.0 270.0

untuk membacanya, berikan perintah:

[kec arah] = textread(‘arus.dat’,’%f %f’,’headerlines’,1)

Jika pada file yang kita miliki terdapat jumlah kolom yang tak seragam, misalnya kita memiliki data pengukuran elevasi pasang surut setiap jam ‘pasut.dat’ dengan format:

1.0 1.1 1.2
1.3 1.4 1.5 1.6
1.7 1.9

maka untuk membacanya dapat digunakan perintah:

elevasi = textread(‘pasut.dat’, ‘%f’)

atau jika pada baris pertama terdapat header seperti di bawah ini:

data pasut Riau
1.0 1.1 1.2
1.3 1.4 1.5 1.6
1.7 1.9

maka dapat dugunakan perintah:

elevasi = textread(‘pasut.dat’, ‘%f’, ‘headerlines’,1 )

Jika kita memiliki data ‘arus1.dat’ yang merupakan data pengamatan arus dengan selang pengukuran 1 jam yang dipisahkan oleh tanda koma (format csv, biasanya hasil upload dari alat pengukuran sepertu current meter menggunakan format ini) seperti:

data arus Riau
dd/mm/yyyy,hh:mm:ss,kec,arah
1/1/2000,00:00:00,1.0,180.0
1/1/2000,01:00:00,1.5,220.0
1/1/2000,02:00:00,2.0,270.0

maka kita dapat membaca harga kecepatan dan arah arus dengan perintah:

arus = csvread(‘arus1.dat’, 2, 2);

dimana angka 2 yang pertama pada perintah di atas berarti kita mulai pembacaan data setelah 2 baris pertama diabaikan, dan angka 2 yang kedua berarti kita mulai pembacaan data setelah 2 kolom pertama diabaikan.

selanjutnya kita dapat mendefiniskan kecepatan dan arah arus dengan perintah:

kec = arus(:,1);
arah = arus(:,2);

Selamat mencoba! Dalam kesempatan lain, jika ada cukup waktu, saya akan mencoba untuk menulis artikel tentang bagaimana membaca file data dalam format biner/unformatted.