%% access_Earth2014_grids1min.m
%
% This script shows how to access the data grids (e.g., surface
% topography, bedrock or shape models) of the Earth2014 models in Matlab. 
% 
% In detail, the script :
% 
%  a) reads the 1arc-min Earth2014 grid files
%  b) visualizes the data
%
% matlab version:  8.2.0.701 (R2013b)
% 
% authors : Christian Hirt, Moritz Rexer
% date : 2015 / 03 / 24
% contact: m.rexer@tum.de
% -------------------------------------------------------------------------
% All input grids are 10800 x 21600 pts in size (1 arc min)

clear all; close all

input_dir= '../data_1min/'

% Earth2014 models' file-names (.bshc)
mnames = [                                     % id
  'topo_grids/Earth2014.SUR2014.1min.geod.bin                     '; % 1  : SUR - topo grids
  'topo_grids/Earth2014.BED2014.1min.geod.bin                     '; % 2  : BED
  'topo_grids/Earth2014.ICE2014.1min.geod.bin                     '; % 3  : ICE
  'topo_grids/Earth2014.RET2014.1min.geod.bin                     '; % 4  : RET
  'topo_grids/Earth2014.TBI2014.1min.geod.bin                     '; % 5  : TBI
  'shape_grids/Earth2014Shape_minus_6371000m.SUR2014.1min.geod.bin'; % 6  : BED - shape grids
  'shape_grids/Earth2014Shape_minus_6371000m.BED2014.1min.geod.bin'; % 7  : ICE
  'shape_grids/Earth2014Shape_minus_6371000m.RET2014.1min.geod.bin'; % 8  : RET
  'shape_grids/Earth2014Shape_minus_6371000m.TBI2014.1min.geod.bin'; % 9  : TBI
  'masks/MSK2014_inputdatasets.1min.geod.bin                      '; % 10 : Input data mask
  'masks/MSK2014_landtypes.1min.geod.bin                          '; % 11 : Land type mask
  ];

% grid definitions
ylimvec = [-90 90]; xlimvec = [-180 180];
res_deg = 1/60;
vec_lat = -90+res_deg/2   :res_deg:90-res_deg/2;
vec_lon = -180+res_deg/2  :res_deg:180-res_deg/2;
n = length(vec_lat); m = length(vec_lon);

for i = 1:11 % use desired model ids here (1 to 11)
 
   % read file 
   filename = strcat(strcat(input_dir,strtrim(mnames(i,:))));
   fid = fopen(filename,'r');
   [data, cnt] = fread(fid, inf, 'int16','ieee-be');
   fclose(fid); data = reshape(data,m,n)'; 

   if (i>=6 && i<=9)  % these are the shape files
      data = data + 6371000; % add constant to data to obtain radii
   end
   
   figure(1)
   imagesc(vec_lon,vec_lat,data);
   set(gca,'XDir','normal','YDir','normal');
   colorbar; grid on; xlabel('Longitude [deg]');ylabel('Latitude [deg]')
   title(strcat('File: ',strrep(filename,'_',' ')))
   print(strrep(filename,'bin','png'),'-dpng','-r300');

end

