%% Spatiotemporal correlations
%{
Creator: Austin M. Weber
Instituion: Byrd Polar & Climate Research Center, Columbus, OH

The spatiotemporal data in this study were acquired as NetCDF (.nc) files
==============================================================================================
PARAMETER     || DATA SET    || LINK
==============================================================================================
Sea Surface   || ERSST       || https://www.ncei.noaa.gov/products/extended-reconstructed-sst
Temperature   ||             ||
              ||             ||
500 hPa       || ERA5        || https://cds.climate.copernicus.eu		
Temperature   ||             ||
              ||             ||
Precipitation || GPCCv2020   || https://psl.noaa.gov/data/gridded/data.gpcc.html
              ||             ||	
===============================================================================================

Only the October through April months were acquired, for the years 1960
through 2019 (420 months total).

This script requires the user to have installed the Climate Data Toolbox (https://www.mathworks.com/matlabcentral/fileexchange/70338-climate-data-toolbox-for-matlab)
%}

%
% SCRIPT BEGINS HERE
%

% Extract data from .nc files
% Adjust strings as necessary depending on the file
filename = 'data.nc';
lat = ncread(filename,'latitude'); % May be 'lat', depending on the .nc file
lon = ncread(filename,'longitude'); % May be 'lon', depending on the .nc file
climate_monthly = ncread(filename,'climate parameter'); % e.g., 'sst' for sea surface temperatures

% Convert spatiotemporal data from monthly to annual
num_months = size(climate_monthly,3);
num_years = numel(1960:2019);
num_months_per_season = 7;
start_month = linspace(1, num_months-6, num_years);
end_month = linspace(num_months_per_season, num_months, num_years);
climate_annual = zeros(length(lon), length(lat), num_years);
for k = 1:num_years
	temp = mean(climate_monthly(lon, lat, start_month(k):end_month(k), 3));
	climate_annual(:,:,k) = temp;
end

% Smooth data and detrend
% Spatiotemporal data
climate_annual_3YRM = movmean(climate_annual,3,3);
climate_annual_3YRM_dt = detrend3(climate_annual_3YRM); 
% Ice core data (e.g., d18O)
ice_core_data_3YRM = movmean(ice_core_data,3);
ice_core_data_3YRM_dt = detrend(ice_core_data_3YRM);

% Correlations between annual ice core data and the spatiotemporal data
% NOTE: Make sure that ice core data and spatiotemporal data are both
% temporally aligned prior to this step (e.g., both data sets start 1960
% and end 2019).
r_vals = zeros(length(lon),length(lat));
p_vals = zeros(length(lon),length(lat));
for i = 1:length(lon)
for j = 1:length(lat)
	[r_temp,p_temp] = corr(ice_core_3YRM_dt, climate_annual_3YRM_dt(i,j,:));
	r_vals(i,j) = r_temp;
	p_vals(i,j) = p_temp;
end
end

% Plot the statistically significant (p <= 0.05) data
p_idx = p_vals <= 0.05;
p_idx = double(p_idx);
p_idx = standardizeMissing(p_idx,0);
r_significant = r_vals;
r_significant = r_significant.*p_idx;
figure
imagescn(lat,lon,r_significant)
cmocean('balance','pivot') 
cb=colorbar;
cb.Label.String='Correlation (r)';
axis equal tight
view([90,-90])
% Use the borders() function from the Climate Data Toolbox to plot world
% boundaries


