Contoh prog ekstrak netcdf

18 08 2005

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.


Actions

Information

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s




%d bloggers like this: