function D=gshhs_read(optn,file,varargin);
%GSHHS_READ   load GSHHS database
%
%  GSHHS_READ  Loads various GSHHS coastline databases and does some preliminary
%              processing to get things into the form desired by the patch-filling
%              algorithm.
%
% Thus function is based on the function mu_coast.m from the m_map toolbox
% (http://www.eos.ubc.ca/~rich/map.html) of Rich Pawlowicz.
%
%See also: M_COAST (m_map), GSHHS ($ Matlab mapping toolbox)

   % use 10^6 degrees internally
   % choose to cover entire world
   
   llim =   190.5e6; %180e6;
   rlim =   190.5e6; %180e6;
   blim =   -90e6;
   tlim =    90e6;
   
   mrlim = -169.5; %-180;
   mllim = -169.5; %-180;
   mblim =  -90;
   mtlim =   90;
   
%% pre-allocate

   D.GSHHS_DATA_VERSION = 6;      % /* For v1.5 data set */
   D.GSHHS_PROG_VERSION = '1.9';  %
   D.GSHHS_SCL          = 1.0e-6; % /* COnvert micro-degrees to degrees */
  %D.area               = [];     % /* Area of polygon in 1/10 km^2 */
  %D.lon                = [];     %/ * Each lon, lat pair is stored in micro-degrees in 4-byte integer format */
  %D.lat                = [];
   
   % decfac is for decimation of areas outside the lat/long bdys.
   % Sizes updated for gshhs v1.3
   switch optn(1),  %makre sure pre-allocation is enough
     case 'f',D.lon=NaN+zeros(1,11340934);D.lat=NaN+zeros(1,11340934);D.area=zeros(1,188613);D.k=ones(1,188614);decfac=12500;
     case 'h',D.lon=NaN+zeros(1, 2063513);D.lat=NaN+zeros(1, 2063513);D.area=zeros(1,153545);D.k=ones(1,153546);decfac= 2500;
     case 'i',D.lon=NaN+zeros(1,  493096);D.lat=NaN+zeros(1,  493096);D.area=zeros(1, 41529);D.k=ones(1, 41530);decfac=  500;
     case 'l',D.lon=NaN+zeros(1,  101207);D.lat=NaN+zeros(1,  101207);D.area=zeros(1, 10775);D.k=ones(1, 10776);decfac=  100;
     case 'c',D.lon=NaN+zeros(1,   14884);D.lat=NaN+zeros(1,   14884);D.area=zeros(1,  1870);D.k=ones(1,  1871);decfac=   20;
   end;
   
   D.id                 = zeros(size(D.area)); % /* Unique polygon id number, starting at 0 */
   D.n                  = zeros(size(D.area)); % /* Number of points in this polygon */
                                               % /* = level + version << 8 + greenwich << 16 + source << 24 */
                                               % /* flag contains 4 items, one in each byte, as follows:
   D.level              =  int8(size(D.area)); %  * low byte:	level     =  flag        & 255: Values: 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
   D.version            =  int8(size(D.area)); %  * 2nd byte:	version   = (flag >> 8 ) & 255: Values: Should be 4 for GSHHS version 1.4
   D.greenwich          =  int8(size(D.area)); %  * 3rd byte:	greenwich = (flag >> 16) & 255: Values: Greenwich is 1 if Greenwich is crossed
   D.source             =  int8(size(D.area)); %  * 4th byte:	source    = (flag >> 24) & 255: Values: 0 = CIA WDBII, 1 = WVS
                                               %  */
   D.west               = zeros(size(D.area)); %  /* min/max extent in micro-degrees */
   D.east               = zeros(size(D.area));        
   D.south              = zeros(size(D.area));        
   D.north              = zeros(size(D.area));
   
%% load   
   
   fid   = fopen(file,'r','ieee-be');
   area2 = D.area;
   
   % Read the File header
   %%[A,cnt]=fread(fid,9,'int32');
    [A,cnt]=get_gheader(fid);
   
   
   nu = 0; % number of patches in file
   l  = 0; % number of patches in specified axes box
   while cnt>0,
   
   nu = nu + 1;
   
    % A: 1:ID, 2:num points, 3:land/lake etc., 4-7:w/e/s/n, 8:area, 9:greenwich crossed.
    
    w = A(4);
    e = A(5);
    s = A(6);
    n = A(7);
    
    C=fread(fid,A(2)*2,'int32'); % Read all points in the current segment.
   
    %For Antarctica the lime limits are 0 to 360 (exactly), thus c==0 and the
    %line is not chosen for (e.g. a conic projection of part of Antarctica)
    % Fix 30may/02
    if e==360e6, e=e-1; end;
    
    a = rlim>llim;  % Map limits cross longitude jump? (a==1 is no)
    b = A(9)<65536; % Cross boundary? (b==1 if no).
    c = llim<rem(e+360e6,360e6); 
    d = rlim>rem(w+360e6,360e6);
    e = tlim>s & blim<n;
    
    % This test checks whether the lat/long box containing the line overlaps that of
    % the map. There are various cases to consider, depending on whether map limits
    % and/or the line limits cross the longitude jump or not.
    
   %% if e & (  a&( b&c&d | ~b&(c|d)) | ~a&(~b | (b&(c|d))) ),
    if e & (  (a&( (b&c&d) | (~b&(c|d)) )) | (~a&(~b | (b&(c|d))) ) ),
    
      l=l+1;
    
      x=C(1:2:end)*1e-6;y=C(2:2:end)*1e-6;
    %%  plot(x,y);pause;
      
      %  make things continuous (join edges that cut across 0-meridian)
   
      dx=diff(x);
   %%fprintf('%f %f\n', max(dx),min(dx))
   %   if A(9)>65536 | any(dx)>356 | any(dx<356),
   %  if ~b | any(dx>356) | any(dx<-356),
        x=x-360*cumsum([x(1)>180;(dx>356) - (dx<-356)]);
   %   end;
   
      % Antarctic is a special case - extend contour to make nice closed polygon
      % that doesn't surround the pole.   
      if abs(x(1))<1 & abs(y(1)+68.9)<1,
        y=[-89.9;-78.4;y(x<=-180);y(x>-180);   -78.4;-89.9*ones(18,1)];
        x=[  180; 180 ;x(x<=-180)+360;x(x>-180);-180; [-180:20:160]'];
      end;
   
      % First and last point should be the same.
      
      if x(end)~=x(1) | y(end)~=y(1), x=[x;x(1)];y=[y;y(1)]; end;
   
      % get correct curve orientation for patch-fill algorithm.
      
       area2(l)=sum( diff(x).*(y(1:(end-1))+y(2:end))/2 );
      D.area(l)=A(8)/10;
   
      if rem(A(3),2)==0; 
        D.area(l)=-abs(D.area(l)); 
        if area2(l)>0, x=x(end:-1:1);y=y(end:-1:1); end;
      else
        if area2(l)<0, x=x(end:-1:1);y=y(end:-1:1); end; 
      end;
   
      % Here we try to reduce the number of points.
      
      xflag=0;
      if max(x)>180, % First, save original curve for later if we anticipate
        sx=x;sy=y;   % a 180-problem.
        xflag=1;
      end;
      
      % Look for points outside the lat/long boundaries, and then decimate them
      % by a factor of about 'decfac' (don't get rid of them completely because that
      % can sometimes cause problems when polygon edges cross curved map edges).
      
      tol=.2;   
     
      % Do y limits, then x so we can keep corner points.
      
      nn=(y>mtlim+tol) | (y<mblim-tol);
        % keep one extra point when crossing limits, also the beginning/end point.
      nn=logical(nn-min(1,([0;diff(nn)]>0)+([diff(nn);0]<0)));
      nn([1 end])=0;
        % decimate vigorously
      nn=nn & rem(1:length(nn),decfac)'~=0;
      x(nn)=[];y(nn)=[];
            
      if mrlim>mllim,  % no wraparound
          % sections of line outside lat/long limits
        nn=(x>mrlim+tol | x<mllim-tol) & y<mtlim & y>mblim;
       else            % wraparound case
        nn=(x>mrlim+tol & x<mllim-tol ) & y<mtlim & y>mblim;
      end;
      nn=logical(nn-min(1,([0;diff(nn)]>0)+([diff(nn);0]<0)));nn([1 end])=0;
      nn=nn & rem(1:length(nn),decfac)'~=0;
      x(nn)=[];y(nn)=[];
      
      % Move all points "near" to map boundaries.
      % I'm not sure about the wisdom of this - it might be better to clip
      % to the boundaries instead of moving. Hmmm. 
         
      y(y>mtlim+tol)=mtlim+tol;
      y(y<mblim-tol)=mblim-tol;
      if mrlim>mllim,   % Only clip long bdys if I can tell I'm on the right
                        % or left (i.e. not in wraparound case)
        x(x>mrlim+tol)=mrlim+tol;
        x(x<mllim-tol)=mllim-tol;
      end;   
      
    %% save patch data in struct
    
      D.n (l)        = length(x)+1;
      D.id(l)        = A(1);
      D.version(l)   = bitand(bitshift(A(3), -8),255);
      D.level(l)     = bitand(         A(3)     ,255);
      D.greenwich(l) = bitand(bitshift(A(3),-16),255);
      D.source(l)    = bitand(bitshift(A(3),-24),255);
      
                     D.k(l+1)   = D.k(l)+length(x)+1;
      D.lon(D.k(l)+1:D.k(l+1)-1)= x;
      D.lat(D.k(l)+1:D.k(l+1)-1)= y;
      D.lon(         D.k(l+1)  )= NaN;
      D.lat(         D.k(l+1)  )= NaN;
      
      % This is a little tricky...the filling algorithm expects data to be in the
      % range -180 to 180 deg long. However, there are some land parts that cut across
      % this divide so they appear at +190 but not -170. This causes problems later...
      % so as a kludge I replicate some of the problematic features at 190-360=-170.
      % Small islands are just duplicated, for the Eurasian landmass I just clip
      % off the eastern part.
      
      if xflag,
        l=l+1;D.area(l)=D.area(l-1); 
        if abs(D.area(l))>1e5,
          nn=find(sx>180);nn=[nn;nn(1)];
                         D.k(l+1)   =D.k(l)+length(nn)+1;
          D.lon(D.k(l)+1:D.k(l+1)-1)=sx(nn)-360;
          D.lat(D.k(l)+1:D.k(l+1)-1)=sy(nn);
        else   % repeat the island at the other edge.
                         D.k(l+1)   =D.k(l)+length(sx)+1;
          D.lon(D.k(l)+1:D.k(l+1)-1)=sx-360;
          D.lat(D.k(l)+1:D.k(l+1)-1)=sy;
        end;
        D.lon(D.k(l+1))=NaN;
        D.lat(D.k(l+1))=NaN;
      end;
      
    end;
    
    
    %%[A,cnt]=fread(fid,9,'int32');
    [A,cnt]=get_gheader(fid);
    
   end;
   
   fclose(fid);
   
%% clean up

   D.n_unique = nu;
   
   D.lon       ((D.k(l+1)+1):end)=[];  % get rid of unused part of data matrices
   D.lat       ((D.k(l+1)+1):end)=[];  % get rid of unused part of data matrices
   D.area      (    (l+1)   :end)=[];
   D.id        (    (l+1)   :end)=[];
   D.n         (    (l+1)   :end)=[];
   D.level     (    (l+1)   :end)=[];
   D.version   (    (l+1)   :end)=[];
   D.greenwich (    (l+1)   :end)=[];
   D.source    (    (l+1)   :end)=[];

   D.west      (    (l+1)   :end)=[];
   D.east      (    (l+1)   :end)=[];
   D.south     (    (l+1)   :end)=[];
   D.north     (    (l+1)   :end)=[];
 
   D.k         (    (l+2)   :end)=[];

%% GET_GHEADER
function [A,cnt]=get_gheader(fid);

% Reads the gshhs file header
% 
% A bit of code added because header format changed with version 1.3.
%
% 17/Sep/2008 - added material to handle latest GSHHS version.
%
% For version 1.1 this is the header ( 9*4 = 36 bytes long)
%
%int id;				/* Unique polygon id number, starting at 0 */
%int n;				/* Number of points in this polygon */
%int level;			/* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
%int west, east, south, north;	/* min/max extent in micro-degrees */
%int area;			/* Area of polygon in 1/10 km^2 */
%short int greenwich;		/* Greenwich is 1 if Greenwich is crossed */
%short int source;		/* 0 = CIA WDBII, 1 = WVS */
%
% For version 1.3 of GMT format was changed to this ( 10*4 = 40 bytes long)
%
%int id;				/* Unique polygon id number, starting at 0 */
%int n;				/* Number of points in this polygon */
%int level;			/* 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake */
%int west, east, south, north;	/* min/max extent in micro-degrees */
%int area;			/* Area of polygon in 1/10 km^2 */
%int version;			/* Polygon version, set to 3
%short int greenwich;		/* Greenwich is 1 if Greenwich is crossed */
%short int source;		/* 0 = CIA WDBII, 1 = WVS */
%
% For version 1.4, we have (8*4 = 32 bytes long)
%
%int id;				/* Unique polygon id number, starting at 0 */
%int n;				/* Number of points in this polygon */
%int flag;			/* level + version << 8 + greenwich << 16 + source << 24 
%int west, east, south, north;	/* min/max extent in micro-degrees */
%int area;			/* Area of polygon in 1/10 km^2 */
%
%Here, level, version, greenwich, and source are
%level:		1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake
%version:	Set to 4 for GSHHS version 1.4
%greenwich:	1 if Greenwich is crossed
%source:		0 = CIA WDBII, 1 = WVS
%

   % Now, in the calling code I have to use A(2),A(3),A(5-7), A(8), A(9) from original.
   
   [A,cnt]=fread(fid,8,'int32');
   
   if cnt<8,  % This gets triggered by the EOF
     return;
   end;
     
   ver=bitand(bitshift(A(3),-8),255);
   
   if ver==0,  % then its an old version
   
     % This works for version 1.2, but not 1.3.
   
     [A2,cnt2]=fread(fid,1,'int32');
     A=[A;A2];
   
     if (cnt+cnt2)==9 & A(9)==3,  % we have version 1.3, this would be one of 0,1,65535,65536 in v 1.2
       % Read one more byte
       A2=fread(fid,1,'int32');
       % This is the easiest way not to break existing code.
       A(9)=A2;   % one of 0,1,65536,65537
     end;
   %%fprintf('%d ',A(9));
   
   else  % a newest versions
   
    level     = bitand(A(3),255);
    greenwich = bitand(bitshift(A(3),-16),255);
    source    = bitand(bitshift(A(3),-24),255);
    A(3)      = level;
    A(9)      = greenwich*65536;
    
   end;

%% EOF