
% module load license_matlab matlab
% matlab -nodesktop

% Script to demonstrate EUMETSAT Hybrid PC Representation for CrIS

% This script reads in SDR granules and computes and saves the Gobal PC scores (coefficients)
% and the Hybrid PCs and scores.
%
%
% David Tobin, CIMSS/SSEC/UW-Madison, 15 Nov 2023

  
addpath PCA_Mfiles/
addpath SDR_Reader/

  
% Load the Global PCs
% This is a sample file for J1/NOAA-20 file
gpcfile = 'globalPCs_version1_20220309.nc';
globalPC.U = ncread(gpcfile,'U');  % PCs (2223 x 150)
globalPC.M = ncread(gpcfile,'M');  % Mean spectrum (2223 x 1)
globalPC.D = ncread(gpcfile,'D');  % Eigenvalues (2223 x 1) only first 150 retained
v = ncread(gpcfile,'v');           % wavenumbers (2223 x 1)
nedn = ncread(gpcfile,'nedn');     % NEDN (2223 x 1) used for noise normalization


% CrIS SDR files.  Sample J1/NOAA-20 GCRSO and SCRIF files
pname = 'SDR_files/';

files = dir([pname 'SCRIF*.h5']);

% Loop over the SDR files, do the processing, and save PC data to .nc files
for iFile = 1:length(files)

  clear data globalCoef localCoef localPC RecRad Res RecRes

  % Load the granule and concatonate spectral bands
  d = concatsdr_pd(pname,files(iFile));
  data = [d.ES_RealLW(:,:);d.ES_RealMW(:,:);d.ES_RealSW(:,:)];

  % Project noise normalized radiance spectra onto the Global PCs
  [~,nPCs] = size(globalPC.U);
  nrec = size(data,2);
  nednm = repmat(nedn,1,nrec);
  globalCoef = pca_projectM(data./nednm,nPCs,globalPC);

  % Represent the Global Coefs (PC-Scores) as scaled integers.  Hard-coded here for 16 bits
  A = globalCoef(:,:);
  minA = min(A(:));
  rngA = range(A(:));
  B = uint16((A-minA)*(2^16-1)/rngA);
  globalCoef_min = minA;
  globalCoef_rng = rngA;
  globalCoef_int16 = B;
  globalCoef = single(globalCoef_int16)*globalCoef_rng/(2^16-1)+globalCoef_min;

  % Reconstruct radiances using the Global Coefs and Global PCs
  RecRad = pca_expandM(globalCoef,globalPC) .* nednm;

  % Compute the reconstruction residuals for the Global analysis/reconstruction
  Res = data - RecRad;

  % Now do PCA on these residuals, noise normalized
  localPC = pca_generate_pcs(Res./nednm);

  % Compute local coefs (pc-scores)
  nPCsLocal = 10;
  localPC.U = localPC.U(:,1:nPCsLocal);
  localCoef = pca_projectM(Res./nednm,nPCsLocal,localPC);

  % Represent the Local Coefs/Scores as scaled integers.  Hard coded here for 16 bits
  A = localCoef(:,:);
  minA = min(A(:));
  rngA = range(A(:));
  B = uint16((A-minA)*(2^16-1)/rngA);
  localCoef_min = minA;
  localCoef_rng = rngA;
  localCoef_int16 = B;
  localCoef = single(localCoef_int16)*localCoef_rng/(2^16-1)+localCoef_min;

  % and reconstruct
  RecRes = pca_expandM(localCoef,localPC) .* nednm;

  % The final reconstructed radiances are the Global reconstruction plus the
  % reconstructed local residuals
  RecRad = ( pca_expandM(globalCoef,globalPC) + pca_expandM(localCoef,localPC) ) .* nednm;

  % Put the reconstructed radiances into (nChannels,nFovs,nScans,nFORs) arrays
  sz = size(d.ES_RealLW);sz(1) = length(v);
  RecRad = reshape(RecRad,sz);

  % Put the Global and Local coefficients into (nChannels,nFovs,nFORs,nScans) arrays
  sz = size(d.ES_RealLW);sz(1) = nPCs;
  globalCoef_int16 = reshape(globalCoef_int16,sz);
  sz = size(d.ES_RealLW);sz(1) = nPCsLocal;
  localCoef_int16 = reshape(localCoef_int16,sz);

  % Save the global Coefficients, and the local residual PCs and corresponding 
  % local residual coefficients for this granule to file

  hpcfile = ['CRISPCAv1' files(iFile).name(6:end-11) 'ssec.nc'];
  [s,w] = unix(['rm ' hpcfile]);
  disp(hpcfile)
  sz = size(globalCoef_int16);
  nccreate(hpcfile,'globalCoef_int16','Dimensions',{'gpcs' sz(1) 'nFOVs' sz(2) 'nFORs' sz(3) 'nScans' sz(4) },...
         'Format','netcdf4','Datatype','uint16');
  nChannels = length(v);
  nccreate(hpcfile,'U','Dimensions',{'channels' nChannels 'pcs' nPCsLocal },'Format','netcdf4','Datatype','single');
  nccreate(hpcfile,'M','Dimensions',{'channels' nChannels },'Format','netcdf4','Datatype','single');
  nccreate(hpcfile,'D','Dimensions',{'channels' nChannels },'Format','netcdf4','Datatype','single');
  nccreate(hpcfile,'globalCoef_min','Dimensions',{'scalar' 1 },'Format','netcdf4','Datatype','single');
  nccreate(hpcfile,'globalCoef_rng','Dimensions',{'scalar' 1 },'Format','netcdf4','Datatype','single');

  sz = size(localCoef_int16);
  nccreate(hpcfile,'localCoef_int16','Dimensions',{'pcs' sz(1) 'nFOVs' sz(2) 'nFORs' sz(3) 'nScans' sz(4) },...
         'Format','netcdf4','Datatype','single');
  nccreate(hpcfile,'localCoef_min','Dimensions',{'scalar' 1 },'Format','netcdf4','Datatype','single');
  nccreate(hpcfile,'localCoef_rng','Dimensions',{'scalar' 1 },'Format','netcdf4','Datatype','single');

  % Uncomment to also save some other variables to the PC file
  %sz = size(d.lat);
  %nccreate(hpcfile,'lat','Dimensions',{'nFOVs' sz(1) 'nFORs' sz(2) 'nScans' sz(3) },'Format','netcdf4','Datatype','single');
  %nccreate(hpcfile,'lon','Dimensions',{'nFOVs' sz(1) 'nFORs' sz(2) 'nScans' sz(3) },'Format','netcdf4','Datatype','single');
  %nccreate(hpcfile,'sol_zen','Dimensions',{'nFOVs' sz(1) 'nFORs' sz(2) 'nScans' sz(3) },'Format','netcdf4','Datatype','single');
  %nccreate(hpcfile,'sat_zen','Dimensions',{'nFOVs' sz(1) 'nFORs' sz(2) 'nScans' sz(3) },'Format','netcdf4','Datatype','single');
  %nccreate(hpcfile,'sat_azi','Dimensions',{'nFOVs' sz(1) 'nFORs' sz(2) 'nScans' sz(3) },'Format','netcdf4','Datatype','single');
  %sz = size(d.asc_flag);
  %nccreate(hpcfile,'asc_flag','Dimensions',{'nScans' sz(1) },'Format','netcdf4','Datatype','single');
  %sz = size(d.obs_time_tai93);
  %nccreate(hpcfile,'obs_time_tai93','Dimensions',{'nFORs' sz(1) 'nScans' sz(2) },'Format','netcdf4','Datatype','single');

  ncwrite(hpcfile,'globalCoef_int16',globalCoef_int16);
  ncwrite(hpcfile,'globalCoef_min',globalCoef_min);
  ncwrite(hpcfile,'globalCoef_rng',globalCoef_rng);
  ncwrite(hpcfile,'U',localPC.U);
  ncwrite(hpcfile,'M',localPC.M);
  ncwrite(hpcfile,'D',localPC.D);

  ncwrite(hpcfile,'localCoef_int16',localCoef_int16);
  ncwrite(hpcfile,'localCoef_min',localCoef_min);
  ncwrite(hpcfile,'localCoef_rng',localCoef_rng);

  %ncwrite(hpcfile,'lat',d.lat);
  %ncwrite(hpcfile,'lon',d.lon);
  %ncwrite(hpcfile,'sol_zen',d.sol_zen);
  %ncwrite(hpcfile,'sat_zen',d.sat_zen);
  %ncwrite(hpcfile,'sat_azi',d.sat_azi);
  %ncwrite(hpcfile,'asc_flag',d.asc_flag);
  %ncwrite(hpcfile,'obs_time_tai93',d.obs_time_tai93);



  % Load data from the resulting PC file and compute reconstructed radiances.
  % These can be compared to the original data contained in structure d.
  [d2,gpc,lpc] = crisSDR_recrad(gpcfile,hpcfile,0);

  % Move PC file into PC_files/
  eval(['!mv ' hpcfile ' PC_files/'])


end

