Visualize Sea Ice Concentrations from GRIB Data

Since R2023b

Visualize sea ice concentrations for an arctic region using data stored in the GRIB file format. GRIB files typically store meteorological, hydrological, or oceanographic data. In many cases, GRIB files store data using multiple bands. This example shows you how to:

• View sea ice concentrations for a single year.

• Animate sea ice concentrations over time.

• Compare sea ice concentrations for multiple years using histograms.

This example uses data from a GRIB file containing sea ice concentrations for 2010, 2012, 2014, 2016, 2018, 2020, and 2022 [1][2]. The file stores the data as a grid of posting point samples referenced to geographic coordinates.

Visualize Concentrations for Single Year

Visualize sea ice concentrations for 2010 by using an axesm-based map.

Get Information About GRIB File

Get information about the GRIB file by creating a RasterInfo object. Extract the latitude limits, longitude limits, reference times, and number of bands from the object. Each data band corresponds to a year.

info = georasterinfo("seaice.grib");
latlim = info.RasterReference.LatitudeLimits;
lonlim = info.RasterReference.LongitudeLimits;
refTime = info.Metadata.ReferenceTime;
numBands = info.NumBands;

Read Data Band

Read the data band for 2010 by using the metadata stored in the RasterInfo object. RasterInfo objects store GRIB metadata using a table, where each row of the table corresponds to a data band.

Find the band that contains data for 2010.

yrs = year(info.Metadata.ReferenceTime);
band2010 = find(yrs == 2010);

Extract the band for 2010. Replace missing data with NaN values by using the standardizeMissing function. Then, convert the concentrations to percentages.

[iceConcentration2010,R2010] = readgeoraster("seaice.grib",Bands=band2010);

m = info.MissingDataIndicator;
iceConcentration2010 = standardizeMissing(iceConcentration2010,m);

percentIce2010 = 100*iceConcentration2010;

Get the reference time for the band.

refTime2010 = refTime(band2010);
refTime2010 = string(refTime2010,"MMMM d, yyyy");

Display Data Band

Load an ice colormap that is appropriate for ice concentration data. The colormap starts at dark blue and transitions to light blue and white.

load iceColormap.mat

Specify ice concentration levels from 0% to 100%. The level increases every 5%.

levels = 0:5:100;

Set up an axesm-based map for an arctic region by using the setUpArcticAxesm helper function. The function creates a map in a polar stereographic projection, zooms into the arctic region, sets the colormap, and adds a labeled color bar. Specify the colormap as the ice colormap. Specify the tick labels for the color bar using the ice concentration levels.

figure
setUpArcticAxesm(iceColormap,levels)

Provide geographic context for the map by displaying world land areas in gray. Then, display the sea ice concentrations as a surface. To make the surface appear under the land areas, display a matrix of zeros and apply color using the ice concentrations.

land = readgeotable("landareas.shp");
landColor = [0.8 0.8 0.8];
geoshow(land,FaceColor=landColor)

geoshow(zeros(R2010.RasterSize),R2010,DisplayType="surface",CData=percentIce2010)

Add a title and subtitle.

title("Sea Ice Concentrations")
subtitle(refTime2010)

Animate Concentrations Over Time

Visualize sea ice concentrations over time by using an animated GIF.

Read Data for All Bands

Read all the data bands into the workspace. Replace missing data with NaN values by using the standardizeMissing function. Then, convert the concentrations to percentages.

[iceConcentration,R] = readgeoraster("seaice.grib");

m = info.MissingDataIndicator;
iceConcentration = standardizeMissing(iceConcentration,m);

percentIce = 100 * iceConcentration;

Create Animated GIF

Create an animated GIF and save it to a file.

Set up a new axesm-based map by using the setUpArcticAxesm helper function. Avoid displaying the animation by using an invisible figure. Add a title.

f = figure(Color="w",Visible="off");
setUpArcticAxesm(iceColormap,levels)
title("Animation of Sea Ice Concentration")

Display the ice concentrations for 2010 as a surface. To use the same figure for each frame of the GIF, return the surface as a variable. This strategy enables you to display ice concentration data for different years by changing the color data for the surface. Then, provide geographic context for the map by displaying the land areas.

iceSurface = geoshow(zeros(R.RasterSize),R2010,DisplayType="surface",CData=percentIce2010);
geoshow(land,FaceColor=landColor)

Specify a name for the GIF. If a file with that name already exists, delete it.

gifFilename = "seaice_animation.gif";
if exist(gifFilename,"file")
    delete(gifFilename)
end

Create an animation. For each year of data:

• Extract the data band for the year. Change the color data for the surface using the data band.

• Change the subtitle of the map to match the year.

• Write the figure to the GIF by using the writeAnimatedGIF helper function.

for yr = 1:numBands
    % extract data band
    percentIceYr = percentIce(:,:,yr);
    iceSurface.CData = percentIceYr;

    % change subtitle
    refTimeYr = string(refTime(yr),"MMMM d, yyyy");
    subtitle(refTimeYr)
    
    % append map to GIF file
    drawnow
    writeAnimatedGIF(f,"seaice_animation.gif")
end

This image shows the animated GIF.

Compare Years Using Histograms

Quantitatively visualize the data using histograms. Compare the data for two years by displaying two histograms in the same figure. Avoid analyzing ice concentrations below a specified level using a threshold.

