%
% Matlab script to illustrated how plot and use the 
% BYU SCP WorldLandDistMap. dat file
%

% open and read file
fname='WorldLandDistMap.dat';
fid=fopen(fname,'r');
map1=fread(fid,[36000 18000], '*uint8');
fclose(fid);

% plot file  note rotation of axises so that it plots properly
figure(1)
imagesc([-180 180],[-90 90],rot90(map1));
colorbar
title('WorldLandDistMap: Truncated distance to nearest land in km')

%%%%%%%%%%%%%%%%
% show how to use lat/lon to select pixels
lat=[ -10, 0, 20, -80,  80, 14]; % arbitrary locations
lon=[-151, 0, 11, -50,   5, 77];

% show locations on map
hold on;
plot(lon,lat,'r.');
hold off;

NSX=36000;
NSY=18000;

% vectorized computation of indexes into array
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;
iy=round((lat+90)*100);
iy(iy<0)=0;
iy(iy>=NSY)=NSY-1;

% fill land distance array
land_dist=map1(ix+1,iy+1);

% print
fprintf('Distances longer than 255 km are truncated to 255 km\n');
for n=1:length(lon)
  fprintf('land_dist(km)=%3d  lon=%7.2f  lat=%6.2f  %5d %5d\n',land_dist(n),lon(n),lat(n),ix(n),iy(n));
end
