%
% Illustrate how to use BYU SCP WorldLandBitMap.dat file
%

%%%%%%%%%%
% read file
fname='WorldLandBitMapLE.dat'; % little endian case
fid=fopen(fname,'r');
bmap=fread(fid,[36000/32 18000], '*uint32'); % for little endian machine
fclose(fid);

%%%%%%%%%%
% illustrate how to determine land flag for given lat/lon point(s)
lat=[ -10, 0, 20, -80,  80, 14]; % arbitrary locations
lon=[-151, 0, 11, -50,   5, 77];

NSX1=36000/32;
NSX=NSX1*32;
NSY=18000;

% vectorized
lon180=mod(lon+360,360);
lon180(lon180>0)=lon180(lon180>0)-360;
ix=round((lon180+180)*100);
ix(ix<0)=ix(ix<0)+NSX;
ix(ix>=NSX)=ix(ix>=NSX)-NSX;
ix1=floor(ix/32);
ix2=ix-ix1*32;
iy=round((lat+90)*100);
iy(iy<0)=0;
iy(iy>=NSY)=NSY-1;

land_bit=bitget(bmap(ix1+iy*NSX1+1),ix2,'uint32');

fprintf('land_bit: 1=land, 0=ocean\n');
for n=1:length(lon)
  fprintf('land_bit=%d  lon=%7.2f  lat=%6.2f  %8x %5d %4d %2d %5d\n',land_bit(n),lon(n),lat(n),bmap(ix1(n)+iy(n)*NSX1+1),ix(n),ix1(n),ix2(n),iy(n));
end
%%%%%%%%%%

% set an unsigned 32 bit integer to all 1s
all1s=uint32(1);
for k=2:32
  all1s=bitset(all1s,k,'uint32');
end
a32=1:32;

%%%%%%%%%%
% create a byte map to plot the binary mask (slow)
map2=zeros(NSX,NSY,'uint8'); % set byte array to zero (land)
for iy=1:18000
  if mod(iy,2000)==0
    fprintf('iy=%d\n',iy);
  end
  for ix=1:NSX1
    tb=bmap(ix+(iy-1)*NSX1);
    if tb~=0
      if tb==all1s  % all ocean
	map2((ix-1)*32+a32,iy)=1; % set ocean
      else         % only some ocean
	map2((ix-1)*32+a32,iy)=uint8(bitget(tb,a32,'uint32')); % set ocean
      end
    end
  end
end

%%%%%%%%%%
% plot converted byte map
myfigure(2)
imagesc([-180 180],[-90 90],rot90(map2));
title('Bit map');
colorbar;

% save map to figure
%print -dpng WorldLandBitMap.png