Select two years to compare.

yr1 = 2010;  
yr2 = 2022;

For each year, find the data band, read the data band, replace missing data with NaN values, convert the concentrations to percentages, and get the reference time.

bandYr1 = find(yrs == yr1);
bandYr2 = find(yrs == yr2);

[iceConcentrationYr1,RYr1] = readgeoraster("seaice.grib",Bands=bandYr1);
[iceConcentrationYr2,RYr2] = readgeoraster("seaice.grib",Bands=bandYr2);

iceConcentrationYr1 = standardizeMissing(iceConcentrationYr1,m);
iceConcentrationYr2 = standardizeMissing(iceConcentrationYr2,m);

percentIceYr1 = 100*iceConcentrationYr1;
percentIceYr2 = 100*iceConcentrationYr2;

refTimeYr1 = string(refTime(bandYr1),"MMMM d, yyyy");
refTimeYr2 = string(refTime(bandYr2),"MMMM d, yyyy");

Select a threshold for the ice concentrations.

threshold = 30;

For each year, extract the data points that are above the threshold. Find the percentage of data points that are above the threshold, excluding the NaN values.

% yr1
idxAboveThresholdYr1 = percentIceYr1 > threshold;
aboveThresholdYr1 = percentIceYr1(idxAboveThresholdYr1);
percentAboveThresholdYr1 = 100*numel(aboveThresholdYr1)/sum(~isnan(percentIceYr1),"all");

% yr2
idxAboveThresholdYr2 = percentIceYr2 > threshold;
aboveThresholdYr2 = percentIceYr2(idxAboveThresholdYr2);
percentAboveThresholdYr2 = 100*numel(aboveThresholdYr2)/sum(~isnan(percentIceYr2),"all");

Display the ice concentrations using two histograms. Specify the bin edges using the levels. Display the threshold using a vertical line.

figure
histogram(aboveThresholdYr1,levels,FaceAlpha=0.5)
hold on
histogram(aboveThresholdYr2,levels,FaceAlpha=0.5)
xline(threshold,"-",{'Threshold '},HandleVisibility="off")

Display the year and the percentage of data points that are above the threshold in a legend.

labely1 = sprintf("%d \nPoints above threshold: %.2f%%",yr1,percentAboveThresholdYr1);
labely2 = sprintf("%d \nPoints above threshold: %.2f%%",yr2,percentAboveThresholdYr2);
legend(labely1,labely2,Location="southoutside")

Add titles and labels.

title("Comparison of Ice Concentration Levels")
subtitle(refTimeYr1 + " and " + refTimeYr2)
xlabel("Sea Ice Concentration (%)")
ylabel("Number of Points")

When data is stored as a grid of rectangular cells referenced to geographic coordinates, you can find the surface area by using the areamat function. However, most GRIB data is stored as grids of posting point samples. For more information about the difference between cells and posting points, see Spatially Reference Imported Rasters.

References

[1] Hersbach, H., B. Bell, P. Berrisford, G. Biavati, A. Horányi, J. Muñoz Sabater, J. Nicolas, et al. "ERA5 Hourly Data on Single Levels from 1940 to Present." Copernicus Climate Change Service (C3S) Climate Data Store (CDS), 2023. Accessed May 22, 2023. https://doi.org/10.24381/cds.adbb2d47.

[2] Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

Helper Functions

The setupArcticAxesm function sets up an axesm-based map for an arctic region using the colormap cmapIn and the ice concentration levels levelsIn. Within the function:

• Create an axesm-based map that uses a polar stereographic projection. Specify the origin using the coordinates of the North Pole.

• Zoom in to an arctic region by changing the xy-limits of the map. Specifying Cartesian limits instead of geographic limits enables you to display more data within the plot box.

• Set the colormap. Set the color limits using the minimum and maximum values of the ice concentration levels.

• Add a labeled color bar. Specify the ticks using the ice concentration levels.

function setUpArcticAxesm(cmapIn,levelsIn)
    % create axesm-based map
    abm = axesm("stereo",Origin=[90 0 0],Grid="on",GColor="w");

    % zoom in to arctic region
    nlim = 0.45;
    set(abm,XLim=[-nlim nlim],YLim=[-nlim nlim])

    % set colormap and color limits
    colormap(cmapIn)
    clim(abm,[min(levelsIn) max(levelsIn)])

    % add labeled color bar
    cb = colorbar;
    cb.Ticks = levelsIn;
    cb.Label.String = "Sea Ice Concentration (%)";
end

The writeAnimatedGIF function writes the figure figIn to the animated GIF with name filename. Within the function:

• Extract the color data from the figure. Convert the RGB image to an indexed image. Preserve the spatial resolution by turning off dithering.

• Write the image to the GIF. If the GIF does not exist, create a new GIF using the image. If the GIF does exist, append the image to the GIF.

function writeAnimatedGIF(figIn,filename)
    % create indexed image from figure
    data = getframe(figIn);
    RGB = data.cdata;
    [figFrame,figCmap] = rgb2ind(RGB,256,"nodither");
    
    % write image to GIF
    if ~exist(filename,"file")
        imwrite(figFrame,figCmap,filename,"gif",Loopcount=Inf,DelayTime=0.5)
    else
        imwrite(figFrame,figCmap,filename,"gif",WriteMode="append",DelayTime=0.5)
    end
end