
%-------------------------------------------------------------------------
% EurAlpGM2011 read test 
%------------------------------------------------------------------------- 
%
% This Matlab-script reads and displays the content of the 15 EurAlpGM2011
% product files that can be downloaded from the Curtin website.
%
% 
%-------------------------------------------------------------------------
% Dr Christian Hirt
% Western Australian Centre for Geodesy
% Curtin University, Perth
%-------------------------------------------------------------------------
% Version of 27/02/2012
%-------------------------------------------------------------------------
clear all
close all

set(0,'DefaultAxesFontSize',12);
set(0,'DefaultAxesLineWidth',0.8);
set(0,'DefaultLineLineWidth',0.8); 


%-------------------------------------------------------------------------
% Paths (can be customized)
%-------------------------------------------------------------------------
%EurAlpGM2011data_dirname = 'EurAlpGM2011_data_files/'
%EurAlpGM2011img_dirname  = 'EurAlpGM2011_pictures/'
EurAlpGM2011data_dirname = ''
EurAlpGM2011img_dirname  = ''

%-------------------------------------------------------------------------
% Grid definition 
%-------------------------------------------------------------------------

resolution = 0.002               % in deg
minlon =    5+resolution/2       % Longitude min, max
maxlon =   17-resolution/2
minlat =   44+resolution/2       % Latitude min, max
maxlat =   48-resolution/2

vec_lon = minlon:resolution:maxlon;
vec_lat = minlat:resolution:maxlat;
n = length(vec_lat);
m = length(vec_lon);



%-------------------------------------------------------------------------
% Description of the EurAlpGM2011 data files 
%-------------------------------------------------------------------------
% file name convention for publication:
% Productclass_functional_[remarks].[spectralband].unit.resolution[arcmin].bin

% filenames
vec_product_filenames = { 
              'EGM2008_zeta.atTopo.unit_mm.7.2s.int32'
              'EGM2008_gravdist.atTopo.unit_tenthsofmGal.7.2s.int16'
              'EGM2008_NSDOV.atTopo.unit_tenthsofarcsec.7.2s.int16'
              'EGM2008_EWDOV.atTopo.unit_tenthsofarcsec.7.2s.int16'    
              'ERTM2160_zeta.unit_mm.7.2s.int16' 
              'ERTM2160_gravdist.unit_tenthsofmGal.7.2s.int16' 
              'ERTM2160_NSDOV.unit_tenthsofarcsec.7.2s.int16' 
              'ERTM2160_EWDOV.unit_tenthsofarcsec.7.2s.int16'              
              'EurAlpGM2011_zeta.unit_mm.7.2s.int32'
              'EurAlpGM2011_gravdist.unit_tenthsofmGal.7.2s.int16' 
              'EurAlpGM2011_NSDOV.unit_tenthsofarcsec.7.2s.int16' 
              'EurAlpGM2011_EWDOV.unit_tenthsofarcsec.7.2s.int16'          
              'EurAlpGM2011_normalgravity.unit_tenthsofmgal.7.2s.int32'
              'EurAlpGM2011_surfacegravity.unit_tenthsofmgal.7.2s.int32'
              'EurAlpGM2011_totalDOV.unit_tenthsofarcsec.7.2s.int16'          
}


  
% column 1: conversion factor to convert data from file unit to
%           basis units mgal, m, arc seconds
vec_constants = [     
    0.001    % mm to m
    0.1      % to mGal
    0.1      % to sec 
    0.1      % to sec 
    0.001    % mm to m
    0.1      % to mGal
    0.1      % to sec 
    0.1      % to sec 
    0.001    % mm to m
    0.1      % to mGal
    0.1      % to sec 
    0.1      % to sec 
    0.1      % to mGal
    0.1      % to sec 
    0.1      % to sec 
    ]                 

% product class, functional name, unit and file format
vec_class_funct_unit_format = {                          
    'EGM2008'      'zeta'       '[m]'             'int32'
    'EGM2008'      'gravdist'   '[mGal]'          'int16'
    'EGM2008'      'NSDoV'      '[arc sec]'       'int16'
    'EGM2008'      'EWDoV'      '[arc sec]'       'int16'
    'ERTM2160'     'zeta'       '[m]'             'int16'
    'ERTM2160'     'gravdist'   '[mGal]'          'int16'
    'ERTM2160'     'NSDoV'      '[arc sec]'       'int16'
    'ERTM2160'     'EWDoV'      '[arc sec]'       'int16'
    'EurAlpGM2011' 'zeta'       '[m]'             'int32'
    'EurAlpGM2011' 'gravdist'   '[mGal]'          'int16'
    'EurAlpGM2011' 'NSDoV'      '[arc sec]'       'int16'
    'EurAlpGM2011' 'EWDoV'      '[arc sec]'       'int16'
    'EurAlpGM2011' 'gamma'      '[mGal]'          'int32'
    'EurAlpGM2011' 'gravity'    '[mGal]'          'int32'
    'EurAlpGM2011' 'totalDoV'   '[arc sec]'       'int16'
  }                                                      


%-------------------------------------------------------------------------
% Load and display all data files 
%-------------------------------------------------------------------------
list_of_files = [1 2 15]  
%list_of_files = 1:15
% This list can be customized 
% 1  = first file  (EGM2008  zeta evaluated at the topography)
% 2  = second file (EGM2008  gravdist evaluated at the topography)
% 15 = last file   (EurAlpGM2011 total DoV)

for i = 1:length(list_of_files)
    
  ix =  list_of_files(i)
    
  inputfilename  = char(vec_product_filenames(ix))
  fac            = vec_constants(ix,1)
  classname      = char(vec_class_funct_unit_format(ix,1))
  functname      = char(vec_class_funct_unit_format(ix,2))
  unitname       = char(vec_class_funct_unit_format(ix,3))
  formatstring   = char(vec_class_funct_unit_format(ix,4))
  
  % Read data and convert --------------------------
  disp('Read input file')
  fid = fopen(strcat(EurAlpGM2011data_dirname,inputfilename,'.bin'),'rb')
  data = fread(fid,n*m,formatstring,'ieee-be');
  fclose(fid);
  data = data*fac;    %conversion from file unit to basis unit
  data = reshape(data,n,m);

  % Display ----------------------------------------
  figure(1)  
  subplot(2,1,1)
  imagesc(vec_lon,vec_lat,data);
  set(gca,'XDir','normal','YDir','normal');
  xlabel('Longitude [deg]')
  ylabel('Latitude [deg]')
  titlestring = [classname ' ' functname ' ' unitname]
  title(titlestring)
  colorbar
  grid on
  ylim([minlat maxlat])
  xlim([minlon maxlon])
  print(strcat(EurAlpGM2011img_dirname,inputfilename,'.png'),'-r400','-dpng')
   
  
  clear data
  
  close all
end

disp('EurAlpGM2011 read test finished.')












