% read and display Tarfala DEM from binary Golden Software Surfer grid
% peter.jansson&rickard.pettersson@2003-09-29
% Please note that Golden Software Surfer binary files have different
%  formats for different versions o fthe program. this script works with
%   Surfer v. 7 & 8
clear

% open binary file
filename = 'trs_dem.grd';
fid = fopen(filename);

% Reads size of grid to a vector (nR(y),nC(x))
fseek(fid,20,'bof');
gsize = fread(fid,[1 2],'long');

% Reads coordinates for lower left corner of grid to vector (xLL,yLL)
fseek(fid,28,'bof');
LL = fread(fid,2,'double');

% Reads grid cell size (xSize,ySize)
fseek(fid,44,'bof');
cellsize = fread(fid,2,'double');

% Reads range of Z values to a vector (Zmin,Zmax)
fseek(fid,60,'bof');
z_range = fread(fid,2,'double');

% Reads blankvalue used by Surfer
fseek(fid,84,'bof');
blankvalue_surf = fread(fid,1,'double');

% Reads grid data
fseek(fid,100,'bof');
Zgrid = fread(fid,[gsize(2) gsize(1)],'double');
% I need to rotate the grid by -90 deg and flip it left-to-right to plot DEM 
% using contourf and surf commands
Zgrid =fliplr(rot90(Zgrid,-1));
i = find(Zgrid == blankvalue_surf);
Zgrid(i) = nan;

x = (LL(1):cellsize(1):(LL(1)+(gsize(2)-1)*cellsize(1)));
y = (LL(2):cellsize(2):(LL(2)+(gsize(1)-1)*cellsize(2)));
%[Xgrid,Ygrid] = meshgrid(x,y);

fclose(fid);

% Colormap for surface
top =...
[1.0000    1.0000         0;
    0.8784    0.9569         0;
    0.7647    0.9176         0;
    0.6549    0.8745         0;
    0.5569    0.8353         0;
    0.4627    0.7922         0;
    0.3765    0.7529         0;
    0.2941    0.7098         0;
    0.2235    0.6667         0;
    0.1569    0.6275         0;
    0.0980    0.5843         0;
    0.0471    0.5451         0;
         0    0.5020         0;
    0.0471    0.5020         0;
    0.0941    0.5020         0;
    0.1412    0.5020         0;
    0.1882    0.5020         0;
    0.2353    0.5020         0;
    0.2824    0.5020         0;
    0.3294    0.5020         0;
    0.3765    0.5020         0;
    0.4235    0.5020         0;
    0.4706    0.5020         0;
    0.5020    0.4863         0;
    0.5020    0.4392         0;
    0.5020    0.3922         0;
    0.5020    0.3451         0;
    0.5020    0.2980         0;
    0.5020    0.2510         0;
    0.5137    0.2588    0.0235;
    0.5255    0.2627    0.0471;
    0.5373    0.2745    0.0745;
    0.5490    0.2824    0.0980;
    0.5608    0.2941    0.1255;
    0.5686    0.3059    0.1569;
    0.5804    0.3216    0.1843;
    0.5922    0.3373    0.2157;
    0.6039    0.3529    0.2471;
    0.6157    0.3725    0.2784;
    0.6275    0.3922    0.3137;
    0.6392    0.4157    0.3490;
    0.6510    0.4392    0.3843;
    0.6627    0.4667    0.4196;
    0.6745    0.4941    0.4588;
    0.6863    0.5216    0.4980;
    0.6941    0.5569    0.5373;
    0.7059    0.5922    0.5804;
    0.7176    0.6275    0.6196;
    0.7294    0.6667    0.6627;
    0.7412    0.7098    0.7059;
    0.7529    0.7529    0.7529;
    0.7725    0.7725    0.7725;
    0.7922    0.7922    0.7922;
    0.8118    0.8118    0.8118;
    0.8275    0.8275    0.8275;
    0.8471    0.8471    0.8471;
    0.8667    0.8667    0.8667;
    0.8863    0.8863    0.8863;
    0.9059    0.9059    0.9059;
    0.9255    0.9255    0.9255;
    0.9412    0.9412    0.9412;
    0.9608    0.9608    0.9608;
    0.9804    0.9804    0.9804;
    1.0000    1.0000    1.0000];


% make DEM plot
% Establish screen values for plots
screenrect   = get(0,'screensize');
screenwidth  = screenrect(3);
screenheight = screenrect(4);
figwidth     = floor(0.9*screenwidth);
figheight    = floor(0.9*screenheight);

% Plot topographic map
figure('position', [(screenwidth/2-figwidth/2) (screenheight/2-figheight/2) figwidth figheight],...
        'units','pixels','NumberTitle','off','Name','Tarfala DEM topo');
 orient landscape
  set(gcf, 'PaperUnits', 'centimeters');
   set(gcf, 'PaperType', 'A4');

% Set contour intervals
v = 500:25:2200;
contourf(x,y,Zgrid,v)
 colormap(top)

% Plot perspective view surface 
figure('position', [(screenwidth/2-figwidth/2) (screenheight/2-figheight/2) figwidth figheight],...
        'units','pixels','NumberTitle','off','Name','Tarfala DEM surface');
 orient landscape
  set(gcf, 'PaperUnits', 'centimeters');
   set(gcf, 'PaperType', 'A4');

 h = surf(x,y,Zgrid);
  colormap(top)
   shading interp
   %set(h,'EdgeColor','k')
    light('Position',[-2,2,20])
     lighting phong

