% Test writing and reading NetCDF4 with chunking and compression using
% Native Matlab routines (tested in Matlab 2010b)
% Rich Signell (rsignell@usgs.gov)

filename='netcdf4_test.nc'
if 0
    % read data from OpenDAP using NJ Toolbox (http://njtbx.sourceforge.net)
    url='http://coast-enviro.er.usgs.gov/thredds/dodsC/bathy/etopo1_bed_g2';
    nc=mDataset(url);
    topo=nc{'topo'}(1:12:end,1:12:end);
    g=nc{'topo'}(1:12:end,1:12:end).grid;
    topo(topo<0)=0;
    lon=g.lon;
    lat=g.lat;
    save topo.mat topo lon lat
end
% or load previously save mat file
load topo.mat
[ny,nx]=size(topo);
subplot(211);pcolor(lon,lat,double(topo));shading flat;caxis([0 5000])
title('Topo from Mat file');
%%
% write NetCDF4 with chunking & compression (deflation)

ncid = netcdf.create(filename,'NETCDF4');
latdimid = netcdf.defDim(ncid,'lat',ny);
londimid = netcdf.defDim(ncid,'lon',nx);
varid = netcdf.defVar(ncid,'topo','short',[latdimid londimid]);
lonid = netcdf.defVar(ncid,'lon','float',[londimid]);
latid = netcdf.defVar(ncid,'lat','float',[latdimid]);

netcdf.defVarChunking(ncid,varid,'CHUNKED',[180 360]);
netcdf.defVarDeflate(ncid,varid,true,true,5);
netcdf.putAtt(ncid,latid,'units','degrees_north');
netcdf.putAtt(ncid,lonid,'units','degrees_east');
netcdf.putAtt(ncid,varid,'units','m');
%netcdf.putAtt(ncid,varid,'missing_value',int16(-32767));

netcdf.putVar(ncid,lonid,[0],[nx],lon(1:nx));
netcdf.putVar(ncid,latid,[0],[ny],lat(1:ny));
netcdf.putVar(ncid,varid,[0 0],[ny nx],topo(1:ny,1:nx));

netcdf.close(ncid);

% read NetCDF4 file
ncid=netcdf.open(filename,'nowrite');
varid=netcdf.inqVarID(ncid,'topo');
topo2=netcdf.getVar(ncid,varid,[0 0],[ny nx]);
netcdf.close(ncid);
subplot(212);pcolor(lon,lat,double(topo2));shading flat;caxis([0 5000])
title('Topo from NetCDF4 file');
%%