proj_PS_EG_main

%% project netCDF with sea ice parameter and its corresponding lat/lon matrix into
% PolarStereographic and EASE_Grid2.0 with nearest sampling, marginal NaN as landmask
% NaN value specification: MissingValue = -1, FillValue = -999;
% author: Teng Li, litengbnu@foxmail.com 20170726
% =========================================================
% 20170905 modified to accommodate Antarctic and thickness colormap;
% =========================================================
% 20170910 modified to extract cartography snippet to new function;
% Projection specification according to EPSG code(all based on WGS-84 ellipsoid):
% WKID 3413 for NSIDC-PS-North, which is (-45, 70)
% WKID 3976 for NSIDC-PS-Soutn, which is (0, -70)
% WKID 3973 for EASE-Grid-North, which is (0, 90)
% WKID 3974 for EASE-Grid-South which is (0, -90)
% =========================================================
% 20171005 modified according to Mr. WenHuan Chen's suggestion:
% (1) add 'try...catch' to tolerate error in EVERY(?) function;
% (2) =!=alter input 'nc_path' as filename instead of full path;
% (2.1)#####################delete the second#####################
% (3) chose only indispensable function from m_map library.
% =========================================================
% 20171010 modified to combine with YiFan Ding's demo.m together.
% delete the 'try..catch' structure and output flag and file-name.
% 20171023: add mode = carto/normal to convinience of testing by comment.
function proj_PS_EG_main(output_space, inter_nc_path, resolution, mode)
% input: NetCDF file generated by YiFan-Ding, with sea-ice parameter and lat/lon
% output: (1) message to flag execution successful or failed(enter into 'catch-block')
% (2) two NetCDF file, projected into PS and EASE-Grid2.0 of 6250 meter with variable structure respectively.
% one GeoTIFF(optional) and its PNG file, which is projected into PS for cartography
% example: proj_PS_EG('D:\save-path', 'S_AMSR2_20170711_SIC.nc', 6250);
%% construct output filenames, including NetCDF, GeoTIFF, and PNG
[~, nc_file_name] = fileparts(inter_nc_path);
nc_PS_name = [nc_file_name, '_PS', '.nc'];
nc_EG_name = [nc_file_name, '_EG', '.nc'];
% tiff_PS_name = [nc_file_name, '_PS', '.tif'];
% tiff_EG_name = [nc_file_name, '_EG', '.tif'];
png_name = [nc_file_name, '_PS', '.png'];
% seach for variable which does not start with 'lat' or 'lon'
% assume NetCDF file has only 3 variable: lat / lon / parameter in format
var_name = ncinfo(inter_nc_path);
var_name = extractfield(var_name.Variables, 'Name');
if numel(var_name) == 3
    for nn = 1:3
        candi = var_name{nn};
        if ~(strcmp(candi(1:3), 'lat') || strcmp(candi(1:3), 'lon'))
            parameter = candi;
        end
    end
end
% construct date string from input filename for cartography
file_name_split = strsplit(nc_file_name, '_');
polar = file_name_split{1};
sensor = file_name_split{2};
date_hour = datenum(file_name_split{3}, 'yyyymmdd');
date_string = datestr(date_hour);
%% read sea ice parameter and lat/lon of NetCDF from YiFan-Ding
% padding marginal -999(replace originally NaN) to avoid resampling incorrectly.
[mrg_lat, mrg_lon, mrg_par, scale] = margin_nan(inter_nc_path, parameter);
% calculate orginal irregular grids of PS and EG coordinate from original lat/lon
major_axis = 6378137.0;
eccentricity = 0.08181919;
if strcmp(polar, 'N')
    central_meridian = -45;
    true_latitude = 70;
else
    central_meridian = 0;
    true_latitude = -70;
end
[x_PS_irr, y_PS_irr] = polarstereo_fwd(mrg_lat, mrg_lon, major_axis, eccentricity, true_latitude, central_meridian);
if ~strcmp(mode, 'carto')
    [x_EG_irr, y_EG_irr] = latlon2ease_isprs(mrg_lat, mrg_lon, polar);
end
%% generate regular grids of 6250m in PS and EG projection, respectively.
%  !!Attention to the Order of Row & Column!!
[x_PS_reg, y_PS_reg] = regular_grid(x_PS_irr, y_PS_irr, resolution);
if ~strcmp(mode, 'carto')
    [x_EG_reg, y_EG_reg] = regular_grid(x_EG_irr, y_EG_irr, resolution);
end
% reproject and resample by griddata of nearest method, longest waiting!
par_PS = griddata(x_PS_irr, y_PS_irr, mrg_par, x_PS_reg, y_PS_reg, 'nearest');
if ~strcmp(mode, 'carto')
    par_EG = griddata(x_EG_irr, y_EG_irr, mrg_par, x_EG_reg, y_EG_reg, 'nearest');
end

%% write NetCDF and GeoTiff file with geoferenced data and coordinate dimention
if ~exist(output_space, 'dir')
    mkdir(output_space)
end
% check then delete first if there is file with same name(!!dangerous operation!!)
if exist(fullfile(output_space, nc_PS_name), 'file')
    delete(fullfile(output_space, nc_PS_name));
end
if exist(fullfile(output_space, nc_EG_name), 'file')
    delete(fullfile(output_space, nc_EG_name));
end
% construct a struct to store the attributes of the parameter, which is
% then transfered to write_netcdf() as argument for final NetCDF product.
attri_nc.long_name = ncreadatt(inter_nc_path, ['/',parameter], 'LongName');
attri_nc.unit = ncreadatt(inter_nc_path, ['/',parameter], 'Unit');
attri_nc.valid_range = ncreadatt(inter_nc_path, ['/',parameter], 'ValidRange');
attri_nc.scale = scale;
attri_nc.polar = polar;
if ~strcmp(mode, 'carto')
    write_netcdf(fullfile(output_space, nc_PS_name), parameter, x_PS_reg, y_PS_reg, par_PS, attri_nc, 'PS');
    write_netcdf(fullfile(output_space, nc_EG_name), parameter, x_EG_reg, y_EG_reg, par_EG, attri_nc, 'EG');
end
% write GeoTIFF File for cartography, only in PS45
% 20170910 modified: no need to export GeoTIFF although function is already done.
% write_geotiff(fullfile(output_space,tiff_PS_name), polar, 'PS', x_PS_reg, y_PS_reg, par_PS)
% write_geotiff(fullfile(output_space,tiff_EG_name), polar, 'EG', x_EG_reg, y_EG_reg, par_EG)
%% construct attributes for cartography then call carto_png
attri_carto.polar = polar;
attri_carto.sensor = sensor;
attri_carto.parameter = parameter;
attri_carto.unit = attri_nc.unit;
attri_carto.date_string = date_string;
attri_carto.resolution = resolution;
fig = carto_png(x_PS_reg, y_PS_reg, par_PS, attri_carto);
print(fig, fullfile(output_space, png_name), '-dpng', '-r300');
delete(fig)
end