filewant = 'cafe.2007001.hdf';
info = hdfinfo(filewant);
% This is how you can investigate the file structure...
info.SDS
%% note the names
info.SDS.Attributes.Name
%% note the values of the attributes
info.SDS.Attributes.Value
%% read the npp variable
npp = double(hdfread(filewant, 'npp' ));
npp(npp==-9999) = NaN;
% figure;
% imagesc(npp)
npp = flipud(npp); % because npp is flipped upside down wrt lat/lon for map
%%
iLat = linspace(-90,90,1080);
iLon = linspace(-180,180,2160);
[Lat, Lon] = ndgrid(iLat,iLon);
figure
load coastlines % coastineData.mat
latlim=[floor(min(min(Lat))),ceil(max(max(Lat))) ];
lonlim=[floor(min(min(Lon))),ceil(max(max(Lon))) ];
ax = worldmap(latlim, lonlim);
surfacem(Lat, Lon, npp);
geoshow(coastlat, coastlon,'Color', 'k' )
colormap; set(gcf,'Color','white')
colormap jet
map2 = colormap;
map2( 1, : ) = 1; colormap(map2);
colorbar
title('global dataset')
disp('Next we will try to crop out the area of interest')
%% This is an example of how to crop out to the area of interest...
iLat2 = linspace(-15,45,1080);
iLon2 = linspace(-90,0,2160);
[Lat2, Lon2] = ndgrid(iLat2,iLon2);
F = griddedInterpolant(Lat,Lon,npp); % read doc griddedInterpolant
npp_crop = F(Lat2,Lon2);
%
figure
load coastlines % coastineData.mat
latlim=[floor(min(min(Lat2))),ceil(max(max(Lat2))) ];
lonlim=[floor(min(min(Lon2))),ceil(max(max(Lon2))) ];
ax = worldmap(latlim, lonlim);
surfacem(Lat2, Lon2, npp_crop);
geoshow(coastlat, coastlon,'Color', 'k' )
colormap; set(gcf,'Color','white')
colormap jet
map2 = colormap;
map2( 1, : ) = 1; colormap(map2);
colorbar
title('NPP (mgC m-2 day-1)')
Alternative code that I have seen:
https://hdfeos.org/zoo/OTHER/npp_2010361_hdf.m
