function shapewrite(S, filename, fieldInfo)
%SHAPEWRITE Write a geographic data stucture to a shapefile.
%   SHAPEWRITE(S, FILENAME) writes a Version 2 geographic data stucture (a
%   geostruct2) to disk in shapefile format.  S must be a valid geostruct2
%   with specific restrictions on its attribute fields:
%
%   * Each attribute field value must be either a real, finite, scalar
%     double or a character string.
%
%   * The type of a given attribute must be consistent across all features.
%
%   FILENAME must be a character string specifying the output file name and
%   location.  If an extension is included, it must be '.shp' or '.SHP'.
%   SHAPEWRITE creates three output files,
%
%                        [BASENAME '.shp']
%                        [BASENAME '.shx']
%                        [BASENAME '.dbf']
%
%   where BASENAME is FILENAME without its extension.
%
%   If a given attribute is integer-valued for all features, then it is
%   written to the [BASENAME '.dbf'] file as an integer.  If an attribute
%   is a non-integer for any feature, then it is written as a fixed point
%   decimal value with six digits to the right of the decimal place.
%
%   Example
%   -------
%   % Derive a shapefile from concord_roads.shp in which roads of CLASS 5
%   % and greater are omitted.  Note the use of the 'Selector' option in
%   % shaperead, together with an anonymous function, to read only the main
%   % roads from the original shapefile.
%   shapeinfo('concord_roads')  % 609 features
%   S = shaperead('concord_roads', 'Selector', ...
%                 {@(roadclass) roadclass < 4, 'CLASS'});
%   shapewrite(S, 'main_concord_roads.shp')
%   shapeinfo('main_concord_roads')  % 107 features
%
%   See also SHAPEINFO, SHAPEREAD, UPDATEGEOSTRUCT.

%   Copyright 2003-2004 The MathWorks, Inc.
%   $Revision: 1.1.10.1 $  $Date: 2004/12/18 07:46:35 $

basename = validateFilename(filename);
validate(S)

if nargin == 2
    fieldInfo = dbffieldinfo(S);
end

[shapeType, boundingBox, index] = writeSHP(S,basename);
writeSHX(shapeType, boundingBox, index, basename);
if ~isempty(fieldInfo)
    dbfwrite(S, basename, fieldInfo)
end

%--------------------------------------------------------------------------

function basename = validateFilename(filename)

checkinput(filename, {'char'}, {'vector'}, mfilename, 'FILENAME', 2);
[pathstr, name, ext] = fileparts(filename);
if ~any(strcmpi(ext,{'','.shp'}))
    eid = sprintf('%s:%s:invalidExtension',getcomp,mfilename);
    msg = sprintf(...
        'Function %s expected filename\n  %s\nto have the extension: ''shp''.',...
        mfilename, filename);
    error(eid,msg)
end
basename = fullfile(pathstr,name);

%----------------------------------------------------------------------------

function validate(S)

% Make sure that S is a non-empty structure
if ~isstruct(S) || isempty(S)
    eid = sprintf('%s:%s:notAStruct',getcomp,mfilename);
    msg = sprintf('Geostruct S must be a non-empty structure array.');
    error(eid,msg);
end    
    
% Make sure there's a geometry field
if ~isfield(S, 'Geometry')
    eid = sprintf('%s:%s:noGeometry',getcomp,mfilename);
    msg = sprintf('Geostruct S must have a ''Geometry'' field.');
    error(eid,msg);
end

% Make sure the geometry field of the first feature is a string
if ~ischar(S(1).Geometry)
    eid = sprintf('%s:%s:nonStrGeometry',getcomp,mfilename);
    msg = sprintf('The ''Geometry'' field of geostruct S must contain a string.');
    error(eid,msg);
end

% Make sure the geometry of the first feature has a valid value
validGeometries = {'Point', 'MultiPoint', 'Line', 'PolyLine', 'Polygon'};
if ~any(strcmp(S(1).Geometry,validGeometries))
    eid = sprintf('%s:%s:invalidGeometryString',getcomp,mfilename);
    msg = sprintf('%s\n%s',...
        'The ''Geometry'' field of geostruct S must contain one of the following',...
        'strings: ''Point'', ''MultiPoint'', ''Line'', ''PolyLine'', or ''Polygon''.');
    error(eid,msg);
end

