function [shape,SHAPE,data]=readshp(FName)
% [shape,SHAPE]=readshp(FName) ---- use no extension
% Reads ESRI shape file, specified in
% 'http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf'
% A shapefile consists of a mainfile an strmatchi file and a dbase file
% e.g. basename.shp, basename.shx, basename.dbf, basename the same, fnames 8.3 convention
%
% shape is a struct array of shapes which includes the shape data extracted
% from the dbf file
% SHAPE is overall info over all shapes
% data is the contents of the dbf file (not needed, also per record in shape-struct)
%
% dbfread must be in the search path 
% 
% See also WRITESHP, READEXP, WRITEEXP, PLOTSHP
%
% TO 090729


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later

%%
fprintf('\n\nreadshp --- %s\n',datestr(now));

if strcmpi(FName(end-3:end),'.shp'), FName=FName(1:end-4); end  % cut off .shp if necessary

headerlen=100;  % for shp and shx files

%% Legal shape types
shapetypes={
    0, 'Null Shape'
    1, 'Point'
    3, 'PolyLine'
    5, 'Polygon'
    8, 'MultiPoint'
    11,'PointZ'
    13,'PolyLineZ'
    15,'PolygonZ'
    18,'MultiPointZ'
    21,'PontM'
    23,'PolyLineM'
    25,'PolygonM'
    28,'MultiPointM'
    31,'MultiPatch'};

%% PartTypes defined for MultiPatch
parttypes={
    0,'TriangleStrip'
    1,'TriangleFan'
    2,'OuterRing'
    3,'InnerRing'
    4,'FirstRing'
    5,'Ring'};

%% Skip the header because it is idential to that of the .shp file, where
%% we'll read it

shxB =fopen([FName,'.shx'],'r','b');  % big endian bite ordering (Sun or Motorola)
shpB =fopen([FName,'.shp'],'r','b');  % big endian bite ordering (Sun or Motorola)
shpL =fopen([FName,'.shp'],'r','l');  % liggle endian (Intel or PC)

if shxB<1, error('Can''t open file ''%s.shx''\n',FName); end
if shpB<1, error('Can''t open file ''%s.shp''\n',FName); end
if shpL<1, error('Can''t open file ''%s.shp''\n',FName); end

shxHdr=fread(shxB,100,'int8');
shpHdr=fread(shpB,100,'int8');
if ~all(shxHdr([1:25 29:100])==shpHdr([1:25 29:100])), error('Header in .shx and .shp file differ, must be the same\n'); end

%% Header shx file
fseek(shxB, 0,'bof'); fprintf('FileCode  = %d\n',fread(shxB,1,'uint'));
fseek(shxB,24,'bof'); fprintf('FileLenShx= %d\n',fread(shxB,1, 'int'));  % number of bytes
fseek(shxB,28,'bof'); fprintf('Version   = %d\n',fread(shxB,1, 'int'));

%% Header shp file
fseek(shpB, 0,'bof'); fprintf('FileCode  = %d\n',fread(shpB,1,'uint'));
fseek(shpB,24,'bof'); fprintf('FileLenSp = %d\n',fread(shpB,1, 'int')); % in 16 bit words
fseek(shpB,28,'bof'); fprintf('Version   = %d\n',fread(shpB,1, 'int'));

fclose(shpB); % no longer needed

%% Get number of shape records from strmatchi file

SHAPE.name=FName;

fseek(shxB,0,'eof');
SHAPE.NShapes=(ftell(shxB)-headerlen)/8;

%% read bounding box
fseek(shpL,36,'bof');
    SHAPE.Xmin=fread(shpL,1,'double');
    SHAPE.Ymin=fread(shpL,1,'double');
    SHAPE.Xmax=fread(shpL,1,'double');
    SHAPE.Ymax=fread(shpL,1,'double');
    SHAPE.Zmin=fread(shpL,1,'double');
    SHAPE.Zmax=fread(shpL,1,'double');
    SHAPE.Mmin=fread(shpL,1,'double');
    SHAPE.Mmax=fread(shpL,1,'double');

%% Reading data

