clear all; close all; clc;
if isfolder('rlr_monthly'); delete('rlr_monthly'); end
delete *mat *txt

disp('Downloading newest PSMSL RLR-monthly data ...');
file = urlwrite('http://www.psmsl.org/data/obtaining/rlr.monthly.data/rlr_monthly.zip','rlr_monthly.zip');
unzip(file)
delete(file)
addpath(genpath('rlr_monthly'));

cd rlr_monthly

datev=datevec(date);
outputfile = ['rlr_',num2str(datev(1:3),'%d-%d-%d'),'.mat'];

disp('Loading PSMSL data files...');
data = readRlrMonthly('.');

disp('Preparing data matrices...');
outdata = struct;

num_gauges = length(data);
id_max = max([data.id]);
[id_sorted, id_inds] = sort([data.id]);

monthnumbers_required = [];
for i = 1:length(data)
    monthnumbers_required = union(monthnumbers_required, 12*data(i).year + data(i).month);
end

%convert to row vector, sort, and fill any gaps
monthnumbers_all = min(monthnumbers_required):max(monthnumbers_required);
num_timesteps = length(monthnumbers_all);

%initialize
outdata.id = NaN([num_gauges,1]); %(1:id_max)';
outdata.latitude = NaN([num_gauges,1]);
outdata.longitude = NaN([num_gauges,1]);
outdata.name = cell(num_gauges,1);
outdata.coastline = NaN([num_gauges,1]);
outdata.stationcode = NaN([num_gauges,1]);
outdata.stationflag = NaN([num_gauges,1]);
outdata.year = floor((monthnumbers_all-1)/12);
outdata.month = mod(monthnumbers_all-1,12)+1;
outdata.time = datenum(outdata.year, outdata.month, 15); %slightly rough hack
outdata.height = NaN([num_gauges,num_timesteps]);
outdata.missing = NaN([num_gauges,num_timesteps]);
outdata.dataflag = NaN([num_gauges,num_timesteps]);

for row = 1:num_gauges
%     id = data(i).id;
%     id = id_sorted(i);
    i = id_inds(row);
    monthnumbers = 12*data(i).year + data(i).month;
    [~,~,temporal_inds] = intersect(monthnumbers, monthnumbers_all);
    
    outdata.id(row) = data(i).id;
    outdata.latitude(row) = data(i).latitude;
    outdata.longitude(row) = data(i).longitude;
    outdata.name{row} = data(i).name;
    outdata.coastline(row) = data(i).coastline;
    outdata.stationcode(row) = data(i).stationcode;
    outdata.stationflag(row) = data(i).stationflag;
    outdata.height(row,temporal_inds) = (data(i).height - 7000) / 1000;
    outdata.missing(row,temporal_inds) = data(i).missing;
    outdata.dataflag(row,temporal_inds) = data(i).dataflag;
end

cd ..
save(outputfile,'-struct','outdata','-v7');

%% DEFINING LIMITS FOR TIDE GAUGE SELECTION
disp('Selecting TG-data for Arctic')

tg_data = load(outputfile);

lat_limit = tg_data.latitude > 60;
time_limit = tg_data.time > datenum('01-01-2000');
missing_limit = sum(isnan(tg_data.height(:,time_limit)),2)<100;

time_select = tg_data.time > datenum('01-01-1990');
id_select = missing_limit & lat_limit;

time_readable = datevec(tg_data.time(time_select));

tg_good = struct;

tg_good.id = tg_data.id(id_select);
tg_good.name = tg_data.name(id_select);
tg_good.height = tg_data.height(id_select,time_select);
tg_good.lat = tg_data.latitude(id_select);
tg_good.lon = tg_data.longitude(id_select);
tg_good.time = tg_data.time(time_select)/365.25;
tg_good.outliers = isoutlier(tg_good.height, 'mean');
tg_good.height_nan_outliers = tg_good.height;
tg_good.height_nan_outliers(tg_good.outliers) = NaN;
tg_good.height_outliers_and_gaps_filled = fillmissing(tg_good.height_nan_outliers,'pchip');

%% write to NETcdf file
disp('Write netCDF file and station list..')
ncfile = ['TG_DATA_',num2str(datev(1:3),'%d-%d-%d'),'.nc'];

if isfile(ncfile); delete(ncfile); end

sizex = size(tg_good.height,1);
sizey = size(tg_good.height,2);

nccreate(ncfile,'id','Dimensions',{'tg_locations',sizex});
ncwrite(ncfile,'id',tg_good.id)

nccreate(ncfile,'height','Dimensions',{'tg_locations',sizex,'months',sizey});
ncwrite(ncfile,'height',tg_good.height)

nccreate(ncfile,'lat','Dimensions',{'tg_locations',sizex});
ncwrite(ncfile,'lat',tg_good.lat)

nccreate(ncfile,'lon','Dimensions',{'tg_locations',sizex});
ncwrite(ncfile,'lon',tg_good.lon)

nccreate(ncfile,'time','Dimensions',{'months',sizey});
ncwrite(ncfile,'time',tg_good.time)

nccreate(ncfile,'outliers','Dimensions',{'tg_locations',sizex,'months',sizey});
ncwrite(ncfile,'outliers',single(tg_good.outliers))

nccreate(ncfile,'height_nan_outliers','Dimensions',{'tg_locations',sizex,'months',sizey});
ncwrite(ncfile,'height_nan_outliers',tg_good.height_nan_outliers)

nccreate(ncfile,'height_outliers_and_gaps_filled','Dimensions',{'tg_locations',sizex,'months',sizey});
ncwrite(ncfile,'height_outliers_and_gaps_filled',tg_good.height_outliers_and_gaps_filled)

T = cell2table(tg_good.name);
writetable(T,'tg_location_names.txt','WriteVariableNames',0)