The answer depends on exactly what you want to display. Do you want 12 different maps, to show salinity distribution for each month? If so, this would do it:
for k = 1:12
subplot(3,4,k)
imagesc(x,y,sisali(:,:,k))
axis xy
title(['Month ',num2str(k)])
caxis([lowerlimit upperlimit]) % for consistency across all plots
cmocean haline % optional colormap
end
Or do you simply want a single map showing average salinity for the year? If so,
imagesc(x,y,mean(sisali,3))
axis xy
cmocean haline
You could similarly subtract mean salinity from each month to show monthly anomalies, like this.
for k = 1:12
subplot(3,4,k)
imagesc(x,y,sisali(:,:,k)-mean(sisali,3))
axis xy
title(['Month ',num2str(k)])
caxis([-1 1]*limit) % for consistency across all plots
cmocean balance % divergent colormap
end
h = imagesc(x,y,sisali(:,:,k));
axis xy
caxis([lowerlimit upperlimit])
cmocean haline
title('Month 1')
gif('amundsen_sea_salinity.gif')
for k = 2:12
h.CData = sisali(:,:,k);
gif
end
Or if you simply want to plot a time series of the average salinity in a given region, you could create a mask for that region that is true for any grid cells you're interested in. An easy way to do that might be to grid up your 1D arrays of coordinates and then use my isopenocean function to determine which of those grid cells correspond to ocean:
[X,Y] = meshgrid(x,y);
mask = isopenocean(X,Y);
Then use my local function to average all the grid cells in that region, for each month. Then you'll end up with a time series of 12 values, containing the mean salinity within the masked region for each month of the year:
sal = local(sisali,mask);
plot(1:12,sal)
xlabel 'month of year'
ylabel 'average salinity within region of interest'
