Plotting global flux data with colormap

I have global flux data and latitude longitude data with dimension: flux=9*30*4, lat=9*30*4, lon=9*30*4. How can I get a colormap that will show the value (color) of the flux data on a global map including the color axis next to it? I have the mapping toolbox. I am looking to plot image like the one below. I tried surfacem, but my matlab crashes as the data is too big to load.

Answer 1

[Kelly Kearney]

Your example image looks very Ferret-y, so I'm assuming your data is currently in netCDF files. Can you give us some more details about those files? Are you sure about those variable sizes? They seem awfully small for what looks like satellite data (and a handful of 9 x 30 x 4 matrices are nowhere near large enough to cause Matlab memory problems).

Can you show us some code you've tried, and the full text of the error message you receive when you try surfacem?

[Chowdhury Nazmi]

Hi Kelly,

Thanks for responding. Yes these are netCDF files. Actually the dimension is only for one file. There are about 500 .nc files that I have to read through to get the whole globe's data. The files I am using are from the link below: <ftp://ftp-npp.class.ngdc.noaa.gov/20160428/NDE-L2/NUCAPS-Outgoing-Longwave-Radiation/> and I am using the first file : NDE-L2_NUCAPS-Outgoing-Longwave-Radiation_20160428_06475.tar and then when I extracted it , there are about 500 files. The code is attached below.

clear all
close all
nucaps_olr_directory_path=['C:\Users\nazmi.chowdhury\Documents\MATLAB\NDE-L2_NUCAPS-Outgoing-Longwave-Radiation_20160428_06475\OLR_NCF'];
nucaps_olr_directory_info=dir(nucaps_olr_directory_path);
lat_store=[];
lon_store=[];
flux_store=[];
figure
latlim=[-90 90];
lonlim=[-180 180]; 
ax=worldmap(latlim,lonlim);
oceanColor = [.5 .7 .9];
setm(ax, 'FFaceColor', oceanColor)
land = shaperead('landareas', 'UseGeoCoords', true);
geoshow(ax, land, 'FaceColor', [0.5 0.7 0.5])
lakes = shaperead('worldlakes', 'UseGeoCoords', true);
geoshow(lakes, 'FaceColor', 'blue')
hold on
for i=3:length(nucaps_olr_directory_info)
i
file_name = nucaps_olr_directory_info(i).name; 
file_full_name = [nucaps_olr_directory_path file_name]; 
lat=ncread(file_name,'LAT');
lat=reshape(lat,1080,1);
lon=ncread(file_name,'LON');
lon=reshape(lon,1080,1);
flux=ncread(file_name,'FLUX');
flux=reshape(flux,1080,1);
lat_store=[lat_store;lat];
lon_store=[lon_store;lon];
flux_store=[flux_store;flux];    
end
surfacem(lat_store,lon_store,flux_store); hold on;
h2=colorbar('vert');

