Thursday, November 19, 2020

Plotting NPP hdf data

A sample code to plot the NPP data
The data source:
http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php 

 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

 

No comments: