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)
% -------------------------------------------------------------------------
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))