Finger-rafted sea ice, seen from the helicopter.

Photo: Nick Cobbing

function [lats,lons,C_mon,e_lat,e_lon] = SSMI_Arctic_download_process_plot(months_to_plot,~)

% -------------------------------------------------------------------------

% This script downloads, processes, and plots monthly sea ice concentration 

% from Nimbus-7 SMMR and DMSP SSM/I-SSMIS which was pre-processed using the

% Nasa-Team algorithm.

%

% Documentation for this product is found here:

% http://nsidc.org/data/NSIDC-0051

%

% There is 1 necessary input argument and 1 optional argument.

%

% The first argument specifies the months you want to plot. If you want 

% plots of the March Maximum and the September Minimum, run

%

% [lats,lons,C_mon,e_lat,e_lon] = SSMI_Arctic_download_process_plot([3 9])

%

% The second, optional, argument should ONLY be set the first time you run

% the script, since it will download the data from the NOAA ftp server. 

% It's about 50MB so it shouldn't take long. For this, run (e.g.)

%

% [lats,lons,C_mon,e_lat,e_lon] = SSMI_Arctic_download_process_plot([3 9],1) 

%

% This also downloads the lat/lon grid files that you'll need.

%

% The outputs give the following:

%       lats: latitudes

%       lons: longitudes

%       C_mon: a 12x1 cell array with climatological mean sea ice

%              concentrations for each month

%       e_lat: a 12x1 cell array with the latitude values of the ice edges 

%              corresponding to C_mon

%       e_lon: a 12x1 cell array with the longitude values of the ice edges 

%              corresponding to C_mon

%

% The plotted figures give the climatological mean sea ice concentration

% over the period 1979-2015 (years can easily be changed in the code).

% They further show (in red) the (marine) ice edge. 

% As per convention, the ice edge is set to be where the concentration 

% is <15%. 

%

% I give NO guarantee that any of this is accurate.

% Feel free to pass this code on to whoever could use it.

%

% With any questions, and particularly suggestions for improvements please

% email me!

% Till Wagner (Sept 2016)

% tjwagner@ucsd.edu

% -------------------------------------------------------------------------

 

if nargin > 1 %if the download flag is set, this will download the nasateam monthly concentrations

    NOAA = ftp('sidads.colorado.edu');

    for yr = 1978:2015

        ms = 1; me = 12;

        if yr == 1978

            ms = 11;

        end

        for m = ms:me

            direct = sprintf('/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/north/monthly/');

            ttl = sprintf('nt_%d%02d_*_v1.1_n.bin',yr,m);

            cd(NOAA,direct)

            mget(NOAA,ttl)

        end

    end

    %now download the lat/lon grids

    cd(NOAA,'/pub/DATASETS/seaice/polar-stereo/tools/')

    mget(NOAA,'psn25lons_v3.dat')

    mget(NOAA,'psn25lats_v3.dat')

end

 

A_dims = [304,448]; %arctic SSMI grid dimensions

yr_min = 1979;      % first full year of data

yr_max = 2015;      % last full year of data

files = dir('nt_*_n.bin');

l =  length(files);

C_mon = cell(12,1);                 %set up cell arrays

n = zeros(12,1);

 

% From the data description:

% 0 - 250   Sea ice concentration (fractional coverage scaled by 250)

% 251   Circular mask used in the Arctic to cover the irregularly-shaped data gap around the pole (caused by the orbit inclination and instrument swath)

% 252   Unused

% 253   Coastlines

% 254   Superimposed land mask

% 255   Missing data

 

for i = 1:l                         %go through all months in data set

    fid = fopen(files(i).name);

    

    a = strsplit(files(i).name,'_'); aa =  a{2}; mon = str2double(aa(5:6));

    yr = str2double(aa(1:4));       %get year and month

    

    if yr >= yr_min && yr <= yr_max   

        fseek(fid, 300, 'bof');     %the first 300 lines are header

        A = fread(fid,A_dims);

        A2 = A;

        A2(A2>251)=nan;             %set everything that's not ocean to nan

        A2 = A2/250;                %rescale such that sea ice conc. is [0 1]

        

        n(mon) = n(mon)+1;

        if n(mon)==1                %save individual months in cell array

            C_mon{mon} = A2;

        else

            C_mon{mon} = C_mon{mon}+A2;

        end

    end

    

    if i == 3

        A_coast = A;                %get array with land mask

        A_coast(A_coast~=253)=0;

        A_coast(A_coast==253)=1;

        

        A_pole = A;                 %get array with pole mask

        A_pole(A_pole~=251)=0;

        A_pole(A_pole==251)=1;

        A_pole = logical(A_pole);

    end

    

    fclose(fid);

end

 

se  = strel('disk',1);

A_coast_nan = imdilate(A_coast,se);  %Creating a slightly expanded land mask

A_coast_nan(A_coast_nan==1)=nan;     %to remove land-ice interface from edge

 

remove_small = 10;                   %remove small blobs of ice edge

 

fid1 = fopen('psn25lons_v3.dat');

lons = fread(fid1,A_dims,'long')/1e5;

fclose(fid1);

fid0 = fopen('psn25lats_v3.dat');

lats = fread(fid0,A_dims,'long')/1e5;

fclose(fid0);

 

months = 1:12;

e_lon = cell(12,1); e_lat = e_lon;

 

for i = months

    % first get the mean concentration fields for each month

    a = C_mon{i}/n(i);                          %compute mean over period

    a2 = a; a2(a2<0.15) = nan; a2(A_pole) = nan;%set everything <15% to nan

    C_mon{i} = a2;

    % now compute the ice edge

    a(a>0.15)=1; a(a<0.15)=0; a(isnan(a))=0;    %make binary ice/no ice map

    a = imdilate(imerode(a,se),se);             %clean up map

    e = A_coast_nan+edge(a); e(isnan(e))=0;     %remove ice-coast edge

    e = bwareaopen(e,remove_small);             %get edge

    e_lon{i} = lons(e==1);                      %turn edge array to lat/lon

    e_lat{i} = lats(e==1);                      %positions

end

 

% now plot everything nicely ----------------------------------------------

figure(1); clf

mon_ttl = {'Jan','Feb','March','Apr','May','Jun','Jul','Aug','September','Oct','Nov','Dec'};

 

plot_mon= months_to_plot;

lp = length(plot_mon);

 

pole=[90 0 0];          %plot origin

latlim = [55,90];       %plot limits

 

load coastlines

for i = 1:lp %loop over chosen months

    subplot(max(1,round(lp/3)),min(lp,3),i)

    

    %make pretty subpanels

    axesm('stereo','MapLatLimit',latlim,'Origin',pole)

    fra = framem('on'); set(fra,'clipping','on');

    mlabel off; plabel off; gridm off

    set(gca,'PlotBoxAspectRatio','default')

    set(gca,'color','none'); set(gca,'XColor','none'); set(gca,'YColor','none')

    setm(gca,'FFaceColor',[.85 .95 1])      %ocean color

    

    ind = plot_mon(i);

    surfacem(lats,lons,C_mon{ind});   %plot concentration

    plotm(e_lat{ind},e_lon{ind},'r.','markersize',12)  %plot ice edge

    geoshow(coastlat, coastlon, 'DisplayType', 'polygon','facecolor',[.8 .8 .8]);

    title(mon_ttl{ind})

    

    caxis([0 1])

    

    if i ==lp

        pos = get(gca,'position');

        cc = colorbar; ylabel(cc,'sea ice concentration')

        set(gca,'position',pos)

    end

end

suptitle(sprintf('Climatological Mean Sea Ice Concentration, %d - %d',yr_min,yr_max))