% Get the attribute field names
attributeNames = fields(S);
[t1,fIndex,t2] = setxor(attributeNames,{'Geometry','BoundingBox','X','Y','Lat','Lon'});
attributeNames = attributeNames(sort(fIndex));

% Determine the data classes for the attributes of the first feature
dataClasses = cell(numel(attributeNames),1);
for j = 1:numel(attributeNames)
    dataClasses{j} = class(S(1).(attributeNames{j}));
end

% Make sure that the rest of the features are consistent with the
% first, for numerical (N) and character (C) attributes.  Ignore
% the others, which will be filtered out by dbffieldinfo.  In the
% same loop, check for consistent geometry.
for k = 1:numel(S)
    if ~strcmp(S(k).Geometry,S(1).Geometry)
        eid = sprintf('%s:%s:inconsistentGeometry',getcomp,mfilename);
        msg = 'All features in geostruct S must have the same geometry.';
        error(eid,msg);
    end
    for j = 1:numel(attributeNames)
        v = S(k).(attributeNames{j});
        switch(dataClasses{j})
            case 'char'
                if ~ischar(v)
                    eid = sprintf('%s:%s:inconsistentAttrClassC',getcomp,mfilename);
                    msg = sprintf('Inconsistent data classes in attribute field: %s.',...
                                  attributeNames{j});
                    error(eid,msg);
                end
            case 'double'
                if ~isa(v,'double')
                    eid = sprintf('%s:%s:inconsistentAttrClassN',getcomp,mfilename);
                    msg = sprintf('Inconsistent data classes in attribute field: %s.',...
                                  attributeNames{j});
                    error(eid,msg);
                end
                if isempty(v)
                    eid = sprintf('%s:%s:emptyNumericalAttribute',getcomp,mfilename);
                    msg = sprintf('Numerical attributes of geostruct S must be non-empty.');
                    error(eid,msg);
                end
                if (numel(v) ~= 1)
                    eid = sprintf('%s:%s:nonScalarAttribute',getcomp,mfilename);
                    msg = sprintf('Numerical attributes of geostruct S must be scalar.');
                    error(eid,msg);
                end
                if ~isfinite(v) || ~isreal(v)
                    eid = sprintf('%s:%s:attributeNotFiniteReal',getcomp,mfilename);
                    msg = sprintf('Numerical attributes of geostruct S must be finite and real.');
                    error(eid,msg);
                end
        end
    end
end

% Check the coordinate field names
hasXY     = (isfield(S,'X') && isfield(S,'Y'));
hasLatLon = (isfield(S,'Lon') && isfield(S,'Lat'));
if ~hasXY && ~hasLatLon
    eid = sprintf('%s:%s:missingCoordinateFields',getcomp,mfilename);
    msg = sprintf('Geographic data structure S must have coordinate fields:\n%s',...
          'either ''X'' and ''Y'' or ''Lat'' and ''Lon''.');
    error(eid,msg)
end    
if (hasXY && hasLatLon)...
        || (hasXY     && (isfield(S,'Lat') || isfield(S,'Lon')))...
        || (hasLatLon && (isfield(S,'X')   || isfield(S,'Y')))
    eid = sprintf('%s:%s:redundantCoordinateFields',getcomp,mfilename);
    msg = sprintf('Geographic data structure S must have coordinate fields:\n%s',...
          'either ''X'' and ''Y'' or ''Lat'' and ''Lon'', but not both.');
    error(eid,msg)
end    

% Check the coordinate values
if hasXY
    cfield1 = 'X';
    cfield2 = 'Y';
else
    cfield1 = 'Lon';
    cfield2 = 'Lat';
end
for k = 1:numel(S)
    % Make sure they're doubles
    if ~isa(S(k).(cfield1),'double') || ~isa(S(k).(cfield2),'double')
        eid = sprintf('%s:%s:nonDoubleCoordinate',getcomp,mfilename);
        msg = sprintf('Coordinates in geographic data structure S must be doubles.');
        error(eid,msg)
    end
    
    % Make sure they match in length
    
    % Make sure they have matching NaNs, in any
    
    % Make sure they're real and finite
    
end

%----------------------------------------------------------------------------
function [shapeType, boundingBox, index] = writeSHP(S,basename)
% Write the main (SHP) file.

% Open the SHP file.
fid = fopen([basename '.shp'],'w','ieee-be');
if fid < 0
    eid = sprintf('%s:%s:failedToOpenSHP',getcomp,mfilename);
    msg = sprintf('Unable to open %s for writing.',[basename '.shp']);
    error(eid,msg);
end

% Get the shape type and a handle to function to write individual records.
shapeType = getShapeType(S(1).Geometry);
switch(S(1).Geometry)
    case 'Point'
        writefcn = @shpWritePoint;
    case 'MultiPoint'
        writefcn = @shpWriteMultiPoint;
    case {'Line','PolyLine','Polygon'}
        writefcn = @shpWritePoly;
end

% Write 100 bytes of zeros to reserve room for the file header.
fwrite(fid, uint8(zeros(1,100)), 'uint8');

% Write an SHP record for each element in S.
% Accumulate an index array and bounding box.
headerLengthInWords = 50;
fileLengthInWords = headerLengthInWords;
boundingBox = [Inf Inf -Inf -Inf];
index = zeros(2,numel(S));
for k = 1:numel(S)
    [fileLengthInWords, boundingBox, index(:,k)] = ...
        shpWriteRecord(fid, writefcn, shapeType, S(k), k, fileLengthInWords, boundingBox);
end

% Back up to the beginning of the file and write the header into the first 100 bytes.
fseek(fid, 0, 'bof');
shpWriteHeader(fid,shapeType,fileLengthInWords,boundingBox);

% Close the SHP file.
fclose(fid);

%----------------------------------------------------------------------------
function [fileLengthInWords, boundingBox, index] = shpWriteRecord(fid,...
      writefcn, shapeType, s, recordNumber, fileLengthInWords, boundingBox)
% Write an SHP record for each structure element s. Obtain coordinates from
% the 'X' and 'Y' fields, if present.  Otherwise use the 'Lon' and 'Lat'
% fields.  Update file length and bounding box.  Return index entry.

if isfield(s,'X') && isfield(s,'Y')
    xField = 'X';
    yField = 'Y';
elseif isfield(s,'Lon') && isfield(s,'Lat')
    xField = 'Lon';
    yField = 'Lat';
else
    eid = sprintf('%s:%s:missingCoordinateField',getcomp,mfilename);
    msg = 'Geographic data must have coordinate fields: either X and Y or Lat and Lon.';
    error(eid,msg)
end    

recordHeaderLengthInWords = 4;
offsetInWords = ftell(fid) / 2;
[contentLengthInWords, recordBoundingBox] ...
    = feval(writefcn, fid, recordNumber, shapeType, s.(xField), s.(yField));
fileLengthInWords = fileLengthInWords ...
    + recordHeaderLengthInWords + contentLengthInWords;
boundingBox(1:2) = min([boundingBox(1:2); recordBoundingBox(1:2)]);
boundingBox(3:4) = max([boundingBox(3:4); recordBoundingBox(3:4)]);
index = [offsetInWords; contentLengthInWords];

%----------------------------------------------------------------------------
function writeSHX(shapeType, boundingBox, index, basename)
% Write the SHX file.

% Open the SHX file.
fid = fopen([basename '.shx'],'w','ieee-be');
if fid < 0
    eid = sprintf('%s:%s:failedToOpenSHX');
    msg = sprintf('Unable to open %s for writing.',[basename '.shx']);
    error(eid,msg);
end

% Write the header.
headerLengthInWords = 50;
fileLengthInWords = headerLengthInWords + 4 * size(index,2);
shpWriteHeader(fid,shapeType,fileLengthInWords,boundingBox);

% Write the index records.
count = fwrite(fid,int32(index),'int32','ieee-be');

% Close the SHX file.
fclose(fid);

%----------------------------------------------------------------------------
function shapeType = getShapeType(geometry)

switch(geometry)
    case 'Point'
        shapeType = 1;
    case 'MultiPoint'
        shapeType = 8;
    case {'Line','PolyLine'}
        shapeType = 3;
    case {'Polygon'}
        shapeType = 5;
    otherwise
        shapeType = 0;
end

%----------------------------------------------------------------------------
function shpWriteHeader(fid, shapeType, fileLengthInWords, boundingBox)

fileCode = 9994;
version  = 1000;

bytes0thru27 = int32([fileCode 0 0 0 0 0 fileLengthInWords]);
count = fwrite(fid, bytes0thru27, 'int32', 'ieee-be');
% Check that count = 7

bytes28thru35 = int32([version shapeType]);
count = fwrite(fid, bytes28thru35, 'int32', 'ieee-le');
% Check that count = 2

bytes36thru99 = [boundingBox 0 0 0 0];
count = fwrite(fid, bytes36thru99, 'double', 'ieee-le');
% Check that count = 8

%----------------------------------------------------------------------------
function [contentLengthInWords, boundingBox] ...
    = shpWritePoint(fid,recordNumber,shapeType,X,Y)
% Write a Point record. Return the record content length measured in
% 16-bit words and the bounding box for the current record.

% Assign a degenerate bounding box for use in calculating the bounding box
% for the overall file.
boundingBox = [X Y X Y]; 

% The content length (the length of the record excluding its header) is
% fixed at 20 bytes (10 words) for Point shapes.
contentLengthInWords = 10;

% Write the record header
count = fwrite(fid, int32([recordNumber contentLengthInWords]), 'int32', 'ieee-be');

% Write the shape type.
count = fwrite(fid, int32(shapeType), 'int32', 'ieee-le');

% Write the point coordinates.
count = fwrite(fid, [X Y], 'double', 'ieee-le');

%----------------------------------------------------------------------------
function [contentLengthInWords, boundingBox] ...
    = shpWriteMultiPoint(fid,recordNumber,shapeType,X,Y)
% Write a MultiPoint record. Return the record content length measured in
% 16-bit words and the bounding box for the current record.

% Calculate the 2-by-numPoints points array and the bounding box.
points = [X(:)'; Y(:)'];
boundingBox = [min(X) min(Y) max(X) max(Y)]; 

% Calculate the content length (the length of the record excluding its header).
contentLengthInWords = (40 + 8 * numel(points)) / 2;

% Write the record header
count = fwrite(fid, int32([recordNumber contentLengthInWords]), 'int32', 'ieee-be');

% Write the shape type.
count = fwrite(fid, int32(shapeType), 'int32', 'ieee-le');

% Write the bounding box.
count = fwrite(fid, boundingBox, 'double', 'ieee-le');

% Write the number of points.
count = fwrite(fid, size(points,2), 'int32', 'ieee-le');

% Write the point coordinates.
count = fwrite(fid, points, 'double', 'ieee-le');

%----------------------------------------------------------------------------
function [contentLengthInWords, boundingBox] ...
    = shpWritePoly(fid,recordNumber,shapeType,X,Y)
% Write a PolyLine or Polygon record. Return the record content length measured in
% 16-bit words and the bounding box for the current record.

% Calculate the 1-by-numParts parts index, the 2-by-numPoints points array,
% and the bounding box.
[parts, points] = constructPartsIndexAndPointsArray(X,Y);
boundingBox = [min(X) min(Y) max(X) max(Y)]; 

% Calculate the content length (the length of the record excluding its header).
contentLengthInWords = (44 + 4 * numel(parts) + 8 * numel(points)) / 2;

% Write the record header
count = fwrite(fid, int32([recordNumber contentLengthInWords]), 'int32', 'ieee-be');

% Write the shape type.
count = fwrite(fid, int32(shapeType), 'int32', 'ieee-le');

% Write the bounding box.
count = fwrite(fid, boundingBox, 'double', 'ieee-le');

% Write the part and point counts, and the parts index.
count = fwrite(fid, int32([numel(parts) size(points,2) parts]), 'int32','ieee-le');

% Write the point coordinates.
count = fwrite(fid, points, 'double', 'ieee-le');

%----------------------------------------------------------------------------
function [parts, points] = constructPartsIndexAndPointsArray(X,Y)
% Calculate the zero-based index array PARTS giving the offset to the where
% the start of each non-NaN sequence in vector X will be after all the NaNs
% were removed from X.  Organize the non-NaN elements of X and Y into a
% 2-by-numPoints array POINTS with X-coordinates in row 1 and Y-coordinates
% in row 2.

% Check that isequal(isnan(X),isnan(Y)). If not, then throw an error here
% and catch it above so we have a chance to close open files.  Also require
% that X and Y each have at least one non-NaN element.

% Determine where the NaNs are and add a (virtual) terminating NaN.
nanplaces = [find(isnan(X(:))); 1+numel(X)];

% Compute the parts index. Note that in both the initializion of parts and
% in subsequent updates, we calculate the offset to the next part before
% determining if that part exists.  This is fine, except that at the end we
% end up with parts(end) containing the offset of a non-existent part, so
% we must remove it.
parts = 0;
indexOfLastNaN = 0;
for k = 1:numel(nanplaces)
    partLength = nanplaces(k) - (indexOfLastNaN + 1);
    if partLength > 0
        parts(end+1) = parts(end) + partLength;
    end
    % If partLength is zero, then the current NaN was preceded
    % immediately by another NaN, and should hence be ignored.
    indexOfLastNaN = nanplaces(k);