The error that I am getting is : Error using repmat Requested 540000x540000 (2172.6GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

Error in meshgrid (line 58) xx = repmat(xrow,size(ycol));

Error in meshgrat (line 123) [lon,lat] = meshgrid(lon,lat);

Error in surfacem>parseInputs (line 110) [lat,lon] = meshgrat(lat,lon);

Error in surfacem (line 70) [lat, lon, alt, grid, pvPairs] = parseInputs(varargin{:});

Error in CrIS_OLR (line 43) surfacem(lat_store,lon_store,flux_store); hold on;

[Kelly Kearney]

The surfacem routine is designed for plotting gridded data. Your data isn't gridded; it's scattered. (Ideally, the surfacem function would throw an error or warning before it tried to plot this... as it is, it saw that you passed vectors as the first two inputs and assumed that you wanted a 540000 x 540000 grid, which would of course be far too large a dataset.)

In order to plot this properly, you have two choices. You can either plot the raw data as a scatter plot, or regrid the data. The scatter option is probably a bad idea here... 540000 points is definitely overkill for plotting purposes (and Matlab's scatter function is one of its more inefficient plotting routines, even after some recent improvements). So I recommend regridding (see bottom example):

% Read data

folder = '~/Downloads/NDE-L2_NUCAPS-Outgoing-Longwave-Radiation_20160428_06475';
Ncfiles = dir(fullfile(folder, '*.nc'));
nfile = length(Ncfiles);

[lat,lon,flux] = deal(nan(1080,nfile));

for ii = 1:length(Ncfiles)
    lat(:,ii)  = reshape(ncread(fullfile(folder, Ncfiles(ii).name), 'LAT'), [], 1);
    lon(:,ii)  = reshape(ncread(fullfile(folder, Ncfiles(ii).name), 'LON'), [], 1);
    flux(:,ii) = reshape(ncread(fullfile(folder, Ncfiles(ii).name), 'FLUX'), [], 1);
end

% Regrid

latlim=[-90 90];
lonlim=[-180 180]; 

F = scatteredInterpolant(lon(:), lat(:), flux(:), 'linear', 'none');
latg = linspace(latlim(1), latlim(2), 180);
long = linspace(lonlim(1), lonlim(2), 360);
fluxg = F({long, latg});

% Plot scatter

oceanColor = [.5 .7 .9];
land = shaperead('landareas', 'UseGeoCoords', true);
lakes = shaperead('worldlakes', 'UseGeoCoords', true);

figure('color', 'w');
subplot(2,1,1);
ax(1) = worldmap(latlim,lonlim);
geoshow(land, 'FaceColor', [0.5 0.7 0.5])
geoshow(lakes, 'FaceColor', 'blue');
scatterm(lat(:), lon(:), 2, flux(:), 'filled');
colormap(jet);
title('scatter');

% Plot surface

subplot(2,1,2);
ax(2) = worldmap(latlim,lonlim);
geoshow(land, 'FaceColor', [0.5 0.7 0.5])
geoshow(lakes, 'FaceColor', 'blue');
surfacem(latg, long, fluxg');
colormap(jet);
title('regridded');

(I only unzipped 5 files...)

[Chowdhury Nazmi]

Hi Kelly, Thank you very much. This worked. I really appreciate your help.

[Chowdhury Nazmi]

Hi Kelly, Just plotted the whole globe using all the files, there are some NaN values in the lat/lon data, which is causing the scatteredInterpolant() to fail. It gives error: Error using scatteredInterpolant The coordinates of the input points must be finite values; Inf and NaN are not permitted.

Error in Untitled2 (line 17) F = scatteredInterpolant(lon(:), lat(:), flux(:), 'linear', 'none');

[Kelly Kearney]

This means you need to strip out any invalid data point before building the interpolant:

isn = isnan(lon) | isnan(lat) | isnan(flux);
F = scatteredInterpolant(lon(~isn), lat(~isn), flux(~isn), 'linear', 'none');

[Silver]

Hi Kelly !

I have an "out of memory" error while reading chlorophyl data from an netcdf file. The file contains :

• latitude
• longitude
• time
• CHL ( chl_a)

when I use ncread for latitude longitude and time , it works without any issue , but when I use the ame function for chlorophyl read, it shows the out of memory error.

I used this code :

ncdisp('dataset-oc-med-chl-multi-l4-interp_1km_daily-rt-v02_1541025401266.nc'); % to display the nc file
lon2 = ncread('dataset-oc-med-chl-multi-l4-interp_1km_daily-rt-v02_1541025401266.nc','lon') ; %to read longitude
lat2 = ncread('dataset-oc-med-chl-multi-l4-interp_1km_daily-rt-v02_1541025401266.nc','lat') ; %to read latitude
time2 = ncread('dataset-oc-med-chl-multi-l4-interp_1km_daily-rt-v02_1541025401266.nc','time') ; %to read time
chl = ncread('dataset-oc-med-chl-multi-l4-interp_1km_daily-rt-v02_1541025401266.nc','CHL') ; %to read chlorophyl , the error shows in this line

the file size is 955 Mo.

any ideas ? Thks in advance!

[Chad Greene]

Silver: You posed this question both as a comment on an answer to an old question, and as an answer to the old question. I deleted the answer version, and I'll offer some guidance here, but you'll generally have better luck starting a new question whenever you have a new question.

If you're running out of memory I recommend just loading a small portion of the data when you call ncread. To do that, use the start and count options. If you're only interested in a few years of data, just load those years. If you only need a small region of the world, just load that region's worth of data. If you need all the data from all times and all locations, you might need to find a way to do your processing in chunks.