function D=gshhs_read(optn,file,varargin);
%GSHHS_READ   extract number of segments and number of points in a 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.
%
%
% D = gshhs_index('f','gshhs_f.b')
% D = gshhs_index('h','gshhs_h.b')
% D = gshhs_index('i','gshhs_i.b')
% D = gshhs_index('l','gshhs_l.b')
% D = gshhs_index('c','gshhs_c.b')
%
%    n_unique:   188603
%    n_points: 10340934
%
%    n_unique:   153535
%    n_points:  1906488
%
%    n_unique:    41518
%    n_points:   450663
%
%    n_unique:    10768
%    n_points:    90239
%
%    n_unique:     1866
%    n_points:    12967
%
%See also: M_COAST, GSHHS

%% load   
   
   fid   = fopen(file,'r','ieee-be');
   
   % Read the File header
    [A,cnt]=get_gheader(fid);
   
   nu = 0; % number of patches in file
   np = 0; % number of patches in file
   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.

    C=fread(fid,A(2)*2,'int32'); % Read all points in the current segment.
    
     np = np + A(2);

    %%[A,cnt]=fread(fid,9,'int32');
    [A,cnt]=get_gheader(fid);
    
   end;
   
   fclose(fid);
   
%% clean up

D.n_unique = nu;
D.n_points = np;
 
%% 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