end
parts(end) = [];

% Calculate the points array:
%    [X(1) X(2) ...;
%     Y(1) Y(2) ...]
% after removing NaNs.
x = X(:)'; x(isnan(x)) = [];
y = Y(:)'; y(isnan(y)) = [];
points = [x; y];

%----------------------------------------------------------------------------
function dbfwrite(S, basename, fieldInfo)
% Write the DBF file.

% Open the DBF file.
fid = fopen([basename '.dbf'],'w','ieee-le');
if fid < 0
    eid = sprintf('%s:%s:failedToOpenDBF');
    msg = sprintf('Unable to open %s for writing.',[basename '.dbf']);
    error(eid,msg);
end

dbfWriteTableFileHeader(fid, numel(S), fieldInfo);
for k = 1:numel(S)
    dbfWriteTableRecord(fid, S(k), fieldInfo);
end
dbfMarkEOF(fid);

% Close the DBF file.
fclose(fid);

%----------------------------------------------------------------------------
function dbfWriteTableFileHeader(fid, numberOfRecords, fieldInfo)

version = 3;

dateVector = datevec(date);
year  = dateVector(1);
month = dateVector(2);
day   = dateVector(3);

nFields = numel(fieldInfo);
headerLength = 32 + 32 * nFields + 1;
recordLength = 1 + sum([fieldInfo.FieldLength]);

% Bytes 0-3: dBASE version and today's date
fwrite(fid, [version year-1900 month day], 'uint8');

% Bytes 4-7: Number of records in the table
fwrite(fid, numberOfRecords, 'uint32');

% Bytes 8-9: Number of bytes in the header
fwrite(fid, headerLength, 'uint16');

% Bytes 10-11: Number of bytes per record
fwrite(fid, recordLength, 'uint16');

% Bytes 12-31: Reserved (zero-fill)
fwrite(fid, zeros(1,20), 'uint8');

% Bytes 32-n: Table field descriptors (n = 32 + 32 * nFields)
for k = 1:nFields
    dbfWriteTableFieldDescriptor(fid, fieldInfo(k));
end

% Byte n+1: Field terminator
fwrite(fid, hex2dec('0D') ,'uint8');

%----------------------------------------------------------------------------
function dbfWriteTableFieldDescriptor(fid, fieldInfo)

% Bytes 0-10: Field name in ASCII
%   (truncate or zero-fill to precisely 11 bytes)
fwrite(fid, truncateOrFill(fieldInfo.FieldName,11), 'uchar');

% Byte 11: Field type in ASCII ('C' or 'N')
fwrite(fid, fieldInfo.FieldType, 'uchar');

% Bytes 12-15: Field data address (zero-fill)
fwrite(fid, zeros(1,4), 'uint8');

% Byte 16: Field length
fwrite(fid, fieldInfo.FieldLength, 'uint8');

% Byte 17: Field decimal count
%   (number of digits to the right of the decimal place)
fwrite(fid, fieldInfo.FieldDecimalCount, 'uint8');

% Bytes 18-31: Reserved or zero
%   ("work area ID, SET FIELDS flag", ".MDX field flag")
fwrite(fid, zeros(1,14), 'uint8');

%----------------------------------------------------------------------------
function str = truncateOrFill(str,len)
% Ensure a string containing precisely LEN characters,
% padded with zeros if necessary.

str(len+1:end) = [];
str(end+1:len) = 0;

%----------------------------------------------------------------------------
function dbfWriteTableRecord(fid, s, fieldInfo)
% Write a single record, including the preceding SPACE character.

% Precede the record with a SPACE (0x20)
spaceChar = hex2dec('20');
fwrite(fid, spaceChar, 'uint8');

% Write each field
for k = 1:numel(fieldInfo)
    value = s.(fieldInfo(k).AttributeName);
    str = sprintf(fieldInfo(k).Format,value);
    fwrite(fid,str,'uchar');
end

%----------------------------------------------------------------------------
function dbfMarkEOF(fid)
% Mark end of dBASE file.

eofMarker = hex2dec('1A');
fwrite(fid, eofMarker, 'uint8');