Thursday, August 18, 2005

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'}(:);
v_wind=nc{'vwnd',1}(1,:,:); % contoh hanya diambil jam pertama saja

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

% 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

%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';

% 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

% 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

lat_awal = 8.5713081359863; % sel baris ke-34
lat_akhir = -12.3807764053345; % sel baris ke-42
lon_awal = 93.75; % sel kolom ke-39
lon_akhir = 142.5; % sel kolom ke-58

% 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.

0 Comments:

Post a Comment

<< Home