fseek(shxB,headerlen,'bof');  % Position pointer in shx strmatchi file
shape(SHAPE.NShapes).Type=0; % allocate
for iShape=1:SHAPE.NShapes
    % reading offset of record in shp file from record in shx file, skip
    % rest of header in skx file
    shpOffset =fread(shxB,1,'int')*2;    fread(shxB,1,'int');
   
    % position pointer in shp file read using big endian read offset in shp
    % file and skip rest of record header
    fseek(shpL,shpOffset,'bof');    fread(shpL,2,'int');
    
    % reading record contents, using shpL Little endian pointer
    shape(iShape).Type=fread(shpL,1,'int');
    j=find(shape(iShape).Type==vertcat(shapetypes{:,1}));
    if isempty(j), error('Illegal fileshapetype %d in file ''%s.shp''\n,',FName);
    else
       shape(iShape).TypeName=shapetypes{j,2};
    end
  
   switch shape(iShape).Type
       case 0  % Point, skip, only type was read
       case 1  % Point
          shape(iShape).Points=fread(shpL,[1,2],'double');
       case {3,5}  % 'PolyLine' 'Polygon'
          shape(iShape).Box      =fread(shpL,4,'double');
          shape(iShape).NumParts =fread(shpL,1,'int');
          shape(iShape).NumPoints=fread(shpL,1,'int');
          shape(iShape).Parts    =fread(shpL,shape(iShape).NumParts,'int');
          shape(iShape).Points   =fread(shpL,[2,shape(iShape).NumPoints],'double')';
       case 8, % 'MultiPoint'
          shape(iShape).Box         =fread(shpL,4,'double');
          shape(iShape).NumParts    =fread(shpL,1,'int');
          shape(iShape).NumPoints   =fread(shpL,1,'int');
          shape(iShape).Points      =fread(shpL,[2,shape(iShape).NumPoints],'double')';
          shape(iShape).MRange      =fread(shpL,[1,2],'double');
          shape(iShape).MArray      =fread(shpL,shape(iShape).NumPoints,'double');
       case 11, % 'PointZ'  [x y z Measure]
          shape(iShape).Point     =fread(shpL,[4,shape(iShape).NumPoints],'double')';         
       case {13,15} % 'PolyLineZ' 'PolygonZ'
          shape(iShape).Box      =fread(shpL,4,'double');
          shape(iShape).NumParts =fread(shpL,1,'int');
          shape(iShape).NumPoints=fread(shpL,1,'int');
          shape(iShape).Parts    =fread(shpL,shape(iShape).NumParts,'int')/2+1;
          shape(iShape).Points   =fread(shpL,[2,shape(iShape).NumPoints],'double')';
          shape(iShape).ZRange    =fread(shpL,[1,2],'double');
          shape(iShape).ZArray    =fread(shpL,shape(iShape).NumPoints,'double');
          shape(iShape).MRange    =fread(shpL,[1,2],'double');
          shape(iShape).MArray    =fread(shpL,shape(iShape).NumPoints,'double');
       case 18, % 'MultiPointZ'
          shape(iShape).Box      =fread(shpL,4,'double');
          shape(iShape).NumPoints=fread(shpL,1,'int');
          shape(iShape).Points   =fread(shpL,[2,shape(iShape).NumPoints],'double')';
          shape(iShape).ZRange    =fread(shpL,[1,2],'double');
          shape(iShape).ZArray    =fread(shpL,shape(iShape).NumPoints,'double');
          shape(iShape).MRange    =fread(shpL,[1,2],'double');
          shape(iShape).MArray    =fread(shpL,shape(iShape).NumPoints,'double');
       case 21, % 'PontM'          
          shape(iShape).Points   =fread(shpL,[1,2],'double')';
          shape(iShape).MArray   =fread(shpL,1,'double')';
       case {23,25} % 'PolyLineM' and 'PolygonM'
          shape(iShape).Box      =fread(shpL,4,'double');
          shape(iShape).NumParts =fread(shpL,1,'int');
          shape(iShape).NumPoints=fread(shpL,1,'int');
          shape(iShape).Parts    =fread(shpL,shape(iShape).NumParts,'int')/2+1;
          shape(iShape).Points   =fread(shpL,[2,shape(iShape).NumPoints],'double')';
          shape(iShape).MRange      =fread(shpL,[1,2],'double');
          shape(iShape).MArray      =fread(shpL,shape(iShape).NumPoints,'double');
       case 28, % 'MultiPointM'
          shape(iShape).Box      =fread(shpL,4,'double');
          shape(iShape).NumPoints=fread(shpL,1,'int');
          shape(iShape).Points    =fread(shpL,[2,shape(iShape).NumPoints],'double')';
          shape(iShape).ZRange    =fread(shpL,[1,2],'double');
          shape(iShape).ZArray    =fread(shpL,shape(iShape).NumPoints,'double');
          shape(iShape).MRange    =fread(shpL,[1,2],'double');
          shape(iShape).MArray    =fread(shpL,shape(iShape).NumPoints,'double');
       case 31, % 'MultiPatch'};
          shape(iShape).Box      =fread(shpL,4,'double');
          shape(iShape).NumParts =fread(shpL,1,'int');
          shape(iShape).NumPoints=fread(shpL,1,'int');
          shape(iShape).Parts    =fread(shpL,shape(iShape).NumParts,'int')/2+1;
          shape(iShape).PartTypes=fread(shpL,shape(iShape).NumParts,'int');
            for j=1:size(parttypes,1)
                if ~any(shape(iShape).PatTypes==parttypes(j,1))
                    error('Illegal parttypes found [outsed 0..5] in Multipatch, shape type 31\n');
                end
            end
          shape(iShape).Points   =fread(shpL,[2,shape(iShape).NumPoints],'double')';
          shape(iShape).ZRange    =fread(shpL,[1,2],'double');
          shape(iShape).ZArray    =fread(shpL,shape(iShape).NumPoints,'double');
          shape(iShape).MRange    =fread(shpL,[1,2],'double');
          shape(iShape).MArray    =fread(shpL,shape(iShape).NumPoints,'double');
       otherwise
           fprintf('Non implemented shape file %d skipped\n',shape.type)
   end
end
fclose(shpL);

%plot shapes see plotshp
% for iShape=1:length(shape)
%     for iParts=1:length(shape(iShape).NumParts)
%         first=shape(iShape).Parts(iParts)+1;
%         if iParts<shape(iShape).NumParts
%             last=shape(iShape).Parts(iParts+1);
%         else
%             last=shape(iShape).NumPoints;
%         end
%         range=first:last;
%         box=shape(iShape).Box;
%         plot(box([1 3 3 1 1]),box([2 2 4 4 2]),'r'); hold on
%         plot(shape(iShape).Points(range,1),shape(iShape).Points(range,2));
%         hold on
%     end
% end

%% DBF data
data=dbfread(FName)'; %PK Transformatie extra toegepast

for iShape=1:SHAPE.NShapes
    for iField=1:length(data)
                shape(iShape).(data(iField).fieldname)=data(iField).values(iShape);
    end
end
